Changeset 22798


Ignore:
Timestamp:
05/21/18 13:26:18 (7 years ago)
Author:
seroussi
Message:

CHG: fixed ocean heat flux used in ISSM for ocean coupling

Location:
issm/trunk-jpl/src/c
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r22784 r22798  
    170170                        if(VerboseSolution()) _printf0_("   ocean coupling: exchanging information\n");
    171171                        int my_rank;
    172                         int oceangridnxsize,oceangridnysize,ngrids_ocean;
    173                         IssmDouble *oceanmelt;
    174                         IssmDouble *icebase;
    175                         IssmDouble *oceangridx;
    176                         IssmDouble *oceangridy;
    177                         IssmDouble oceantime;
    178172                        ISSM_MPI_Comm tomitgcmcomm;
    179173                        ISSM_MPI_Status status;
    … …  
    184178                        tomitgcmcomm=parcom->GetParameterValue();
    185179                        if(my_rank==0){
    186                                 ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
    187                                 ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
    188                                 if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
    189 
    190                                 /*Recover inputs needed*/
    191                                 femmodel->parameters->FindParam(&oceangridnxsize,OceanGridNxEnum);
    192                                 femmodel->parameters->FindParam(&oceangridnysize,OceanGridNyEnum);
    193                                 femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
    194                                 femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
    195 
    196                                 /*Exchange information*/
    197                                 oceanmelt = xNew<IssmDouble>(ngrids_ocean);
    198                                 ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
    199                                 icebase = xNew<IssmDouble>(ngrids_ocean);
    200                                 for(int i=0;i<oceangridnxsize;i++){
    201                                         for(int j=0;j<oceangridnysize;j++){
    202                                                 icebase[i*oceangridnysize+j]=9999.;
    203                                         }
    204                                 }
    205                                 ISSM_MPI_Send(icebase,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
    206 
    207                                 /*Interpolate melt onto mesh*/
     180                                int ngrids_ocean;
     181                                IssmDouble oceantime;
     182                                IssmDouble rho_ice;
     183                                IssmDouble *oceanmelt;
     184                                IssmDouble *base_oceangrid;
     185                                IssmDouble *oceangridx;
     186                                IssmDouble *oceangridy;
    208187                                IssmDouble* x_ice = NULL;
    209188                                IssmDouble* y_ice = NULL;
    … …  
    211190                                IssmDouble* lon_ice = NULL;
    212191                                IssmDouble* z_ice = NULL;
     192                                IssmDouble* icebase= NULL;
    213193                                IssmDouble* melt_mesh = NULL;
    214194                                int*        index_ice= NULL;
    … …  
    216196                                int         nels_ocean;
    217197                                int         ngrids_ice=femmodel->vertices->NumberOfVertices();
     198                                int         nels_ice=femmodel->elements->NumberOfElements();
    218199                                lat_ice= xNew<IssmDouble>(ngrids_ice);
    219200                                lon_ice= xNew<IssmDouble>(ngrids_ice);
     201                                icebase= xNew<IssmDouble>(ngrids_ice);
    220202                                melt_mesh= xNew<IssmDouble>(ngrids_ice);
     203                               
     204                                /*Recover mesh and inputs needed*/
     205                                femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
    221206                                femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
     207                                femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
     208                                femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
    222209                                BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
     210                                base_oceangrid= xNew<IssmDouble>(ngrids_ocean);
     211                                oceanmelt = xNew<IssmDouble>(ngrids_ocean);
     212
    223213                                femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
     214
     215                                /*Interpolate ice base onto ocean grid*/
     216//                              InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
     217                                InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
     218                                                        icebase,ngrids_ice,1,
     219                                                        oceangridx,oceangridy,ngrids_ocean,NULL);   
     220                               
     221                                /*Send and receive data*/
     222                                ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
     223                                ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
     224                                if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
     225                                ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
     226                                ISSM_MPI_Send(base_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
     227
     228                                /*Interp melt onto ice grid*/
    224229                                InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
    225230                                                        oceanmelt,ngrids_ocean,1,
    226231                                                        lon_ice,lat_ice,ngrids_ice,NULL);   
     232
     233                                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
    227234                                InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
     235
     236                                /*Delete*/
    228237                                xDelete<IssmDouble>(lat_ice);
    229238                                xDelete<IssmDouble>(lon_ice);
    … …  
    237246                                xDelete<IssmDouble>(oceanmelt);
    238247                                xDelete<IssmDouble>(icebase);
     248                                xDelete<IssmDouble>(base_oceangrid);
    239249                        }
    240250                        #else

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r22749 r22798  
    105105                parameters->AddObject(iomodel->CopyConstantObject("md.transient.iscoupler",TransientIscouplerEnum));
    106106                parameters->AddObject(iomodel->CopyConstantObject("md.transient.isoceancoupling",TransientIsoceancouplingEnum));
     107                parameters->AddObject(iomodel->CopyConstantObject("md.materials.rho_ice",MaterialsRhoIceEnum));
    107108                parameters->AddObject(iomodel->CopyConstantObject("md.transient.amr_frequency",TransientAmrFrequencyEnum));
    108109                parameters->AddObject(iomodel->CopyConstantObject("md.gia.cross_section_shape",GiaCrossSectionShapeEnum));

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r22787 r22798  
    200200        MasstransportRequestedOutputsEnum,
    201201        MasstransportStabilizationEnum,
     202        MaterialsRhoIceEnum,
    202203        MeltingOffsetEnum,
    203204        MeshAverageVertexConnectivityEnum,
    … …  
    834835        MaterialsRheologyLawEnum,
    835836        MaterialsRhoFreshwaterEnum,
    836         MaterialsRhoIceEnum,
    837837        MaterialsRhoSeawaterEnum,
    838838        MaterialsTemperateiceconductivityEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r22787 r22798  
    208208                case MasstransportRequestedOutputsEnum : return "MasstransportRequestedOutputs";
    209209                case MasstransportStabilizationEnum : return "MasstransportStabilization";
     210                case MaterialsRhoIceEnum : return "MaterialsRhoIce";
    210211                case MeltingOffsetEnum : return "MeltingOffset";
    211212                case MeshAverageVertexConnectivityEnum : return "MeshAverageVertexConnectivity";
    … …  
    838839                case MaterialsRheologyLawEnum : return "MaterialsRheologyLaw";
    839840                case MaterialsRhoFreshwaterEnum : return "MaterialsRhoFreshwater";
    840                 case MaterialsRhoIceEnum : return "MaterialsRhoIce";
    841841                case MaterialsRhoSeawaterEnum : return "MaterialsRhoSeawater";
    842842                case MaterialsTemperateiceconductivityEnum : return "MaterialsTemperateiceconductivity";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r22787 r22798  
    211211              else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum;
    212212              else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum;
     213              else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    213214              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    214215              else if (strcmp(name,"MeshAverageVertexConnectivity")==0) return MeshAverageVertexConnectivityEnum;
    … …  
    259260              else if (strcmp(name,"SealevelriseTidalLoveK")==0) return SealevelriseTidalLoveKEnum;
    260261              else if (strcmp(name,"SealevelriseTransitions")==0) return SealevelriseTransitionsEnum;
    261               else if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
     265              if (strcmp(name,"SealevelriseUElastic")==0) return SealevelriseUElasticEnum;
     266              else if (strcmp(name,"SettingsIoGather")==0) return SettingsIoGatherEnum;
    266267              else if (strcmp(name,"SettingsOutputFrequency")==0) return SettingsOutputFrequencyEnum;
    267268              else if (strcmp(name,"SettingsRecordingFrequency")==0) return SettingsRecordingFrequencyEnum;
    … …  
    382383              else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum;
    383384              else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
    384               else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
    385385         else stage=4;
    386386   }
    387387   if(stage==4){
    388               if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
     388              if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
     389              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
    389390              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    390391              else if (strcmp(name,"Base")==0) return BaseEnum;
    … …  
    505506              else if (strcmp(name,"Sealevel")==0) return SealevelEnum;
    506507              else if (strcmp(name,"SealevelriseDeltathickness")==0) return SealevelriseDeltathicknessEnum;
    507               else if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
    508508         else stage=5;
    509509   }
    510510   if(stage==5){
    511               if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
     511              if (strcmp(name,"SedimentHead")==0) return SedimentHeadEnum;
     512              else if (strcmp(name,"SedimentHeadOld")==0) return SedimentHeadOldEnum;
    512513              else if (strcmp(name,"SedimentHeadResidual")==0) return SedimentHeadResidualEnum;
    513514              else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum;
    … …  
    628629              else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
    629630              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
    630               else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
    631631         else stage=6;
    632632   }
    633633   if(stage==6){
    634               if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
     634              if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
     635              else if (strcmp(name,"WeightsSurfaceObservation")==0) return WeightsSurfaceObservationEnum;
    635636              else if (strcmp(name,"InputsEND")==0) return InputsENDEnum;
    636637              else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
    … …  
    751752              else if (strcmp(name,"GaussTetra")==0) return GaussTetraEnum;
    752753              else if (strcmp(name,"GaussTria")==0) return GaussTriaEnum;
    753               else if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
     757              if (strcmp(name,"GenericOption")==0) return GenericOptionEnum;
     758              else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
    758759              else if (strcmp(name,"GiadWdt")==0) return GiadWdtEnum;
    759760              else if (strcmp(name,"GiaIvinsAnalysis")==0) return GiaIvinsAnalysisEnum;
    … …  
    856857              else if (strcmp(name,"MaterialsRheologyLaw")==0) return MaterialsRheologyLawEnum;
    857858              else if (strcmp(name,"MaterialsRhoFreshwater")==0) return MaterialsRhoFreshwaterEnum;
    858               else if (strcmp(name,"MaterialsRhoIce")==0) return MaterialsRhoIceEnum;
    859859              else if (strcmp(name,"MaterialsRhoSeawater")==0) return MaterialsRhoSeawaterEnum;
    860860              else if (strcmp(name,"MaterialsTemperateiceconductivity")==0) return MaterialsTemperateiceconductivityEnum;
Note: See TracChangeset for help on using the changeset viewer.