FloatingiceMeltingRatex.cpp File Reference

#include "./FloatingiceMeltingRatex.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
#include "./../../classes/Inputs2/DatasetInput2.h"

Go to the source code of this file.

Functions
void	FloatingiceMeltingRatex (FemModel *femmodel)
void	LinearFloatingiceMeltingRatex (FemModel *femmodel)
void	SpatialLinearFloatingiceMeltingRatex (FemModel *femmodel)
void	MismipFloatingiceMeltingRatex (FemModel *femmodel)
void	FloatingiceMeltingRateIsmip6x (FemModel *femmodel)
void	BeckmannGoosseFloatingiceMeltingRatex (FemModel *femmodel)

Function Documentation

FloatingiceMeltingRatex()

void FloatingiceMeltingRatex

(

FemModel *

femmodel

)

Definition at line 10 of file FloatingiceMeltingRatex.cpp.

                                                {/*{{{*/
 
   /*Intermediaties*/
   int  basalforcing_model;
 
   /*First, get BMB model from parameters*/
   femmodel->parameters->FindParam(&basalforcing_model,BasalforcingsEnum);
 
   /*branch to correct module*/
   switch(basalforcing_model){
      case FloatingMeltRateEnum:
      case MantlePlumeGeothermalFluxEnum:
         /*Nothing to be done*/
         break;
      case LinearFloatingMeltRateEnum:
         if(VerboseSolution())_printf0_("   call Linear Floating melting rate module\n");
         LinearFloatingiceMeltingRatex(femmodel);
         break;
      case MismipFloatingMeltRateEnum:
         if(VerboseSolution())_printf0_("   call Mismip Floating melting rate module\n");
         MismipFloatingiceMeltingRatex(femmodel);
         break;
      case SpatialLinearFloatingMeltRateEnum:
         if(VerboseSolution())_printf0_("   call Spatial Linear Floating melting rate module\n");
         SpatialLinearFloatingiceMeltingRatex(femmodel);
         break;
      case BasalforcingsPicoEnum:
         if(VerboseSolution())_printf0_("   call Pico Floating melting rate module\n");
         FloatingiceMeltingRatePicox(femmodel);
         break;
      case BasalforcingsIsmip6Enum:
         if(VerboseSolution())_printf0_("   call ISMIP 6 Floating melting rate module\n");
         FloatingiceMeltingRateIsmip6x(femmodel);
         break;
      case BeckmannGoosseFloatingMeltRateEnum:
         if(VerboseSolution())_printf0_("        call BeckmannGoosse Floating melting rate module\n");
         BeckmannGoosseFloatingiceMeltingRatex(femmodel);
         break;
      default:
         _error_("Basal forcing model "<<EnumToStringx(basalforcing_model)<<" not supported yet");
   }
 
}/*}}}*/

LinearFloatingiceMeltingRatex()

void LinearFloatingiceMeltingRatex

(

FemModel *

femmodel

)

Definition at line 54 of file FloatingiceMeltingRatex.cpp.

                                                      {/*{{{*/
 
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->LinearFloatingiceMeltingRate();
   }
 
}/*}}}*/

SpatialLinearFloatingiceMeltingRatex()

void SpatialLinearFloatingiceMeltingRatex

(

FemModel *

femmodel

)

Definition at line 62 of file FloatingiceMeltingRatex.cpp.

                                                             {/*{{{*/
 
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->SpatialLinearFloatingiceMeltingRate();
   }
 
}/*}}}*/

MismipFloatingiceMeltingRatex()

void MismipFloatingiceMeltingRatex

(

FemModel *

femmodel

)

Definition at line 70 of file FloatingiceMeltingRatex.cpp.

                                                      {/*{{{*/
 
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->MismipFloatingiceMeltingRate();
   }
}

FloatingiceMeltingRateIsmip6x()

void FloatingiceMeltingRateIsmip6x

(

FemModel *

femmodel

)

