Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24915) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 24916) @@ -4684,5 +4684,5 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop) { /*{{{*/ /*serialized vectors:*/ @@ -4706,7 +4706,9 @@ /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ for(int i=0;iSize();i++){ + IssmDouble area; Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - oceanarea_cpu += element->OceanArea(); - eartharea_cpu += element->GetAreaSpherical(); + area=element->GetAreaSpherical(); + eartharea_cpu += area; + if (element->IsOceanInElement()) oceanarea_cpu += area; } ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); @@ -4737,15 +4739,14 @@ /*Assign output pointers:*/ + *peartharea=eartharea; + *poceanarea=oceanarea; *peustatic=eustatic; } /*}}}*/ -void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ +void FemModel::SealevelriseNonEustatic(Vector* pRSLgo, Vector* pRSLg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, bool verboseconvolution,int loop){/*{{{*/ /*serialized vectors:*/ IssmDouble* RSLg_old=NULL; - - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; IssmDouble* RSLgo = NULL; @@ -4762,13 +4763,4 @@ RSLg_old=pRSLg_old->ToMPISerial(); - /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - eartharea_cpu += element->GetAreaSpherical(); - } - - ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - /*Call the sea level rise core: */ for(int i=0;iSize();i++){ @@ -4786,10 +4778,8 @@ } /*}}}*/ -void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ +void FemModel::SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius){/*{{{*/ /*serialized vectors:*/ IssmDouble* RSLg_old=NULL; - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; IssmDouble tide_love_h, tide_love_k, fluid_love, moi_e, moi_p, omega, g; IssmDouble load_love_k2 = -0.30922675; //degree 2 load Love number @@ -4799,12 +4789,4 @@ /*Serialize vectors from previous iteration:*/ RSLg_old=pRSLg_old->ToMPISerial(); - - /*First, figure out the area of the ocean, which is needed to compute the eustatic component: */ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - eartharea_cpu += element->GetAreaSpherical(); - } - ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); IssmDouble moi_list[3]={0,0,0}; @@ -4872,11 +4854,8 @@ } /*}}}*/ -void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ +void FemModel::SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pRSLg, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/ /*serialized vectors:*/ IssmDouble* RSLg=NULL; - - IssmDouble eartharea=0; - IssmDouble eartharea_cpu=0; IssmDouble* Up = NULL; @@ -4898,12 +4877,4 @@ indices=xNew(gsize); for (int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - eartharea_cpu += element->GetAreaSpherical(); - } - ISSM_MPI_Reduce (&eartharea_cpu,&eartharea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&eartharea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - /*Call the sea level rise core: */ for(int i=0;iSize();i++){ @@ -4929,9 +4900,8 @@ } /*}}}*/ -IssmDouble FemModel::SealevelriseOceanAverage(Vector* RSLg) { /*{{{*/ +IssmDouble FemModel::SealevelriseOceanAverage(Vector* RSLg,IssmDouble oceanarea) { /*{{{*/ IssmDouble* RSLg_serial=NULL; IssmDouble oceanvalue,oceanvalue_cpu; - IssmDouble oceanarea,oceanarea_cpu; /*Serialize vectors from previous iteration:*/ @@ -4940,17 +4910,11 @@ /*Initialize:*/ oceanvalue_cpu=0; - oceanarea_cpu=0; /*Go through elements, and add contribution from each element and divide by overall ocean area:*/ for(int i=0;iSize();i++){ Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - oceanarea_cpu += element->OceanArea(); oceanvalue_cpu += element->OceanAverage(RSLg_serial); } - ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - - //_printf0_("Ocean area: " << oceanarea << "\n"); - + ISSM_MPI_Reduce (&oceanvalue_cpu,&oceanvalue,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); ISSM_MPI_Bcast(&oceanvalue,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24915) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 24916) @@ -162,9 +162,9 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseEustatic(Vector* pSgi, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); - void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); - void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); - void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); - IssmDouble SealevelriseOceanAverage(Vector* Sg); + void SealevelriseEustatic(Vector* pSgi, IssmDouble* peartharea, IssmDouble* poceanarea, IssmDouble* peustatic, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,int loop); + void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, IssmDouble eartharea, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius,bool verboseconvolution,int loop); + void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius); + void SealevelriseElastic(Vector* pUp, Vector* pNorth, Vector* pEast, Vector* pSg_old, IssmDouble eartharea,IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz); + IssmDouble SealevelriseOceanAverage(Vector* Sg,IssmDouble oceanarea); #endif void HydrologyEPLupdateDomainx(IssmDouble* pEplcount); Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 24915) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 24916) @@ -55,7 +55,7 @@ void geodetic_core(FemModel* femmodel); void steric_core(FemModel* femmodel); -Vector* sealevelrise_core_eustatic(FemModel* femmodel); -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic); -void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg); +Vector* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea); +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea); +void sealevelrise_core_elastic(Vector** pU_radial, Vector** pU_north,Vector** pU_east,FemModel* femmodel,Vector* RSLg,IssmDouble eartharea); void sealevelrise_core_viscous(Vector** pU_gia,Vector** pN_gia,FemModel* femmodel,Vector* RSLg); void sealevelrise_diagnostics(FemModel* femmodel,Vector* RSLg); Index: /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24915) +++ /issm/trunk-jpl/src/c/cores/sealevelrise_core.cpp (revision 24916) @@ -92,4 +92,5 @@ int geodetic=0; IssmDouble dt; + IssmDouble eartharea,oceanarea; /*Should we even be here?:*/ @@ -145,11 +146,11 @@ /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */ - RSLg_eustatic=sealevelrise_core_eustatic(femmodel); + RSLg_eustatic=sealevelrise_core_eustatic(femmodel,&eartharea,&oceanarea); /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */ - RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic); + RSLg=sealevelrise_core_noneustatic(femmodel,RSLg_eustatic,eartharea,oceanarea); /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */ - sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg); + sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,eartharea); /*compute viscosus (GIA) geodetic signatures:*/ @@ -303,5 +304,5 @@ } /*}}}*/ -Vector* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/ +Vector* sealevelrise_core_eustatic(FemModel* femmodel,IssmDouble* peartharea,IssmDouble* poceanarea){ /*{{{*/ /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/ @@ -317,4 +318,5 @@ IssmDouble *radius = NULL; int loop; + IssmDouble eartharea,oceanarea; /*outputs:*/ @@ -334,8 +336,8 @@ /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(RSLgi,&eustatic, latitude, longitude, radius,loop); //this computes + femmodel->SealevelriseEustatic(RSLgi,&eartharea, &oceanarea,&eustatic, latitude, longitude, radius,loop); //this computes /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi); + RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,oceanarea); /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the @@ -350,7 +352,11 @@ xDelete(longitude); xDelete(radius); + + /*Assign output pointers and return: */ + *peartharea=eartharea; + *poceanarea=oceanarea; return RSLgi; }/*}}}*/ -Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic){ /*{{{*/ +Vector* sealevelrise_core_noneustatic(FemModel* femmodel,Vector* RSLg_eustatic,IssmDouble eartharea,IssmDouble oceanarea){ /*{{{*/ /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5. @@ -419,5 +425,5 @@ /*call the non eustatic module: */ - femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop); + femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, eartharea,latitude,longitude, radius,verboseconvolution,loop); /*assemble solution vector: */ @@ -428,5 +434,5 @@ /*call rotational feedback module: */ RSLgo_rot = new Vector(gsize); RSLgo_rot->Assemble(); - femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius); + femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,eartharea,latitude,longitude,radius); RSLgo_rot->Assemble(); @@ -440,5 +446,5 @@ /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */ - RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo); + RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,oceanarea); /*RSLg is the sum of the eustatic term, and the ocean terms: */ @@ -477,5 +483,5 @@ return RSLg; } /*}}}*/ -void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg){ /*{{{*/ +void sealevelrise_core_elastic(Vector** pU_esa, Vector** pU_north_esa,Vector** pU_east_esa,FemModel* femmodel,Vector* RSLg,IssmDouble eartharea){ /*{{{*/ Vector *U_esa = NULL; @@ -515,5 +521,5 @@ /*call the elastic main modlule:*/ - femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz); + femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,eartharea,latitude,longitude,radius,xx,yy,zz,loop,horiz); /*Assign output pointers:*/