Changeset 22825

Timestamp:

06/04/18 17:32:11 (7 years ago)

Author:

seroussi

Message:

CHG: finished exchange of ice base in the initializatio

File:

• issm/trunk-jpl/src/c/cores/transient_core.cpp (modified) (3 diffs)

issm/trunk-jpl/src/c/cores/transient_core.cpp

@@ -77 +77 @@
         #endif
 
-        if(isoceancoupling){ {{{
+        if(isoceancoupling){ /*{{{*/
+                #ifndef _HAVE_ADOLC_
                 if(VerboseSolution()) _printf0_("   ocean coupling: initialization \n");
                 int my_rank;
-                int oceangridnxsize,oceangridnysize;
-                IssmDouble *oceangridx;
-                IssmDouble *oceangridy;
-                IssmDouble *icebase;
-                IssmDouble coupling_time,oceantime;
                 ISSM_MPI_Comm tomitgcmcomm;
                 ISSM_MPI_Status status;
 
                 my_rank=IssmComm::GetRank();
-                femmodel->parameters->FindParam(&coupling_time,TimesteppingCouplingTimeEnum);
                 GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(ToMITgcmCommEnum));
                 if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator");
                 tomitgcmcomm=parcom->GetParameterValue();
+
                 if(my_rank==0){
+                        int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean;
+                        IssmDouble  oceantime,coupling_time;
+                        IssmDouble *oceangridx;
+                        IssmDouble *oceangridy;
+                        IssmDouble *base_oceangrid;
+                        IssmDouble* x_ice = NULL;
+                        IssmDouble* y_ice = NULL;
+                        IssmDouble* lat_ice = NULL;
+                        IssmDouble* lon_ice = NULL;
+                        IssmDouble* z_ice = NULL;
+                        IssmDouble* icebase= 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);
                         ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm);
                         ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status);
@@ -100 +113 @@
                         femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
                         femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
-                        oceangridx = xNew<IssmDouble>(oceangridnxsize*oceangridnysize);
-                        ISSM_MPI_Recv(oceangridx,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
-                        oceangridy = xNew<IssmDouble>(oceangridnxsize*oceangridnysize);
-                        ISSM_MPI_Recv(oceangridy,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
-                        femmodel->parameters->SetParam(oceangridx,oceangridnxsize*oceangridnysize,OceanGridXEnum);
-                        femmodel->parameters->SetParam(oceangridy,oceangridnxsize*oceangridnysize,OceanGridYEnum);
+                        ngrids_ocean=oceangridnxsize*oceangridnysize;
+                        oceangridx = xNew<IssmDouble>(ngrids_ocean);
+                        ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
+                        oceangridy = xNew<IssmDouble>(ngrids_ocean);
+                        ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
+                        femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);
+                        femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum);
 
                         /*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);
-                        icebase = xNew<IssmDouble>(oceangridnxsize*oceangridnysize);
-                        for(int i=0;i<oceangridnxsize;i++){
-                                for(int j=0;j<oceangridnysize;j++){
-                                        icebase[i*oceangridnysize+j]=2*oceangridx[i*oceangridnysize+j];
-                                }
-                        }
-                        ISSM_MPI_Send(icebase,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
+
+                        /*Interpolate ice base onto ocean grid*/
+                        femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
+                        BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
+                        femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
+                        GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum);
+                        InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
+                                                        icebase,ngrids_ice,1,
+                                                        oceangridx,oceangridy,ngrids_ocean,NULL);   
+
+                        ISSM_MPI_Send(base_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
+
+                        /*Delete*/
+                        xDelete<int>(index_ice);
+                        xDelete<int>(index_ocean);
+                        xDelete<IssmDouble>(lat_ice);
+                        xDelete<IssmDouble>(lon_ice);
+                        xDelete<IssmDouble>(x_ice);
+                        xDelete<IssmDouble>(y_ice);
+                        xDelete<IssmDouble>(z_ice);
                         xDelete<IssmDouble>(icebase);
+                        xDelete<IssmDouble>(base_oceangrid);
                         xDelete<IssmDouble>(oceangridx);
                         xDelete<IssmDouble>(oceangridy);
                 }
+                #else
+                _error_("not supported");
+                #endif
         }
-        }}}
+        /*}}}*/
 
                 IssmDouble  output_value;
@@ -219 +250 @@
                                 if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
                                 ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
-                                for(int i=0;i<ngrids_ice;i++) base_oceangrid[i]=0;
+                                //for(int i=0;i<ngrids_ice;i++) base_oceangrid[i]=0;
                                 ISSM_MPI_Send(base_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);