Definition at line 78 of file FloatingiceMeltingRatex.cpp.

                                                      {/*{{{*/
 
   int         num_basins, basinid, num_depths, domaintype;
   IssmDouble  area, tf, base, time;
   bool        islocal;
   IssmDouble* tf_depths = NULL;
 
   femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmip6NumBasinsEnum);
   femmodel->parameters->FindParam(&tf_depths,&num_depths,BasalforcingsIsmip6TfDepthsEnum); _assert_(tf_depths); 
   femmodel->parameters->FindParam(&islocal,BasalforcingsIsmip6IsLocalEnum);
 
   /*Binary search works for vectors that are sorted in increasing order only, make depths positive*/
   for(int i=0;i<num_depths;i++) tf_depths[i] = -tf_depths[i];
 
   IssmDouble* tf_weighted_avg     = xNewZeroInit<IssmDouble>(num_basins);
   IssmDouble* tf_weighted_avg_cpu = xNewZeroInit<IssmDouble>(num_basins);
   IssmDouble* areas_summed        = xNewZeroInit<IssmDouble>(num_basins);
   IssmDouble* areas_summed_cpu    = xNewZeroInit<IssmDouble>(num_basins);
 
   /*Get TF at each ice shelf point - linearly intepolate in depth and time*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element     = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      int      numvertices = element->GetNumberOfVertices();
 
      /*Set melt to 0 if non floating*/
      if(!element->IsIceInElement() || !element->IsFloating() || !element->IsOnBase()){
         IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
         element->AddInput2(BasalforcingsFloatingiceMeltingRateEnum,values,P1DGEnum);
         element->AddInput2(BasalforcingsIsmip6TfShelfEnum,values,P1DGEnum);
         xDelete<IssmDouble>(values);
         continue;
      }
 
      /*Get TF on all vertices*/
      IssmDouble*    tf_test        = xNew<IssmDouble>(numvertices);
      IssmDouble*    depth_vertices = xNew<IssmDouble>(numvertices);
      DatasetInput2* tf_input = element->GetDatasetInput2(BasalforcingsIsmip6TfEnum); _assert_(tf_input);
 
      element->GetInputListOnVertices(&depth_vertices[0],BaseEnum);
 
      Gauss* gauss=element->NewGauss();
      for(int iv=0;iv<numvertices;iv++){
         gauss->GaussVertex(iv);
 
         /*Find out where the ice shelf base is within tf_depths*/
         IssmDouble depth = -depth_vertices[iv]; /*NOTE: make sure we are dealing with depth>0*/
         int offset;
         int found=binary_search(&offset,depth,tf_depths,num_depths); 
         if(!found) _error_("depth not found");
 
         if (offset==-1){
            /*get values for the first depth: */
            _assert_(depth<=tf_depths[0]);
            tf_input->GetInputValue(&tf_test[iv],gauss,0);
         }
         else if(offset==num_depths-1){
            /*get values for the last time: */
            _assert_(depth>=tf_depths[num_depths-1]);
            tf_input->GetInputValue(&tf_test[iv],gauss,num_depths-1);
         }
         else {
            /*get values between two times [offset:offset+1], Interpolate linearly*/
            _assert_(depth>=tf_depths[offset] && depth<tf_depths[offset+1]);
            IssmDouble deltaz=tf_depths[offset+1]-tf_depths[offset];
            IssmDouble alpha2=(depth-tf_depths[offset])/deltaz;
            IssmDouble alpha1=(1.-alpha2);
            IssmDouble tf1,tf2;
            tf_input->GetInputValue(&tf1,gauss,offset);
            tf_input->GetInputValue(&tf2,gauss,offset+1);
            tf_test[iv] = alpha1*tf1 + alpha2*tf2;
         }
      }
 
      element->AddInput2(BasalforcingsIsmip6TfShelfEnum,tf_test,P1DGEnum);
      xDelete<IssmDouble>(tf_test);
      xDelete<IssmDouble>(depth_vertices);
      delete gauss;
   }
 
   if(!islocal) {
      /*Compute sums of tf*area and shelf-area per cpu*/
      for(int i=0;i<femmodel->elements->Size();i++){
         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
         if(!element->IsOnBase() || !element->IsIceInElement() || !element->IsFloating()) continue;
         /*Spawn basal element if on base to compute element area*/
         Element* basalelement = element->SpawnBasalElement();
         Input2* tf_input=basalelement->GetInput2(BasalforcingsIsmip6TfShelfEnum); _assert_(tf_input);
         basalelement->GetInputValue(&basinid,BasalforcingsIsmip6BasinIdEnum);
         Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
         tf_input->GetInputValue(&tf,gauss);
         delete gauss;
         area=basalelement->GetHorizontalSurfaceArea();
         tf_weighted_avg[basinid]+=tf*area;
         areas_summed[basinid]   +=area;  
         /*Delete spawned element if we are in 3D*/
         basalelement->FindParam(&domaintype,DomainTypeEnum);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
      }
 
      /*Syncronize across cpus*/
      ISSM_MPI_Allreduce(tf_weighted_avg,tf_weighted_avg_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
      ISSM_MPI_Allreduce(areas_summed,areas_summed_cpu,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
 
      /*Compute weighted means and save*/
      for(int k=0;k<num_basins;k++){tf_weighted_avg_cpu[k] = tf_weighted_avg_cpu[k]/areas_summed_cpu[k];}
      femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsIsmip6AverageTfEnum,tf_weighted_avg_cpu,num_basins));
   }
 
   /*Compute meltrates*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->Ismip6FloatingiceMeltingRate();
   }
 
   /*Cleanup and return */
   xDelete<IssmDouble>(tf_weighted_avg);
   xDelete<IssmDouble>(tf_weighted_avg_cpu);
   xDelete<IssmDouble>(areas_summed);
   xDelete<IssmDouble>(areas_summed_cpu);
   xDelete<IssmDouble>(tf_depths);
}

BeckmannGoosseFloatingiceMeltingRatex()

void BeckmannGoosseFloatingiceMeltingRatex

(

FemModel *

femmodel

)

Definition at line 200 of file FloatingiceMeltingRatex.cpp.

                                                              {/*{{{*/
 
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->BeckmannGoosseFloatingiceMeltingRate();
   }
}