Changeset 24760


Ignore:
Timestamp:
04/29/20 14:07:15 (5 years ago)
Author:
schlegel
Message:

CHG: Gemb maintain initial layer depth profile throughout the simulation

Location:
issm/trunk-jpl
Files:
13 edited

Legend:

Unmodified
Added
Removed

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

    r24736 r24760  
    14231423        IssmDouble* Zcum=NULL;
    14241424        IssmDouble* dzMin2=NULL;
     1425        IssmDouble* dzMax2=NULL;
    14251426        IssmDouble zY2=zY;
    14261427        bool lastCellFlag = false;
    … …  
    16891690                celldelete(&W,n,D,D_size);
    16901691                celldelete(&d,n,D,D_size);
     1692                celldelete(&dz,n,D,D_size);
    16911693                celldelete(&T,n,D,D_size);
    16921694                celldelete(&a,n,D,D_size);
    … …  
    17101712        Zcum=xNew<IssmDouble>(n);
    17111713        dzMin2=xNew<IssmDouble>(n);
    1712 
     1714        dzMax2=xNew<IssmDouble>(n);
     1715
     1716        X=0;
    17131717        Zcum[0]=dz[0]; // Compute a cumulative depth vector
    1714 
    1715         for (int i=1;i<n;i++){
    1716                 Zcum[i]=Zcum[i-1]+dz[i];
    1717         }
    1718 
    17191718        for (int i=0;i<n;i++){
    1720                 if (Zcum[i]<=zTop+Dtol || i<1){
     1719                if (i==0){
    17211720                        dzMin2[i]=dzMin;
    17221721                }
    17231722                else{
    1724                         dzMin2[i]=zY2*dzMin2[i-1];
    1725                 }
    1726         }
    1727 
    1728         // check if depth is too small
    1729         X = 0;
    1730         for(int i=n-1;i>=0;i--){
    1731                 if(dz[i]<dzMin2[i]-Dtol){
    1732                         X=i;
    1733                         break;
    1734                 }
    1735         }
    1736 
    1737         //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
    1738         if(X==n-1){
    1739                 lastCellFlag = true;
    1740                 X = X-1;
    1741         }
    1742         else{
    1743                 lastCellFlag = false;
    1744         }
    1745 
    1746         for (int i = 0; i<=X;i++){
    1747                 if (dz[i] < dzMin2[i]-Dtol && d[i+1]<dIce-Dtol){  // merge top cells, don't merge with ice
    1748                         // _printf_("dz > dzMin * 2')
     1723                        Zcum[i]=Zcum[i-1]+dz[i];
     1724                        if (Zcum[i]<=zTop+Dtol){
     1725                                dzMin2[i]=dzMin;
     1726                                X=i;
     1727                        }
     1728                        else{
     1729                                dzMin2[i]=zY2*dzMin2[i-1];
     1730                        }
     1731                }
     1732        }
     1733
     1734        // Check to see if any cells are too small and need to be merged
     1735        for (int i=0; i<n; i++){
     1736                if ( (i<=X && dz[i]<dzMin-Dtol) || (i>X && dz[i]<dzMin2[i]-Dtol) ) {
     1737
     1738                        if (i==n-1){
     1739                                X2=i;
     1740                                //find closest cell to merge with
     1741                                for(int j=n-2;j>=0;j--){
     1742                                        if(m[j]!=Delflag){
     1743                                                X1=j;
     1744                                                break;
     1745                                        }
     1746                                }
     1747                        }
     1748                        else{
     1749                                X1=i+1;
     1750                                X2=i;
     1751                        }
     1752
    17491753                        // adjust variables as a linearly weighted function of mass
    1750                         IssmDouble m_new = m[i] + m[i+1];
    1751                         T[i+1] = (T[i]*m[i] + T[i+1]*m[i+1]) / m_new;
    1752                         a[i+1] = (a[i]*m[i] + a[i+1]*m[i+1]) / m_new;
    1753                         re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
    1754                         gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
    1755                         gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
     1754                        IssmDouble m_new = m[X2] + m[X1];
     1755                        T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
     1756                        a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
     1757                        re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
     1758                        gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
     1759                        gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
    17561760
    17571761                        // merge with underlying grid cell and delete old cell
    1758                         dz[i+1] = dz[i] + dz[i+1];                 // combine cell depths
    1759                         d[i+1] = m_new / dz[i+1];                   // combine top densities
    1760                         W[i+1] = W[i+1] + W[i];                     // combine liquid water
    1761                         m[i+1] = m_new;                             // combine top masses
     1762                        dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
     1763                        d[X1] = m_new / dz[X1];                   // combine top densities
     1764                        W[X1] = W[X1] + W[X2];                     // combine liquid water
     1765                        m[X1] = m_new;                             // combine top masses
    17621766
    17631767                        // set cell to -99999 for deletion
    1764                         m[i] = Delflag;
    1765                 }
    1766         }
    1767 
    1768         //If last cell has to be merged
    1769         if(lastCellFlag){
    1770       //find closest cell to merge with
    1771                 for(int i=n-2;i>=0;i--){
    1772                         if(m[i]!=Delflag){
    1773                                 X2=i;
    1774                                 X1=n-1;
    1775                                 break;
    1776                         }
    1777                 }
    1778 
    1779                 // adjust variables as a linearly weighted function of mass
    1780                 IssmDouble m_new = m[X2] + m[X1];
    1781                 T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
    1782                 a[X1] = (a[X2]*m[X2] + a[X1]*m[X1]) / m_new;
    1783                 re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
    1784                 gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
    1785                 gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
    1786 
    1787                 // merge with underlying grid cell and delete old cell
    1788                 dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
    1789                 d[X1] = m_new / dz[X1];                   // combine top densities
    1790                 W[X1] = W[X1] + W[X2];                     // combine liquid water
    1791                 m[X1] = m_new;                             // combine top masses
    1792 
    1793                 // set cell to -99999 for deletion
    1794                 m[X2] = Delflag;
     1768                        m[X2] = Delflag;
     1769                }
    17951770        }
    17961771
    … …  
    18141789        xDelete<int>(D);
    18151790
    1816         // check if any of the top 10 cell depths are too large
     1791        // check if any of the cell depths are too large
    18171792        X=0;
    1818         for(int i=min(9,n-1);i>=0;i--){
    1819                 if(dz[i]> 2.0*dzMin+Dtol){
    1820                         X=i;
    1821                         break;
    1822                 }
    1823         }
    1824 
    1825         int j=0;
    1826         while(j<=X){
    1827                 if (dz[j] > dzMin*2.0+Dtol){
    1828 
     1793        Zcum[0]=dz[0]; // Compute a cumulative depth vector
     1794        for (int i=0;i<n;i++){
     1795                if (i==0){
     1796                        dzMax2[i]=dzMin*2;
     1797                }
     1798                else{
     1799                        Zcum[i]=Zcum[i-1]+dz[i];
     1800                        if (Zcum[i]<=zTop+Dtol){
     1801                                dzMax2[i]=dzMin*2;
     1802                                X=i;
     1803                        }
     1804                        else{
     1805                                dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2);
     1806                        }
     1807                }
     1808        }
     1809
     1810        for (int j=n-1;j>=0;j--){
     1811                if ((j<X && dz[j] > dzMax2[j]+Dtol) || (dz[j] > dzMax2[j]*zY2+Dtol)){
    18291812                        // _printf_("dz > dzMin * 2");
    18301813                        // split in two
    … …  
    18411824                        cellsplit(&gsp, n, j,1.0);
    18421825                        n++;
    1843                         X=X+1;
    1844                 }
    1845                 else j++;
     1826                }
    18461827        }
    18471828

  • issm/trunk-jpl/src/m/classes/SMBgemb.m

    r24735 r24760  
    344344                        fielddisplay(self,'Aini','Initial albedo when restart [-]');
    345345                        fielddisplay(self,'Tini','Initial snow temperature when restart [K]');
    346                         fielddisplay(self,'Sizeini','Initial number of layers when restart [K]');
     346                        fielddisplay(self,'Sizeini','Initial number of layers when restart [-]');
    347347
    348348

  • issm/trunk-jpl/src/m/classes/SMBgemb.py

    r24735 r24760  
    158158        string = "%s\n%s" % (string, fielddisplay(self, 'Aini', 'Initial albedo when restart [-]'))
    159159        string = "%s\n%s" % (string, fielddisplay(self, 'Tini', 'Initial snow temperature when restart [K]'))
    160         string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [K]'))
     160        string = "%s\n%s" % (string, fielddisplay(self, 'Sizeini', 'Initial number of layers when restart [-]'))
    161161
    162162        #additional albedo parameters:

  • issm/trunk-jpl/test/NightlyRun/test243.m

    r24735 r24760  
    99% Use of Gemb method for SMB computation
    1010md.smb = SMBgemb(md.mesh,md.geometry);
    11 md.smb.dsnowIdx = 0;
     11md.smb.dsnowIdx = 1;
    1212
    1313%load hourly surface forcing date from 1979 to 2009:
    … …  
    4646md=solve(md,'Transient');
    4747
    48 nlayers=size(md.results.TransientSolution(end).SmbT,2);
     48nlayers=size(md.results.TransientSolution(1).SmbT,2);
     49for i=2:length(md.results.TransientSolution)
     50   nlayers=min(size(md.results.TransientSolution(i).SmbT,2), nlayers);
     51end
    4952
    5053%Fields and tolerances to track changes
    5154field_names      ={'Layers','SmbDz','SmbT' ,'SmbD' ,'SmbRe','SmbGdn','SmbGsp','SmbA' ,'SmbEC','SmbMassBalance','SmbMAdd','SmbDzAdd','SmbFAC','SmbMeanSHF','SmbMeanLHF','SmbMeanULW','SmbNetLW','SmbNetSW'};
    52 field_tolerances ={1e-12,1e-12,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11};
     55field_tolerances ={1e-12,2e-12,1e-12,1e-11,1e-11,2e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11};
    5356
    5457field_values={...

  • issm/trunk-jpl/test/NightlyRun/test243.py

    r24735 r24760  
    2121md.smb = SMBgemb()
    2222md.smb.setdefaultparameters(md.mesh, md.geometry)
    23 md.smb.dsnowIdx = 0
     23md.smb.dsnowIdx = 1
    2424
    2525#load hourly surface forcing date from 1979 to 2009:
    … …  
    6060md = solve(md, 'Transient')
    6161
    62 nlayers=np.size(md.results.TransientSolution[-1].SmbT,1)
     62nlayers=np.size(md.results.TransientSolution[0].SmbT,1)
     63for i in range(1, np.size(md.results.TransientSolution,0)):
     64    nlayers=np.minimum(np.size(md.results.TransientSolution[i].SmbT,1),nlayers)
    6365
    6466#Fields and tolerances to track changes
    6567field_names = ['Layers','SmbDz', 'SmbT', 'SmbD', 'SmbRe', 'SmbGdn', 'SmbGsp', 'SmbA', 'SmbEC', 'SmbMassBalance', 'SmbMAdd', 'SmbDzAdd', 'SmbFAC', 'SmbMeanSHF', 'SmbMeanLHF', 'SmbMeanULW', 'SmbNetLW', 'SmbNetSW']
    66 field_tolerances = [1e-12,1e-12,1e-12,1e-11,1e-11,1e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11]
     68field_tolerances = [1e-12,2e-12,1e-12,1e-11,1e-11,2e-11,1e-11,1e-12,1e-11,1e-12,1e-12,1e-12,1e-11,2e-11,2e-11,1e-11,9e-10,2e-11]
    6769#shape is different in python solution (fixed using reshape) which can cause test failure:
    6870field_values = [nlayers,md.results.TransientSolution[-1].SmbDz[0, 0:nlayers].reshape(1, -1),

  • issm/trunk-jpl/test/NightlyRun/test252.m

    r24739 r24760  
    99% Use of Gemb method for SMB computation
    1010md.smb = SMBgemb(md.mesh,md.geometry);
    11 md.smb.dsnowIdx = 0;
     11md.smb.dsnowIdx = 4;
    1212
    1313%load hourly surface forcing date from 1979 to 2009:
    … …  
    5454md=solve(md,'Transient');
    5555
    56 nlayers=size(md.results.TransientSolution(end).SmbT,2);
     56nlayers=size(md.results.TransientSolution(1).SmbT,2);
     57for i=2:length(md.results.TransientSolution)
     58   nlayers=min(size(md.results.TransientSolution(i).SmbT,2), nlayers);
     59end
    5760
    5861%Fields and tolerances to track changes
    … …  
    6265        'SmbDz4','SmbT4' ,'SmbD4' ,'SmbRe4','SmbGdn4','SmbGsp4','SmbA4' ,'SmbEC4','SmbMassBalance4','SmbMAdd4','SmbDzAdd4','SmbFAC4'};
    6366field_tolerances ={1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,...
    64                    1e-12,1e-12,1e-11,1e-10,1e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,...
    65                    1e-12,1e-12,2e-12,2e-11,1e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,...
     67                   1e-12,1e-12,1e-11,1e-10,2e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,...
     68                   1e-12,1e-12,2e-12,2e-11,4e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,...
    6669                   1e-11,1e-11,4e-11,4e-11,1e-12,3e-11,1e-12,1e-12,1e-10,1e-12,1e-12,2e-11};
    6770

  • issm/trunk-jpl/test/NightlyRun/test252.py

    r24739 r24760  
    2121md.smb = SMBgemb()
    2222md.smb.setdefaultparameters(md.mesh, md.geometry)
    23 md.smb.dsnowIdx = 0
     23md.smb.dsnowIdx = 4
    2424
    2525#load hourly surface forcing date from 1979 to 2009:
    … …  
    6969md = solve(md, 'Transient')
    7070
    71 nlayers=np.size(md.results.TransientSolution[-1].SmbT,1)
     71nlayers=np.size(md.results.TransientSolution[0].SmbT,1)
     72for i in range(1, np.size(md.results.TransientSolution,0)):
     73    nlayers=np.minimum(np.size(md.results.TransientSolution[i].SmbT,1),nlayers)
    7274
    7375#Fields and tolerances to track changes
    … …  
    7880               'SmbDz4', 'SmbT4', 'SmbD4', 'SmbRe4', 'SmbGdn4', 'SmbGsp4', 'SmbA4', 'SmbEC4', 'SmbMassBalance4', 'SmbMAdd4', 'SmbDzAdd4', 'SmbFAC4']
    7981field_tolerances = [1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,1e-12,
    80                    1e-12,1e-12,1e-11,1e-10,1e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,
    81                    1e-12,1e-12,2e-12,2e-11,1e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,
     82                   1e-12,1e-12,1e-11,1e-10,2e-11,1e-11,1e-12,4e-12,1e-12,1e-12,1e-12,1e-11,
     83                   1e-12,1e-12,2e-12,2e-11,4e-11,1e-11,1e-12,1e-11,1e-11,1e-12,1e-12,1e-11,
    8284                   1e-11,1e-11,4e-11,4e-11,1e-12,3e-11,1e-12,1e-12,1e-10,1e-12,1e-12,2e-11]
    8385

  • issm/trunk-jpl/test/NightlyRun/test253.m

    r24737 r24760  
    5858
    5959nlayers=size(md.results.TransientSolution(1).SmbT,2);
     60for i=2:length(md.results.TransientSolution)
     61        nlayers=min(size(md.results.TransientSolution(i).SmbT,2), nlayers);
     62end
    6063
    6164%Fields and tolerances to track changes

  • issm/trunk-jpl/test/NightlyRun/test253.py

    r24737 r24760  
    7373
    7474nlayers=np.size(md.results.TransientSolution[0].SmbT,1)
     75for i in range(1, np.size(md.results.TransientSolution,0)):
     76    nlayers=np.minimum(np.size(md.results.TransientSolution[i].SmbT,1),nlayers)
    7577
    7678#Fields and tolerances to track changes
Note: See TracChangeset for help on using the changeset viewer.