Changeset 19582


Ignore:
Timestamp:
09/25/15 09:08:47 (9 years ago)
Author:
Eric.Larour
Message:

CHG: fixed the last remaining bugs. Big bug in the ER computation. Went from array to double!

Location:
issm/trunk-jpl/src/c/modules/SurfaceMassBalancex
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r19566 r19582  
    77#include "../../toolkits/toolkits.h"
    88
    9 const double Pi =3.141592653589793238462643383279502884197169399375105820974944592308;
     9const double Pi = 3.141592653589793;
    1010
    1111void Gembx(FemModel* femmodel){  /*{{{*/
    … …  
    100100
    101101        // initialize
    102         IssmDouble F = 0, H=0, G=0, E;
    103 
     102        IssmDouble F = 0, H=0, G=0;
     103        const IssmDouble E = 0.09;        //[mm d-1] model time growth constant E
    104104        // convert T from K to ºC
    105105        T = T - 273.15;
    106106
    107107        // temperature coefficient F
    108         if(T>-6.0) F =  0.7 + ((T/-6) * 0.3);
    109         if(T<=-6.0 && T>-22.0) F =  1 - ((T+6)/-16 * 0.8);
    110         if(T<=-22.0 && T>=-40.0) F =  0.2 - ((T+22)/-18 * 0.2);
     108        if(T>-6.0) F =  0.7 + ((T/-6.0) * 0.3);
     109        if(T<=-6.0 && T>-22.0) F =  1 - ((T+6.0)/-16.0 * 0.8);
     110        if(T<=-22.0 && T>-40.0) F =  0.2 - ((T+22.0)/-18.0 * 0.2);
    111111
    112112        // density coefficient F
    113         if(d<150) H=1.0;
    114 
    115         if(d>=150 & d<400) H = 1 - ((d-150)/250);
     113        if(d<150.0) H=1.0;
     114
     115        if(d>=150.0 & d<400.0) H = 1 - ((d-150.0)/250.0);
    116116
    117117        // temperature gradient coefficient G
    … …  
    121121        if(dT >= 0.5  && dT < 0.7)  G = 0.9 + (((dT - 0.5)/0.2) * 0.1);
    122122        if(dT >= .7              )  G = 1.0;
    123 
    124         // model time growth constant E
    125         E = 0.09;        //[mm d-1]
    126123
    127124        // grouped coefficient Q
    … …  
    182179
    183180        /*Figure out grain size from effective grain radius: */
    184         gsz=xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsz[i]=re[i]*2;
     181        gsz=xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsz[i]=re[i]*2.0;
    185182
    186183        /*Convert dt from seconds to day: */
    187         dt=smb_dt/86400;
     184        dt=smb_dt/86400.0;
    188185
    189186        /*Determine liquid-water content in percentage: */
    … …  
    191188
    192189        //set maximum water content by mass to 9 percent (Brun, 1980)
    193         for(int i=0;i<m;i++)if(lwc[i]>9) lwc[i]=9.0;
     190        for(int i=0;i<m;i++)if(lwc[i]>9.0) lwc[i]=9.0;
    194191
    195192
    … …  
    224221
    225222                        //index for dentricity > 0 and T gradients < 5 degC m-1 and >= 5 degC m-1
    226                         if(dT[i]<5){
     223                        if(dT[i]<5.0){
    227224                                //determine coefficients
    228225                                IssmDouble A = - 2e8 * exp(-6e3 / T[i]) * dt;
    … …  
    242239
    243240                        // wet snow metamorphism
    244                         if(W[i]>0){
     241                        if(W[i]>0.0){
    245242
    246243                                //_printf_("D}ritic wet snow metamorphism\n");
    247244
    248245                                //determine coefficient
    249                                 IssmDouble D = (1/16) * pow(lwc[i],3.0) * dt;
     246                                IssmDouble D = (1.0/16.0) * pow(lwc[i],3.0) * dt;
    250247
    251248                                // new dendricity and sphericity for wet snow
    … …  
    255252
    256253                        // dendricity and sphericity can not be > 1 or < 0
    257                         if (gdn[i]<0)gdn[i]=0.0;
    258                         if (gsp[i]<0)gsp[i]=0.0;
    259                         if (gdn[i]>1)gdn[i]=1.0;
    260                         if (gsp[i]>1)gsp[i]=1.0;
     254                        if (gdn[i]<0.0)gdn[i]=0.0;
     255                        if (gsp[i]<0.0)gsp[i]=0.0;
     256                        if (gdn[i]>1.0)gdn[i]=1.0;
     257                        if (gsp[i]>1.0)gsp[i]=1.0;
    261258
    262259                        // determine new grain size (mm)
    … …  
    276273
    277274                        //Wet snow metamorphism (Brun, 1989)
    278                         if(W[i]>0){
     275                        if(W[i]>0.0){
    279276                                //_printf_("Nond}ritic wet snow metamorphism\n");
    280277                                //wet rate of change coefficient
    281                                 IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400);   // [mm^3 s^-1]
     278                                IssmDouble E = 1.28e-8 + (4.22e-10 * pow(lwc[i],3.0))* (dt *86400.0);   // [mm^3 s^-1]
    282279
    283280                                // calculate change in grain volume and convert to grain size
    284                                 gsz[i] = 2 * pow(3/(Pi * 4)*((4 / 3)*Pi*pow(gsz[i]/2,3.0) + E),1.0/3.0);
     281                                gsz[i] = 2.0 * pow(3.0/(Pi * 4.0)*((4.0/ 3.0)*Pi*pow(gsz[i]/2.0,3.0) + E),1.0/3.0);
    285282
    286283                        }
    287284
    288285                        // grains with sphericity == 1 can not have grain sizes > 2 mm (Brun, 1992)
    289                         if (gsp[i]==1.0 && gsz[i]>2) gsz[i]=2.0;
     286                        if (gsp[i]==1.0 && gsz[i]>2.0) gsz[i]=2.0;
    290287
    291288                        // grains with sphericity == 0 can not have grain sizes > 5 mm (Brun, 1992)
    … …  
    294291
    295292                //convert grain size back to effective grain radius:
    296                 re[i]=gsz[i]/2;
     293                re[i]=gsz[i]/2.0;
    297294        }
    298295       
    … …  
    505502        IssmDouble  dt;
    506503        IssmDouble max_fdt=0;
    507         IssmDouble  Ts;
     504        IssmDouble  Ts=0;
    508505        IssmDouble  L;
    509506        IssmDouble  eS;
    510         IssmDouble  Ri;
     507        IssmDouble  Ri=0;
    511508        IssmDouble  coefM;
    512509        IssmDouble  coefH;
    … …  
    690687
    691688        // PREALLOCATE ARRAYS BEFORE LOOP FOR IMPROVED PERFORMANCE
    692         T0 = xNew<IssmDouble>(m+2);
     689        T0 = xNewZeroInit<IssmDouble>(m+2);
    693690        Tu=xNew<IssmDouble>(m);
    694691        Td=xNew<IssmDouble>(m);
    … …  
    831828
    832829}  /*}}}*/
    833 void shortwave(IssmDouble* swf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
     830void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid){ /*{{{*/
    834831
    835832        // DISTRIBUTES ABSORBED SHORTWAVE RADIATION WITHIN SNOW/ICE
    … …  
    852849        // Outputs
    853850        //   swf     = absorbed shortwave radiation [W m-2]
     851        //
     852       
     853        /*outputs: */
     854        IssmDouble* swf=NULL;
    854855
    855856        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   shortwave module\n");
     857
     858        /*Initialize and allocate: */
     859        swf=xNewZeroInit<IssmDouble>(m);
     860
    856861
    857862        // SHORTWAVE FUNCTION
    … …  
    10201025                }
    10211026        }
     1027        /*Assign output pointers: */
     1028        *pswf=swf;
     1029
    10221030} /*}}}*/
    10231031void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble T_air, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid){ /*{{{*/
    … …  
    10751083        m=*pm;
    10761084
    1077         /*Allocate: */
     1085        // determine initial mass
    10781086        mInit=xNew<IssmDouble>(m);
     1087        for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
     1088        massinit=0; for(int i=0;i<m;i++)massinit+=mInit[i];
    10791089
    10801090        if (P > 0){
    1081                 // determine initial mass
    1082                 for(int i=0;i<m;i++) mInit[i]= d[i] * dz[i];
    1083        
     1091                       
    10841092
    10851093                if (T_air <= 273.15){ // if snow
    … …  
    11481156                // check for conservation of mass
    11491157                mass=0; for(int i=0;i<m;i++)mass+=d[i]*dz[i];
    1150                 massinit=0; for(int i=0;i<m;i++)massinit+=mInit[i];
    11511158
    11521159                mass_diff = mass - massinit - P;
    … …  
    11881195        IssmDouble* F=NULL;
    11891196        IssmDouble* flxDn=NULL;
    1190         IssmDouble* R=NULL;
    1191         IssmDouble* ER=NULL;
     1197        IssmDouble  ER=0;
    11921198        IssmDouble* EI=NULL;
    11931199        IssmDouble* EW=NULL;
    … …  
    12071213        IssmDouble Wi;
    12081214        int        D_size;
     1215        int         i;
    12091216
    12101217        /*outputs:*/
    … …  
    12201227        IssmDouble* gsp=*pgsp;
    12211228        int         n=*pn;
     1229        IssmDouble* R=0;
    12221230       
    12231231        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
    … …  
    12461254
    12471255        // initialize melt and runoff scalars
    1248         R = 0;          // runoff [kg]
     1256        Rsum = 0;       // runoff [kg]
    12491257        sumM = 0;       // total melt [kg]
    12501258        mAdd = 0;       // mass added/removed to/from base of model [kg]
    12511259        addE = 0;       // energy added/removed to/from base of model [J]
    12521260
    1253         // calculate temperature excess above 0 °C
     1261        // calculate temperature excess above 0 deg C
    12541262        exsT=xNewZeroInit<IssmDouble>(n);
    12551263        for(int i=0;i<n;i++) exsT[i]= fmax(0, T[i] - CtoK);        // [K] to [°C]
    … …  
    12641272        // check if any pore water
    12651273        if (cellsum(W,n) > 0){
    1266                 // _printf_("PORE WATER REFREEZE");
     1274                if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("      pore water refreeze\n");
    12671275                // calculate maximum freeze amount, maxF [kg]
    12681276                for(int i=0;i<n;i++) maxF[i] = fmax(0, -((T[i] - CtoK) * m[i] * CI) / LF);
    … …  
    12761284
    12771285                // if pore water froze in ice then adjust d and dz thickness
    1278                 for(int i=0;i<n;i++)if(d[i]>dIce)d[i]=m[i]/d[i];
     1286                for(int i=0;i<n;i++)if(d[i]>dIce)d[i]=dIce;
     1287                for(int i=0;i<n;i++) dz[i]= m[i]/d[i];
    12791288        }
    12801289
    12811290        // squeeze water from snow pack
    1282         exsW=xNew<IssmDouble>(n); for(int i=0;i<n;i++){
     1291        exsW=xNew<IssmDouble>(n);
     1292        for(int i=0;i<n;i++){
    12831293                Wi= (910 - d[i]) * Swi * (m[i] / d[i]);        // irreducible water content [kg]
    12841294                exsW[i] = fmax(0, W[i] - Wi);                  // water "squeezed" from snow [kg]
    … …  
    13051315                        while (cellsum(surpE,n) > 0){
    13061316                                // use surplus energy to increase the temperature of lower cell
    1307                                 T[i+i] = surpE[i] * m[i+1]/CI + T[i+i];
     1317                                T[i+1] = surpE[i] * m[i+1]/CI + T[i+1];
    13081318                                surpT[i+1] = fmax(0, (T[i+1] - CtoK - 159.1342));
    13091319                                surpE[i+1] = surpT[i+1] * CI / m[i+1];
    … …  
    13261336
    13271337                // initialize refreeze, runoff, flxDn and dW vectors [kg]
    1328                 F = xNewZeroInit<IssmDouble>(n);
    1329                 R = xNewZeroInit<IssmDouble>(n);
     1338                IssmDouble* F = xNewZeroInit<IssmDouble>(n);
     1339                IssmDouble* R=xNewZeroInit<IssmDouble>(n);
     1340
    13301341                for(int i=0;i<n;i++)dW[i] = 0;
    13311342                flxDn=xNewZeroInit<IssmDouble>(n+1); for(int i=0;i<n;i++)flxDn[i+1]=F[i];
    … …  
    14361447                celldelete(&gdn,n,D,D_size);
    14371448                celldelete(&gsp,n,D,D_size);
     1449                celldelete(&EI,n,D,D_size);
     1450                celldelete(&EW,n,D,D_size);
     1451                celldelete(&R,n,D,D_size);
    14381452                n=D_size;
    14391453       
    14401454                // calculate new grid lengths
    14411455                for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
     1456
     1457                //calculate Rsum:
     1458                Rsum=cellsum(R,n);
     1459       
    14421460        }
    14431461
    … …  
    14761494        xDelete<int>(D); D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999)D_size++; D=xNew<int>(D_size);
    14771495        D_size=0; for(int i=0;i<n;i++)if(m[i]!=99999){ D[D_size] = i; D_size++;}
     1496
     1497        celldelete(&m,n,D,D_size);
     1498        celldelete(&W,n,D,D_size);
     1499        celldelete(&dz,n,D,D_size);
     1500        celldelete(&d,n,D,D_size);
     1501        celldelete(&T,n,D,D_size);
     1502        celldelete(&a,n,D,D_size);
     1503        celldelete(&re,n,D,D_size);
     1504        celldelete(&gdn,n,D,D_size);
     1505        celldelete(&gsp,n,D,D_size);
     1506        celldelete(&EI,n,D,D_size);
     1507        celldelete(&EW,n,D,D_size);
     1508        n=D_size;
     1509
     1510        // check if any of the top 10 cell depths are too large
     1511        X=0;
     1512        for(int i=9;i>=0;i--){
     1513                if(dz[i]> 2* dzMin){
     1514                        X=i;
     1515                        break;
     1516                }
     1517        }
     1518       
     1519        i=0;
     1520        while(i<=X){
     1521                if (dz [i] > dzMin *2){
     1522
     1523                        // _printf_("dz > dzMin * 2");
     1524                        // split in two
     1525                        cellsplit(&dz, n, i,.5);
     1526                        cellsplit(&W, n, i,.5);
     1527                        cellsplit(&m, n, i,.5);
     1528                        cellsplit(&T, n, i,1.0);
     1529                        cellsplit(&d, n, i,1.0);
     1530                        cellsplit(&a, n, i,1.0);
     1531                        cellsplit(&EI, n, i,1.0);
     1532                        cellsplit(&EW, n, i,1.0);
     1533                        cellsplit(&re, n, i,1.0);
     1534                        cellsplit(&gdn, n, i,1.0);
     1535                        cellsplit(&gsp, n, i,1.0);
     1536                        n++;
     1537                        X=X+1;
     1538                }
     1539                else i++;
     1540        }
    14781541
    14791542        //// CORRECT FOR TOTAL MODEL DEPTH
    … …  
    15161579
    15171580        // calculate final mass [kg] and energy [J]
    1518         ER=xNew<IssmDouble>(n);
    1519         Rsum = cellsum(R,n);
     1581        sumER = Rsum * (LF + CtoK * CI);
    15201582        for(int i=0;i<n;i++)EI[i] = m[i] * T[i] * CI;
    1521         for(int i=0;i<n;i++)ER[i] = R[i] * (LF + CtoK * CI);
    15221583        for(int i=0;i<n;i++)EW[i] = W[i] * (LF + CtoK * CI);
    15231584
    15241585        mSum1 = cellsum(W,n) + cellsum(m,n) + Rsum;
    15251586        sumE1 = cellsum(EI,n) + cellsum(EW,n);
    1526         sumER = cellsum(ER,n);
    15271587
    15281588        dm = round(mSum0 - mSum1 + mAdd);
    15291589        dE = round(sumE0 - sumE1 - sumER +  addE);
    1530 
    1531         if (dm !=0  && dE !=0) _printf_("mass and energy are not conserved in melt equations");
    1532         else if (dm != 0) _printf_("mass is not conserved in melt equations");
    1533         else if (dE != 0) _printf_("energy is not conserved in melt equations");
    1534 
    1535         // W = round(W * 10000)/10000;
    1536         for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations");
    1537 
     1590       
     1591        for(int i=0;i<n;i++) if (W[i]<0) _error_("negative pore water generated in melt equations\n");
     1592
     1593        if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
     1594                        << "dm: " << dm << " dE: " << dE << "\n");
     1595
     1596       
    15381597        /*Free ressources:*/
    15391598        if(m)xDelete<IssmDouble>(m);
    15401599        if(EI)xDelete<IssmDouble>(EI);
     1600        if(EW)xDelete<IssmDouble>(EW);
    15411601        if(maxF)xDelete<IssmDouble>(maxF);
    15421602        if(dW)xDelete<IssmDouble>(dW);
    … …  
    15451605        if(surpT)xDelete<IssmDouble>(surpT);
    15461606        if(surpE)xDelete<IssmDouble>(surpE);
    1547         if(F)xDelete<IssmDouble>(F);
    15481607        if(flxDn)xDelete<IssmDouble>(flxDn);
    15491608        if(D)xDelete<int>(D);
    1550         if(R)xDelete<IssmDouble>(R);
    15511609        if(M)xDelete<IssmDouble>(M);
    15521610       
    … …  
    16301688        IssmDouble* obp = xNew<IssmDouble>(m);
    16311689        obp[0]=0;
    1632         for(int i=1;i<m;i++)obp[i]=cumdz[i]*d[i];
     1690        for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
    16331691       
    16341692        // calculate new snow/firn density for:
    … …  
    16551713                                // common variable
    16561714                                H = exp((-60/(T[i] * 8.314))) * obp[i] / pow(re[i]/1000,2.0);
    1657                                 c0 = 9.2E-9 * H;
    1658                                 c1 = 3.7E-9 * H;
     1715                                c0 = 9.2e-9 * H;
     1716                                c1 = 3.7e-9 * H;
    16591717                                break;
    16601718

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

    r19568 r19582  
    2525void grainGrowth(IssmDouble* pre, IssmDouble* pgdn, IssmDouble* pgsp, IssmDouble* T,IssmDouble* dz,IssmDouble* d, IssmDouble* W,IssmDouble smb_dt,int m,int aIdx, int sid);
    2626void albedo(IssmDouble* a,int aIdx, IssmDouble* re, IssmDouble* d, IssmDouble cldFrac, IssmDouble aIce, IssmDouble aSnow, IssmDouble* T, IssmDouble* W, IssmDouble P, IssmDouble EC, IssmDouble t0wet, IssmDouble t0dry, IssmDouble K, IssmDouble dt,int m, int sid);
    27 void shortwave(IssmDouble* swf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid);
     27void shortwave(IssmDouble** pswf, int swIdx, int aIdx, IssmDouble dsw, IssmDouble as, IssmDouble* d, IssmDouble* dz, IssmDouble* re, int m, int sid);
    2828void thermo(IssmDouble* pEC, IssmDouble* T, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlw, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, int sid);
    2929void accumulation(IssmDouble** pT, IssmDouble** pdz, IssmDouble** pd, IssmDouble** pW, IssmDouble** pa, IssmDouble** pre, IssmDouble** pgdn, IssmDouble** pgsp, int* pm, IssmDouble Ta, IssmDouble P, IssmDouble dzMin, IssmDouble aSnow, int sid);
Note: See TracChangeset for help on using the changeset viewer.