Changeset 27227

Timestamp:

08/25/22 11:19:06 (3 years ago)

Author:

schlegel

Message:

CHG: make GEMB cell merging its own function

Location:

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex

Files:

• Gembx.cpp (modified) (4 diffs)
• SurfaceMassBalancex.h (modified) (1 diff)

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

@@ -1828 +1828 @@
         xDelete<IssmDouble>(F);
 
+        //Manage the layering to match the user defined requirements
+        managelayers(&mAdd, &dz_add, &addE, &m, &EI, &EW, &T, &d, &dz, &W, &a, &adiff, &re, &gdn, &gsp, &n, dzMin, zMax, zMin, zTop, zY);
+
+        //// CHECK FOR MASS AND ENERGY CONSERVATION
+
+        // calculate final mass [kg] and energy [J]
+        sumER = Rsum * (LF + CtoK * CI);
+        for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
+        for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
+
+        mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
+        sumE1 = cellsum(EI,n) + cellsum(EW,n);
+
+        /*checks: */
+        for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
+
+        /*only in forward mode! avoid round in AD mode as it is not differentiable: */
+        #ifndef _HAVE_AD_
+        dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
+        dE = round(sumE0 - sumE1 - sumER +  addE - surplusE);
+        if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
+                        << "dm: " << dm << " dE: " << dE << "\n");
+        #endif
+
+        
+        /*Free resources:*/
+        if(m)xDelete<IssmDouble>(m);
+        if(EI)xDelete<IssmDouble>(EI);
+        if(EW)xDelete<IssmDouble>(EW);
+        if(maxF)xDelete<IssmDouble>(maxF);
+        if(dW)xDelete<IssmDouble>(dW);
+        if(exsW)xDelete<IssmDouble>(exsW);
+        if(exsT)xDelete<IssmDouble>(exsT);
+        if(surpT)xDelete<IssmDouble>(surpT);
+        if(surpE)xDelete<IssmDouble>(surpE);
+        if(flxDn)xDelete<IssmDouble>(flxDn);
+        if(D)xDelete<int>(D);
+        if(M)xDelete<IssmDouble>(M);
+
+        /*Assign output pointers:*/
+        *pMs=Msurf;
+        *pM=sumM;
+        *pR=Rsum;
+        *pF=Fsum;
+        *pmAdd=mAdd;
+        *pdz_add=dz_add;
+
+        *pT=T;
+        *pd=d;
+        *pdz=dz;
+        *pW=W;
+        *pa=a;
+        *padiff=adiff;
+        *pre=re;
+        *pgdn=gdn;
+        *pgsp=gsp;
+        *pn=n;
+
+} /*}}}*/ 
+void managelayers(IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble* paddE, IssmDouble** pm, IssmDouble** pEI, IssmDouble** pEW, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY){ /*{{{*/
+
+        /*intermediary:*/
+        IssmDouble* Zcum=NULL;
+        IssmDouble* dzMin2=NULL;
+        IssmDouble* dzMax2=NULL;
+        int*        D=NULL;
+
+        IssmDouble zY2=zY;
+        IssmDouble X=0.0;
+        int X1=0;
+        int X2=0;
+        int D_size = 0;
+
+        IssmDouble Ztot=0.0;
+        IssmDouble T_bot=0.0;
+        IssmDouble m_bot=0.0;
+        IssmDouble dz_bot=0.0;
+        IssmDouble d_bot=0.0;
+        IssmDouble W_bot=0.0;
+        IssmDouble a_bot=0.0;
+        IssmDouble adiff_bot=0.0;
+        IssmDouble re_bot=0.0;
+        IssmDouble gdn_bot=0.0;
+        IssmDouble gsp_bot=0.0;
+        IssmDouble EI_bot=0.0;
+        IssmDouble EW_bot=0.0;
+        bool       top=false;
+
+        /*outputs:*/
+        IssmDouble  mAdd = 0.0;
+        IssmDouble  addE = 0.0;
+        IssmDouble  dz_add = 0.0;
+        IssmDouble* T=*pT;
+        IssmDouble* d=*pd;
+        IssmDouble* dz=*pdz;
+        IssmDouble* W=*pW;
+        IssmDouble* a=*pa;
+        IssmDouble* adiff=*padiff;
+        IssmDouble* re=*pre;
+        IssmDouble* gdn=*pgdn;
+        IssmDouble* gsp=*pgsp;
+        IssmDouble* m=*pm;
+        IssmDouble* EI=*pEI;
+        IssmDouble* EW=*pEW;
+        int         n=*pn;
+
         //Merging of cells as they are burried under snow.
         Zcum=xNew<IssmDouble>(n);
@@ -2019 +2125 @@
         }
 
