source: issm/oecreview/Archive/24307-24683/ISSM-24507-24508.diff

Last change on this file was 24684, checked in by Mathieu Morlighem, 5 years ago

CHG: added new review
File size: 5.4 KB

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

     
    493493                        for(int i=0;i<m;i++)a[i] -= d_a[i];                            // new albedo
    494494
    495495                        // modification of albedo due to thin layer of snow or solid
    496                         // condensation (deposition) at the surface surface
     496                        // condensation (deposition) at the surface
    497497
    498498                        // check if condensation occurs & if it is deposited in solid phase
    499499                        if ( EC > 0.0 + Dtol && T[0] < 0.0-Ttol) P = P + (EC/dSnow) * 1000.0;  // add cond to precip [mm]
    … …  
    13391339
    13401340        /*outputs:*/
    13411341        IssmDouble  mAdd = 0.0;
     1342        IssmDouble  surplusE = 0.0;
     1343        IssmDouble  surplusT = 0.0;
    13421344        IssmDouble dz_add = 0.0;
    13431345        IssmDouble  Rsum = 0.0;
    13441346        IssmDouble* T=*pT;
    … …  
    14321434
    14331435                        int i = 0;
    14341436                        while (cellsum(surpE,n) > 0.0+Ttol && i<n){
    1435                                 // use surplus energy to increase the temperature of lower cell
    1436                                 T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
    14371437
    1438                                 exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
    1439                                 T[i+1] = min(CtoK, T[i+1]);
     1438                                if (i<n-1){
     1439                                        // use surplus energy to increase the temperature of lower cell
     1440                                        T[i+1] = surpE[i]/m[i+1]/CI + T[i+1];
    14401441
    1441                                 surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
    1442                                 surpE[i+1] = surpT[i+1] * CI * m[i+1];
     1442                                        exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1];
     1443                                        T[i+1] = min(CtoK, T[i+1]);
    14431444
     1445                                        surpT[i+1] = max(0.0, exsT[i+1] - LF/CI);
     1446                                        surpE[i+1] = surpT[i+1] * CI * m[i+1];
     1447                                }
     1448                                else{
     1449                                        surplusT=max(0.0, exsT[i] - LF/CI);
     1450                                        surplusE=surpE[i];
     1451                                        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0){
     1452                                                _printf0_(" WARNING: surplus energy at the base of GEMB column\n");
     1453                                        }
     1454                                }
     1455
    14441456                                // adjust current cell properties (again 159.1342 is the max T)
    14451457                                exsT[i] = LF/CI;
    14461458                                surpE[i] = 0.0;
    … …  
    16481660                        W[i+1] = W[i+1] + W[i];                     // combine liquid water
    16491661                        m[i+1] = m_new;                             // combine top masses
    16501662
    1651                         // set cell to 99999 for deletion
     1663                        // set cell to -99999 for deletion
    16521664                        m[i] = Delflag;
    16531665                }
    16541666        }
    … …  
    16551667
    16561668        //If last cell has to be merged
    16571669        if(lastCellFlag){
    1658          //find closest cell to merge with
     1670      //find closest cell to merge with
    16591671                for(int i=n-2;i>=0;i--){
    16601672                        if(m[i]!=Delflag){
    16611673                                X2=i;
    … …  
    16781690                W[X1] = W[X1] + W[X2];                     // combine liquid water
    16791691                m[X1] = m_new;                             // combine top masses
    16801692
    1681                 // set cell to 99999 for deletion
     1693                // set cell to -99999 for deletion
    16821694                m[X2] = Delflag;
    16831695        }
    16841696
    … …  
    17031715
    17041716        // check if any of the top 10 cell depths are too large
    17051717        X=0;
    1706         for(int i=9;i>=0;i--){
     1718        for(int i=min(9,n-1);i>=0;i--){
    17071719                if(dz[i]> 2.0*dzMin+Dtol){
    17081720                        X=i;
    17091721                        break;
    … …  
    17141726        while(j<=X){
    17151727                if (dz[j] > dzMin*2.0+Dtol){
    17161728
    1717                                 // _printf_("dz > dzMin * 2");
    1718                                 // split in two
    1719                                 cellsplit(&dz, n, j,.5);
    1720                                 cellsplit(&W, n, j,.5);
    1721                                 cellsplit(&m, n, j,.5);
    1722                                 cellsplit(&T, n, j,1.0);
    1723                                 cellsplit(&d, n, j,1.0);
    1724                                 cellsplit(&a, n, j,1.0);
    1725                                 cellsplit(&EI, n, j,.5);
    1726                                 cellsplit(&EW, n, j,.5);
    1727                                 cellsplit(&re, n, j,1.0);
    1728                                 cellsplit(&gdn, n, j,1.0);
    1729                                 cellsplit(&gsp, n, j,1.0);
    1730                                 n++;
    1731                                 X=X+1;
     1729                        // _printf_("dz > dzMin * 2");
     1730                        // split in two
     1731                        cellsplit(&dz, n, j,.5);
     1732                        cellsplit(&W, n, j,.5);
     1733                        cellsplit(&m, n, j,.5);
     1734                        cellsplit(&T, n, j,1.0);
     1735                        cellsplit(&d, n, j,1.0);
     1736                        cellsplit(&a, n, j,1.0);
     1737                        cellsplit(&EI, n, j,.5);
     1738                        cellsplit(&EW, n, j,.5);
     1739                        cellsplit(&re, n, j,1.0);
     1740                        cellsplit(&gdn, n, j,1.0);
     1741                        cellsplit(&gsp, n, j,1.0);
     1742                        n++;
     1743                        X=X+1;
    17321744                }
    17331745                else j++;
    17341746        }
    … …  
    18161828        /*only in forward mode! avoid round in AD mode as it is not differentiable: */
    18171829        #ifndef _HAVE_AD_
    18181830        dm = round((mSum0 - mSum1 + mAdd)*100.0)/100.0;
    1819         dE = round(sumE0 - sumE1 - sumER +  addE);
     1831        dE = round(sumE0 - sumE1 - sumER +  addE - surplusE);
    18201832        if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
    18211833                        << "dm: " << dm << " dE: " << dE << "\n");
    18221834        #endif
    … …  
    19771989                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
    19781990                                c0arth = 0.07 * H;
    19791991                                c1arth = 0.03 * H;
     1992                                //ERA-5
     1993                                //M0 = max(2.3999 - (0.2610 * log(C)),0.25);
     1994                                //M1 = max(2.7469 - (0.3228 * log(C)),0.25);
     1995                                //RACMO
    19801996                                M0 = max(1.6599 - (0.1724 * log(C)),0.25);
    19811997                                M1 = max(2.0102 - (0.2458 * log(C)),0.25);
     1998                                //From Ligtenberg
     1999                                //H = exp((-60000.0/(Tmean * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
    19822000                                //M0 = max(1.435 - (0.151 * log(C)),0.25);
    19832001                                //M1 = max(2.366 - (0.293 * log(C)),0.25);
    19842002                                c0 = M0*c0arth;
    … …  
    19912009                                H = exp((-60000.0/(T[i] * R)) + (42400.0/(Tmean * R))) * (C * 9.81);
    19922010                                c0arth = 0.07 * H;
    19932011                                c1arth = 0.03 * H;
     2012                                // ERA5
     2013                                //M0 = max(1.8920 - (0.1569 * log(C)),0.25);
     2014                                //M1 = max(2.5662 - (0.2274 * log(C)),0.25);
     2015                                // RACMO
     2016                                M0 = max(1.6201 - (0.1450 * log(C)),0.25);
     2017                                M1 = max(2.5577 - (0.2899 * log(C)),0.25);
     2018                                // From Kuipers Munneke
    19942019                                //M0 = max(1.042 - (0.0916 * log(C)),0.25);
    19952020                                //M1 = max(1.734 - (0.2039 * log(C)),0.25);
    1996                                 M0 = max(1.6201 - (0.1450 * log(C)),0.25);
    1997                                 M1 = max(2.5577 - (0.2899 * log(C)),0.25);
    19982021                                c0 = M0*c0arth;
    19992022                                c1 = M1*c1arth;
    20002023                                break;
Note: See TracBrowser for help on using the repository browser.