Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 19285) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 19286) @@ -449,4 +449,9 @@ endif #}}} +#Esmf sources {{{ +if ESMF +issm_sources += ./main/esmfbinders.cpp +endif +#}}} #Wrapper sources Index: /issm/trunk-jpl/src/c/main/esmfbinders.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 19286) +++ /issm/trunk-jpl/src/c/main/esmfbinders.cpp (revision 19286) @@ -0,0 +1,112 @@ +/*!\file: esmfbinder.cpp + * \brief: ESMF binders for ISSM. Binders developed initially for the GEOS-5 framework. + */ + + +#include "./issm.h" + +/*GCM specific declarations:*/ +const int GCMForcingNumTerms = 1; +const int GCMForcingTerms[GCMForcingNumTerms]= { SMBgcmEnum}; + +extern "C" { + + FemModel *femmodel; + + void InitializeISSM(int argc, char** argv, int** pelementsonlocalrank, int* pnumberofelements, ISSM_MPI_Comm comm_init){ /*{{{*/ + + int numberofelements; + int* elementsonlocalrank=NULL; + + /*Initialize femmodel from arguments provided command line: */ + femmodel = new FemModel(argc,argv,comm_init); + + /*Figure out the partition for elements, and return: */ + numberofelements=femmodel->elements->NumberOfElements(); + + elementsonlocalrank=xNewZeroInit(numberofelements); + for (int i=0;ielements->Size();i++){ + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + elementsonlocalrank[element->Sid()]=1; + } + + /*Some specific code here for the binding: */ + femmodel->parameters->SetParam(SMBgcmEnum,SurfaceforcingsEnum); //bypass SMB model, will be provided by GCM! + + /*Assign output pointers: */ + *pnumberofelements=numberofelements; + *pelementsonlocalrank=elementsonlocalrank; + + } /*}}}*/ + void RunISSM(IssmDouble dt, IssmDouble* gcmforcings, IssmDouble* issmoutputs){ /*{{{*/ + + int numberofelements; + IssmDouble yts; + IssmDouble rho_ice; + IssmDouble area; + IssmDouble start_time,final_time; + + /*Figure out number of elements: */ + numberofelements=femmodel->elements->NumberOfElements(); + + /*Fetch some necessary constants: */ + femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Setup gcm forcings as element-wise input: {{{ */ + for (int f=0;felements->Size();i++){ + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + + switch(forcing_type){ + case SMBgcmEnum: + /*{{{*/ + { + + /*Recover rho_ice: */ + rho_ice=element->matpar->GetMaterialParameter(MaterialsRhoIceEnum); + + /*Recover area of element: */ + area=element->SurfaceArea(); + + /*Recover smb forcing from the gcm forcings: */ + IssmDouble smbforcing=*(gcmforcings+f*numberofelements+i); + + /*Convert to SI: */ + smbforcing=smbforcing*yts/(rho_ice*area); + + /*Add into the element as new forcing :*/ + element->inputs->AddInput(new DoubleInput(SurfaceforcingsMassBalanceEnum,smbforcing)); + + } + /*}}}*/ + break; + default: + { _error_("Unknown forcing type " << forcing_type << "\n"); } + break; + } + } + } + + /*}}}*/ + + /*Before running, setup the time interval: */ + femmodel->parameters->FindParam(&start_time,TimesteppingStartTimeEnum); + final_time=start_time+dt; + femmodel->parameters->SetParam(final_time,TimesteppingFinalTimeEnum); //we are bypassing ISSM's initial final time! + + /*Now, run: */ + femmodel->Solve(); + + /*For the next time around, save the final time as start time */ + femmodel->parameters->SetParam(final_time,TimesteppingStartTimeEnum); + } /*}}}*/ + void FinalizeISSM(){ /*{{{*/ + + /*Wrap up: */ + delete femmodel; femmodel=NULL; + } /*}}}*/ + +} Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 19285) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (revision 19286) @@ -59,4 +59,7 @@ if(VerboseSolution())_printf0_(" call smb Melt Components module\n"); SmbMeltComponentsx(femmodel); + break; + case SMBgcmEnum: + /*Nothing to be done*/ break; default: Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19285) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 19286) @@ -382,4 +382,5 @@ SurfaceforcingsMeltEnum, SurfaceforcingsRefreezeEnum, + SMBgcmEnum, /*TO BE REMOVED*/ SurfaceforcingsIspddEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19285) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 19286) @@ -386,4 +386,5 @@ case SurfaceforcingsMeltEnum : return "SurfaceforcingsMelt"; case SurfaceforcingsRefreezeEnum : return "SurfaceforcingsRefreeze"; + case SMBgcmEnum : return "SMBgcm"; case SurfaceforcingsIspddEnum : return "SurfaceforcingsIspdd"; case SurfaceforcingsIssmbgradientsEnum : return "SurfaceforcingsIssmbgradients"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19285) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 19286) @@ -395,4 +395,5 @@ else if (strcmp(name,"SurfaceforcingsMelt")==0) return SurfaceforcingsMeltEnum; else if (strcmp(name,"SurfaceforcingsRefreeze")==0) return SurfaceforcingsRefreezeEnum; + else if (strcmp(name,"SMBgcm")==0) return SMBgcmEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"Matdamageice")==0) return MatdamageiceEnum; else if (strcmp(name,"Matpar")==0) return MatparEnum; - else if (strcmp(name,"Node")==0) return NodeEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; + if (strcmp(name,"Node")==0) return NodeEnum; + else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum; else if (strcmp(name,"NumericalfluxType")==0) return NumericalfluxTypeEnum; else if (strcmp(name,"Param")==0) return ParamEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"StressTensorxy")==0) return StressTensorxyEnum; else if (strcmp(name,"StressTensorxz")==0) return StressTensorxzEnum; - else if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; + if (strcmp(name,"StressTensoryy")==0) return StressTensoryyEnum; + else if (strcmp(name,"StressTensoryz")==0) return StressTensoryzEnum; else if (strcmp(name,"StressTensorzz")==0) return StressTensorzzEnum; else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum; @@ -751,9 +752,9 @@ else if (strcmp(name,"Outputdefinition63")==0) return Outputdefinition63Enum; else if (strcmp(name,"Outputdefinition64")==0) return Outputdefinition64Enum; - else if (strcmp(name,"Outputdefinition65")==0) return Outputdefinition65Enum; else stage=7; } if(stage==7){ - if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum; + if (strcmp(name,"Outputdefinition65")==0) return Outputdefinition65Enum; + else if (strcmp(name,"Outputdefinition66")==0) return Outputdefinition66Enum; else if (strcmp(name,"Outputdefinition67")==0) return Outputdefinition67Enum; else if (strcmp(name,"Outputdefinition68")==0) return Outputdefinition68Enum; @@ -874,9 +875,9 @@ else if (strcmp(name,"Regular")==0) return RegularEnum; else if (strcmp(name,"Scaled")==0) return ScaledEnum; - else if (strcmp(name,"Separate")==0) return SeparateEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"Sset")==0) return SsetEnum; + if (strcmp(name,"Separate")==0) return SeparateEnum; + else if (strcmp(name,"Sset")==0) return SsetEnum; else if (strcmp(name,"Verbose")==0) return VerboseEnum; else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; Index: /issm/trunk-jpl/src/m/enum/EnumDefinitions.py =================================================================== --- /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19285) +++ /issm/trunk-jpl/src/m/enum/EnumDefinitions.py (revision 19286) @@ -378,4 +378,5 @@ def SurfaceforcingsMeltEnum(): return StringToEnum("SurfaceforcingsMelt")[0] def SurfaceforcingsRefreezeEnum(): return StringToEnum("SurfaceforcingsRefreeze")[0] +def SMBgcmEnum(): return StringToEnum("SMBgcm")[0] def SurfaceforcingsIspddEnum(): return StringToEnum("SurfaceforcingsIspdd")[0] def SurfaceforcingsIssmbgradientsEnum(): return StringToEnum("SurfaceforcingsIssmbgradients")[0]