Changeset 24010 for issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp

Timestamp:

06/11/19 11:07:32 (6 years ago)

Author:

tpelle

Message:

NEW: ISMIP6 and PICO basal melting parameterizations working in HO, PICOP coming soon

File:

• issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (modified) (6 diffs)

issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp

@@ -12 +12 @@
         bool isplume;
 
-   /*First, reset all melt to 0 */
-   for(int i=0;i<femmodel->elements->Size();i++){
-              Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-              int numvertices = element->GetNumberOfVertices();
-              IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
-              element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
-              xDelete<IssmDouble>(values);
-           }
+        /*First, reset all melt to 0 */
+        for(int i=0;i<femmodel->elements->Size();i++){
+                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
+                int numvertices = element->GetNumberOfVertices();
+                IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
+                element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
+                xDelete<IssmDouble>(values);
+        }
 
         /*PICO melt rate parameterization (Reese et al., 2018)*/
-   femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
-   UpdateBoxIdsPico(femmodel);
-   ComputeBoxAreasPico(femmodel);
-   for(int i=0;i<maxbox;i++){
-              UpdateBoxPico(femmodel,i);
-              ComputeAverageOceanvarsPico(femmodel,i);
+        femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
+        UpdateBoxIdsPico(femmodel);
+        ComputeBoxAreasPico(femmodel);
+        for(int i=0;i<maxbox;i++){
+                UpdateBoxPico(femmodel,i);
+                ComputeAverageOceanvarsPico(femmodel,i);
         }
 
@@ -103 +103 @@
 void ComputeBoxAreasPico(FemModel* femmodel){/*{{{*/
 
-        int num_basins,maxbox,basinid,boxid;
+        int num_basins,maxbox,basinid,boxid,domaintype;
         IssmDouble dist_max,area;
 
@@ -113 +113 @@
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                if(!element->IsIceInElement() || !element->IsFloating()) continue;
-                element->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
-                element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
-                boxareas[basinid*maxbox+boxid]+=element->GetHorizontalSurfaceArea();
+                if(!element->IsOnBase()) continue;
+                Element* basalelement = element->SpawnBasalElement();
+                if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
+                basalelement->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
+                basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+                boxareas[basinid*maxbox+boxid]+=basalelement->GetHorizontalSurfaceArea();
+                basalelement->FindParam(&domaintype,DomainTypeEnum);
+                if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         }
 
@@ -140 +144 @@
 void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid){/*{{{*/
 
-        int num_basins, basinid, maxbox, M;
+        int num_basins, basinid, maxbox, M, domaintype;
         IssmDouble area, toc, soc, overturning;
         IssmDouble* boxareas=NULL;
@@ -157 +161 @@
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                if(!element->IsIceInElement() || !element->IsFloating()) continue;
+                if(!element->IsOnBase()) continue;
+                Element* basalelement = element->SpawnBasalElement();
+                if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
                 int el_boxid;
-                element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
+                basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
                 if(el_boxid!=boxid) continue;
 
-                Input* tocs_input=element->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input); 
-                Input* socs_input=element->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
-
-                element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
-                Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
+                Input* tocs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input); 
+                Input* socs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
+
+                basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+                Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
                 tocs_input->GetInputValue(&toc,gauss);
                 socs_input->GetInputValue(&soc,gauss);
                 delete gauss;
-                area=element->GetHorizontalSurfaceArea();
+                area=basalelement->GetHorizontalSurfaceArea();
                 toc_weighted_avg[basinid]+=toc*area;
                 soc_weighted_avg[basinid]+=soc*area;
+                basalelement->FindParam(&domaintype,DomainTypeEnum);
+                if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
         }
 
@@ -194 +202 @@
                 for(int i=0;i<femmodel->elements->Size();i++){
                         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
-                        if(!element->IsIceInElement() || !element->IsFloating()) continue;
+                        if(!element->IsOnBase()) continue;
+                        Element* basalelement = element->SpawnBasalElement();
+                        if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
                         int el_boxid;
-                        element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
+                        basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
                         if(el_boxid!=boxid) continue;
 
-                Input* overturnings_input=element->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
-
-                        element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
-                        Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
+                Input* overturnings_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
+
+                        basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
+                        Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
                         overturnings_input->GetInputValue(&overturning,gauss);
                         delete gauss;
-                        area=element->GetHorizontalSurfaceArea();
+                        area=basalelement->GetHorizontalSurfaceArea();
                         overturning_weighted_avg[basinid]+=overturning*area;
+                        basalelement->FindParam(&domaintype,DomainTypeEnum);
+                        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
                 }