Changeset 27491 for issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

Timestamp:

01/02/23 00:15:34 (2 years ago)

Author:

rueckamp

Message:

CHG: MinThickness Mask and minor changes for Debris Analysis

File:

• issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp (modified) (5 diffs)

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

@@ -522 +522 @@
         IssmDouble AccG=0.1;                    // m w.e. /100m
         IssmDouble AccMax=1.;                    // m w.e.
-        IssmDouble ReferenceElevation=2200.;     // m M&L
+        IssmDouble ReferenceElevation; 
         IssmDouble AblationDays=120.;            //
 
@@ -537 +537 @@
         IssmDouble Alphad=0.07;            //              debris albedo                 0.07
         IssmDouble Alphai=0.4;             //              ice ablbedo
+        IssmDouble Alphaeff;
         IssmDouble Ustar=0.16;             // ms^-1        friction velocity             0.16
         IssmDouble Ca=1000.;                // jkg^-1K^-1   specific heat capacity of air
@@ -552 +553 @@
         IssmDouble smb,yts,z,debris;
         IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
+        bool isdebris;
+        int domaintype;
+        femmodel->parameters->FindParam(&isdebris,TransientIsdebrisEnum);
 
         /*Get material parameters and constants */
@@ -558 +562 @@
         femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); 
         PhiD=0.;
-        //if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
+        if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
 
         /* Loop over all the elements of this partition */
@@ -573 +577 @@
                 /* Get inputs */
                 element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
+                element->FindParam(&domaintype,DomainTypeEnum);         
 
                 /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
                 for(int v=0;v<numvertices;v++){
 
-                        /*Get vertex elevation */
-                        z=surfacelist[v];
-
-                        /* Get debris cover */
-                        debris=debriscover[v];
-
-                        /*IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
-                        IssmDouble n=debris/dk;
-                        IssmDouble nmax=1000;
-                        IssmDouble Alphaeff;
-                        if(n>nmax){
-                                Alphaeff=Alphad;
-                        } else {
-                                Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
-                        }*/
-
-                        // M&L
-                       IssmDouble Alphaeff=Alphad;
-
-                        /* compute smb */
-                        for (int ismb=0;ismb<2;ismb++){
-                                if(ismb==0){
-                                        // calc a reference smb to identify accum and melt region; debris only develops in ablation area
-                                        debris=0.;
-                                }else{
-                                        // only in the meltregime debris develops
-                                        if(-MassBalanceCmDayDebris<0) debris=debriscover[v]; 
-                                }
-                                MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ 
-                                    (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ 
-                                    (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
-                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- 
-                                    ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ 
-                                    ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
-                                    WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ 
-                                    K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- 
-                                    ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* 
-                                    Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) 
+                        /*Get vertex elevation */
+                        z=surfacelist[v];
+
+                        /*Get top element*/
+                        //if(domaintype==Domain3DEnum){
+
+                        //}else{
+                        //      Alphaeff=Alphad;
+                        //      ReferenceElevation=2200.;     // m M&L                          
+                        //}
+
+                        /* compute smb */
+                        for (int ismb=0;ismb<2;ismb++){
+                                if(ismb==0){
+                                        // calc a reference smb to identify accum and melt region; debris only develops in ablation area
+                                        debris=0.;
+                                        PhiD=0.;
+                                }else{
+                                        // only in the meltregime debris develops
+                                        if(-MassBalanceCmDayDebris<1e-14) debris=debriscover[v]; 
+                                }
+                                if(debris<=0.) debris=0.;
+                                IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
+                                IssmDouble n=debris/dk;
+                                IssmDouble nmax=1000;
+                                IssmDouble Alphaeff;
+                                if(n>nmax){
+                                        Alphaeff=Alphad;
+                                } else {
+                                        Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
+                                }
+                                ReferenceElevation=3200.;     // m HEF
+
+
+                                Alphaeff=Alphad;
+                                ReferenceElevation=2200.;     // m M&L  
+
+                                MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+ 
+                                                        (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+ 
+                                                        (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
+                                                                        WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z- 
+                                                                                ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+ 
+                                                                        ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)* 
+                                                                                        WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/ 
+                                                                        K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z- 
+                                                                                                                ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)* 
+                                                                                        Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.) 
                                     -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
                         }