FloatingiceMeltingRatePicox.cpp

Go to the documentation of this file.

 
#include "./FloatingiceMeltingRatePicox.h"
#include "../InputDuplicatex/InputDuplicatex.h"
#include "../../shared/shared.h"
#include "../../toolkits/toolkits.h"
 
void FloatingiceMeltingRatePicox(FemModel* femmodel){/*{{{*/
 
   int maxbox;
   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->AddInput2(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);
   }
 
   /*Optional buoyant plume melt rate parameterization (Lazeroms et al., 2018) */
   femmodel->parameters->FindParam(&isplume,BasalforcingsPicoIsplumeEnum);
   if(isplume) ComputeBasalMeltPlume(femmodel);
}/*}}}*/
 
void UpdateBoxIdsPico(FemModel* femmodel){/*{{{*/
 
   int         numvertices,num_basins,maxbox,basinid;
   IssmDouble  dist_max;
   IssmDouble* distances=NULL;
 
   femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
   femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
   IssmDouble* dmax_basin_cpu=xNew<IssmDouble>(num_basins);
 
   InputDuplicatex(femmodel,MaskOceanLevelsetEnum,DistanceToGroundinglineEnum);
   femmodel->DistanceToFieldValue(MaskOceanLevelsetEnum,0.,DistanceToGroundinglineEnum);
 
   InputDuplicatex(femmodel,MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
   femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0.,DistanceToCalvingfrontEnum);
 
   /*find maximum distance to grounding line per domain and per basin*/
   IssmDouble maxdist_cpu=-1.;
   for(int i=0;i<num_basins;i++){dmax_basin_cpu[i]=-1;}
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      if(!element->IsIceInElement() || !element->IsFloating()) continue;
      numvertices = element->GetNumberOfVertices();
      distances=xNew<IssmDouble>(numvertices);
      element->GetInputListOnVertices(&distances[0],DistanceToGroundinglineEnum);
      element->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
      for(int k=0; k<numvertices; k++){
         if(fabs(distances[k])>maxdist_cpu){maxdist_cpu=fabs(distances[k]);}
         if(fabs(distances[k])>dmax_basin_cpu[basinid]){dmax_basin_cpu[basinid]=fabs(distances[k]);}
      }
      xDelete<IssmDouble>(distances);
   }
 
   /*Synchronize across cpus*/
   IssmDouble* dmax_basin=xNew<IssmDouble>(num_basins);
   ISSM_MPI_Allreduce((void*)&maxdist_cpu,(void*)&dist_max,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
   ISSM_MPI_Allreduce(dmax_basin_cpu,dmax_basin,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
 
   /*Define maximum number of boxes per basin*/
   int* nd=xNew<int>(num_basins);
   for(int i=0; i<num_basins;i++){
      IssmDouble val=sqrt(dmax_basin[i]/dist_max)*(maxbox-1);
 
      #ifdef _HAVE_AD_
      _error_("Check the implementation of floor below");
      /*Do not use floor when AD is on*/
      int k=0; while(k<val+.5){k++;}
      nd[i]=k;
 
      #else
      nd[i]= reCast<int>(floor(val));
      #endif
   } 
 
   /*Assign box numbers*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->PicoUpdateBoxid(nd);
   }
 
   /*Cleanup and return */
   xDelete<int>(nd);
   xDelete<IssmDouble>(dmax_basin);
   xDelete<IssmDouble>(dmax_basin_cpu);
 
}/*}}}*/
void ComputeBoxAreasPico(FemModel* femmodel){/*{{{*/
 
   int num_basins,maxbox,basinid,boxid,domaintype;
   IssmDouble dist_max,area;
 
   femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
   femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
 
   IssmDouble* boxareas=xNewZeroInit<IssmDouble>(num_basins*maxbox);
 
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      if(!element->IsOnBase()) continue;
      Element* basalelement = element->SpawnBasalElement();
      if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
      basalelement->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
      basalelement->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
      boxareas[basinid*maxbox+boxid]+=basalelement->GetHorizontalSurfaceArea();
      basalelement->FindParam(&domaintype,DomainTypeEnum);
      if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
   }
 
   /*Synchronize across cpus*/
   IssmDouble* sumareas =xNew<IssmDouble>(num_basins*maxbox);
   ISSM_MPI_Allreduce(boxareas,sumareas,num_basins*maxbox,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
   //if(sumareas[0]==0){_error_("No elements in box 0, basal meltrates will be 0. Consider decreasing md.basalforcings.maxboxcount or refining your mesh!");}
 
   /*Update parameters to keep track of the new areas in future calculations*/
   femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoBoxAreaEnum,sumareas,maxbox*num_basins));
 
   /*Cleanup and return */
   xDelete<IssmDouble>(boxareas);
   xDelete<IssmDouble>(sumareas);
 
}/*}}}*/
void UpdateBoxPico(FemModel* femmodel, int loopboxid){/*{{{*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->PicoUpdateBox(loopboxid);
   }
}/*}}}*/
void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid){/*{{{*/
 
   int num_basins, basinid, maxbox, M, domaintype;
   IssmDouble area, toc, soc, overturning;
   IssmDouble* boxareas=NULL;
   IssmDouble* overturning_weighted_avg=NULL;
 
   femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
   femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
   femmodel->parameters->FindParam(&boxareas,&M, BasalforcingsPicoBoxAreaEnum);
   IssmDouble* toc_weighted_avg           = xNewZeroInit<IssmDouble>(num_basins);
   IssmDouble* soc_weighted_avg           = xNewZeroInit<IssmDouble>(num_basins);
   IssmDouble* toc_sumweightedavg         = xNewZeroInit<IssmDouble>(num_basins);
   IssmDouble* soc_sumweightedavg         = xNewZeroInit<IssmDouble>(num_basins);
   IssmDouble* overturning_sumweightedavg = xNewZeroInit<IssmDouble>(num_basins);
 
   /* Compute Toc and Soc weighted avg (boxes 0 to n-1) */
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      if(!element->IsOnBase()) continue;
      Element* basalelement = element->SpawnBasalElement();
      if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
      int el_boxid;
      basalelement->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
      if(el_boxid!=boxid) continue;
 
      Input2* tocs_input=basalelement->GetInput2(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input); 
      Input2* socs_input=basalelement->GetInput2(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
 
      basalelement->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
      Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
      tocs_input->GetInputValue(&toc,gauss);
      socs_input->GetInputValue(&soc,gauss);
      delete gauss;
      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;};
   }
 
   /*Syncronize across cpus*/
   ISSM_MPI_Allreduce(toc_weighted_avg,toc_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
   ISSM_MPI_Allreduce(soc_weighted_avg,soc_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
 
   for(int k=0;k<num_basins;k++){
      int p=k*maxbox+boxid; 
      if(boxareas[p]==0) continue;  
      toc_sumweightedavg[k] = toc_sumweightedavg[k]/boxareas[p];
      soc_sumweightedavg[k] = soc_sumweightedavg[k]/boxareas[p];
   }
 
   femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoAverageTemperatureEnum,toc_sumweightedavg,num_basins));
   femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoAverageSalinityEnum,soc_sumweightedavg,num_basins));
 
   /* Compute overturning weighted avg (box 0 only) */
   if(boxid==0){ 
      overturning_weighted_avg=xNewZeroInit<IssmDouble>(num_basins);
      for(int i=0;i<femmodel->elements->Size();i++){
         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
         if(!element->IsOnBase()) continue;
         Element* basalelement = element->SpawnBasalElement();
         if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
         int el_boxid;
         basalelement->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
         if(el_boxid!=boxid) continue;
 
         Input2* overturnings_input=basalelement->GetInput2(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
 
         basalelement->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
         Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
         overturnings_input->GetInputValue(&overturning,gauss);
         delete gauss;
         area=basalelement->GetHorizontalSurfaceArea();
         overturning_weighted_avg[basinid]+=overturning*area;
         basalelement->FindParam(&domaintype,DomainTypeEnum);
         if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
      }
 
      /*Syncronize across cpus*/
      ISSM_MPI_Allreduce(overturning_weighted_avg,overturning_sumweightedavg,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
 
      for(int k=0;k<num_basins;k++){
         int p=k*maxbox+boxid;
         if(boxareas[p]==0.) continue;
         overturning_sumweightedavg[k] = overturning_sumweightedavg[k]/boxareas[p];
      }
      femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoAverageOverturningEnum,overturning_sumweightedavg,num_basins));
   }
 
   /*Cleanup and return */
   xDelete<IssmDouble>(overturning_sumweightedavg);
   xDelete<IssmDouble>(toc_sumweightedavg);
   xDelete<IssmDouble>(soc_sumweightedavg);
   xDelete<IssmDouble>(overturning_weighted_avg);
   xDelete<IssmDouble>(toc_weighted_avg);
   xDelete<IssmDouble>(soc_weighted_avg);
   xDelete<IssmDouble>(boxareas);
}/*}}}*/
void ComputeBasalMeltPlume(FemModel* femmodel){/*{{{*/
   for(int i=0;i<femmodel->elements->Size();i++){
      Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
      element->PicoComputeBasalMelt();
   }
}/*}}}*/