-        //// CHECK FOR MASS AND ENERGY CONSERVATION
-
-        // calculate final mass [kg] and energy [J]
-        sumER = Rsum * (LF + CtoK * CI);
-        for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
-        for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
-
-        mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
-        sumE1 = cellsum(EI,n) + cellsum(EW,n);
-
-        /*checks: */
-        for(int i=0;i<n;i++) if (W[i]<0.0-Wtol) _error_("negative pore water generated in melt equations\n");
-
-        /*only in forward mode! avoid round in AD mode as it is not differentiable: */
-        #ifndef _HAVE_AD_
-        dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
-        dE = round(sumE0 - sumE1 - sumER +  addE - surplusE);
-        if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
-                        << "dm: " << dm << " dE: " << dE << "\n");
-        #endif
-
         /*Free resources:*/
-        if(m)xDelete<IssmDouble>(m);
-        if(EI)xDelete<IssmDouble>(EI);
-        if(EW)xDelete<IssmDouble>(EW);
-        if(maxF)xDelete<IssmDouble>(maxF);
-        if(dW)xDelete<IssmDouble>(dW);
-        if(exsW)xDelete<IssmDouble>(exsW);
-        if(exsT)xDelete<IssmDouble>(exsT);
-        if(surpT)xDelete<IssmDouble>(surpT);
-        if(surpE)xDelete<IssmDouble>(surpE);
-        if(flxDn)xDelete<IssmDouble>(flxDn);
-        if(D)xDelete<int>(D);
-        if(M)xDelete<IssmDouble>(M);
         xDelete<IssmDouble>(Zcum);
         xDelete<IssmDouble>(dzMin2);
@@ -2058 +2131 @@
 
         /*Assign output pointers:*/
-        *pMs=Msurf;
-        *pM=sumM;
-        *pR=Rsum;
-        *pF=Fsum;
-        *pmAdd=mAdd;
-        *pdz_add=dz_add;
-
         *pT=T;
         *pd=d;
@@ -2075 +2141 @@
         *pgsp=gsp;
         *pn=n;
+        *pm=m;
+        *pEI=EI;
+        *pEW=EW;
+
+        *pmAdd=mAdd;
+        *paddE=addE;
+        *pdz_add=dz_add;
 
 } /*}}}*/ 

issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

@@ -37 +37 @@
 void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, IssmDouble* pRa, int* pm, int aIdx, int dsnowIdx, IssmDouble Tmean, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, IssmDouble C, IssmDouble V, IssmDouble Vmean, IssmDouble dIce, int sid);
 void melt(IssmDouble* pM, IssmDouble* pMs, IssmDouble* pR, IssmDouble* pF, IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble Ra, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY, IssmDouble dIce, int sid);
+void managelayers(IssmDouble* pmAdd, IssmDouble* pdz_add, IssmDouble* paddE, IssmDouble** pm, IssmDouble** pEI, IssmDouble** pEW, IssmDouble** pT, IssmDouble** pd, IssmDouble** pdz, IssmDouble** pW, IssmDouble** pa, IssmDouble** padiff, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pn, IssmDouble dzMin, IssmDouble zMax, IssmDouble zMin, IssmDouble zTop, IssmDouble zY);
 void densification(IssmDouble** pd,IssmDouble** pdz, IssmDouble* T, IssmDouble* re, int denIdx, int aIdx, int swIdx, IssmDouble adThresh, IssmDouble C, IssmDouble dt, IssmDouble Tmean, IssmDouble dIce, int m, int sid);
 void turbulentFlux(IssmDouble* pshf, IssmDouble* plhf, IssmDouble* pEC, IssmDouble Ta, IssmDouble Ts, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble ds, IssmDouble Ws, IssmDouble Vz, IssmDouble Tz, IssmDouble dIce, int sid);