Changeset 27241


Ignore:
Timestamp:
08/26/22 13:39:13 (3 years ago)
Author:
schlegel
Message:

CHG: incorporate some tolerances into GEMB initial layering

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r27240 r27241  
    101101        }
    102102        #endif
    103         if(dzTop < 0.05){
     103        if(dzTop < 0.05-Dtol){
    104104                _printf_("initial top grid cell length (dzTop) is < 0.05 m\n");
    105105        }
    … …  
    108108        //initialize top grid depth vector
    109109        dzT = xNew<IssmDouble>(gpTop);
    110         for (i=0;i<gpTop;i++)dzT[i]=dzTop;
     110        for (i=0;i<gpTop-Dtol;i++)dzT[i]=dzTop;
    111111
    112112        //bottom grid cell depth = x*zY^(cells from to structure)
    … …  
    115115        z0 = zTop;
    116116        gpBottom = 0;
    117         while (zMax > z0){
     117        while (zMax > z0+Dtol){
    118118                gp0= gp0 * zY;
    119119                z0 = z0 + gp0;
    … …  
    131131        //combine top and bottom dz vectors
    132132        dz = xNew<IssmDouble>(gpTop+gpBottom);
    133         for(i=0;i<gpTop;i++){
     133        for(i=0;i<gpTop-Dtol;i++){
    134134                dz[i]=dzT[i];
    135135        }
    … …  
    215215
    216216        //effective solar zenith angle
    217         IssmDouble x = min(pow(t/(3*cos(Pi*SZA/180.0)),0.5), 1.0);
    218         IssmDouble u = 0.64*x + (1-x)*cos(Pi*SZA/180.0);
     217        IssmDouble x = min(pow(t/(3.0*cos(Pi*SZA/180.0)),0.5), 1.0);
     218        IssmDouble u = 0.64*x + (1.0-x)*cos(Pi*SZA/180.0);
    219219
    220220        // pure snow albedo
    … …  
    243243
    244244                // determine the effective change due to finite depth and soot loading
    245                 IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1-as) - 0.1*c1 - 0.13)));
     245                IssmDouble A = min(1.0, (2.1 * pow(z1,1.35*(1.0-as) - 0.1*c1 - 0.13)));
    246246
    247247                dac =  (as2 + dac2 - as) + A*((as + dac) - (as2 + dac2));
    … …  
    249249
    250250        // change in albedo due to solar zenith angle
    251         IssmDouble dasz = 0.53*as*(1 - (as + dac))*pow(1 - u,1.2);
     251        IssmDouble dasz = 0.53*as*(1.0 - (as + dac))*pow(1.0 - u,1.2);
    252252
    253253        // change in albedo due to cloud (apart from change in diffuse fraction)
    254         IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1 + 1.5*t,as));
     254        IssmDouble dat = (0.1*t*pow(as + dac,1.3)) / (pow(1.0 + 1.5*t,as));
    255255
    256256        // Broadband albedo
    … …  
    551551
    552552        //some constants:
    553         const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
    554         const IssmDouble dPHC = 830;  //Pore closeoff density
     553        const IssmDouble dSnow = 300.0;   // density of fresh snow [kg m-3]       
     554        const IssmDouble dPHC = 830.0;  //Pore closeoff density
    555555        const IssmDouble ai_max = 0.58;  //maximum ice albedo, from Lefebre,2003
    556556        const IssmDouble ai_min = aIce;  //minimum ice albedo
    … …  
    565565                        // clabSnow, IssmDouble clabIce, IssmDouble SZA, IssmDouble COT, int m
    566566                        a[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, SZA, COT, dPHC, m);
    567                         adiff[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, 50, COT, dPHC, m);
     567                        adiff[0]=gardnerAlb(re, dz, d, clabSnow, clabIce, 50.0, COT, dPHC, m);
    568568                }
    569569                else if(aIdx==2){
    … …  
    688688                                        //  before run-off in kg m^-2 (melt per GEMB timestep, i.e. 3 hourly)
    689689                                        IssmDouble M = Msurf+W[0];
    690                                         a[0]=max(ai_min + (ai_max - ai_min)*exp(-1*(M/200)), ai_min);
     690                                        a[0]=max(ai_min + (ai_max - ai_min)*exp(-1.0*(M/200.0)), ai_min);
    691691
    692692                                }
    … …  
    696696
    697697        // Check for erroneous values
    698         if (a[0] > 1) _printf_("albedo > 1.0\n");
    699         else if (a[0] < 0) _printf_("albedo is negative\n");
     698        if (a[0] > 1.0+Ttol) _printf_("albedo > 1.0\n");
     699        else if (a[0] < 0.0-Dtol) _printf_("albedo is negative\n");
    700700        else if (xIsNan(a[0])) _error_("albedo == NAN\n");
    701701
    … …  
    959959
    960960        /* CALCULATE ENERGY SOURCES AND DIFFUSION FOR EVERY TIME STEP [dt]*/
    961         for (IssmDouble i=1;i<=dt0;i+=dt){
     961        for (IssmDouble i=1.0;i<=dt0;i+=dt){
    962962
    963963                // PART OF ENERGY CONSERVATION CHECK
    … …  
    10681068            }
    10691069            //zL = min(zL, 0.5); //Sjoblom, 2014
    1070             zLM=max(zL/Vz*z0,1.e-3);
    1071             zLT=max(zL/Tz*zT,1.e-3);
     1070            zLM=max(zL/Vz*z0,1e-3);
     1071            zLT=max(zL/Tz*zT,1e-3);
    10721072
    10731073            //Ding et al. 2020, from Beljaars and Holtslag (1991)
    … …  
    14661466        // MAIN FUNCTION
    14671467        // specify constants
    1468         IssmDouble dSnow = 150;    // density of snow [kg m-3]
     1468        IssmDouble dSnow = 150.0;    // density of snow [kg m-3]
    14691469        IssmDouble reNew = 0.05;    // new snow grain size [mm]
    14701470        IssmDouble gdnNew = 1.0;     // new snow dendricity
    … …  
    14741474        IssmDouble* mInit=NULL;
    14751475        bool        top=true;
    1476         IssmDouble  mass=0;
    1477         IssmDouble  massinit=0;
    1478         IssmDouble  mass_diff=0;
     1476        IssmDouble  mass=0.0;
     1477        IssmDouble  massinit=0.0;
     1478        IssmDouble  mass_diff=0.0;
    14791479
    14801480        /*output: */
    … …  
    15121512
    15131513                case 1: // Density of Antarctica snow
    1514                         dSnow = 350;
    1515                         //dSnow = 360; //FirnMICE Lundin et al., 2017
     1514                        dSnow = 350.0;
     1515                        //dSnow = 360.0; //FirnMICE Lundin et al., 2017
    15161516                        break;
    15171517
    15181518                case 2: // Density of Greenland snow, Fausto et al., 2018
    1519                         dSnow = 315;
     1519                        dSnow = 315.0;
    15201520                        //From Vionnet et al., 2012 (Crocus)
    15211521                        gdnNew = min(max(1.29 - 0.17*V,0.20),1.0);
    … …  
    15271527                        //dSnow = alpha1 + beta1*T + delta1*C + epsilon1*W
    15281528                        //     7.36x10-2  1.06x10-3  6.69x10-2  4.77x10-3
    1529                         dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK-Ttol) + 6.69e-2*C/1000 + 4.77e-3*Vmean)*1000;
     1529                        dSnow=(7.36e-2 + 1.06e-3*min(Tmean,CtoK-Ttol) + 6.69e-2*C/1000. + 4.77e-3*Vmean)*1000.;
    15301530                        break;
    15311531
    … …  
    17281728        // specify irreducible water content saturation [fraction]
    17291729        const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
    1730         const IssmDouble dPHC = 830;                     //Pore closeoff density
     1730        const IssmDouble dPHC = 830.0;                     //Pore closeoff density
    17311731
    17321732        //// REFREEZE PORE WATER
    17331733        // check if any pore water
    1734         if (cellsum(W,n) >0.0+Wtol){
     1734        if (cellsum(W,n) > 0.0+Wtol){
    17351735                if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("      pore water refreeze\n");
    17361736                // calculate maximum freeze amount, maxF [kg]
    … …  
    18961896                                        dMax = (dIce - d[i])*dz_0;          // maximum refreeze                                             
    18971897                                        IssmDouble maxF2 = min(dMax, maxF[i]-F1);      // maximum refreeze
    1898                                         F2 = min(-1*dW[i], maxF2);            // pore water refreeze
     1898                                        F2 = min(-1.0*dW[i], maxF2);            // pore water refreeze
    18991899                                        m[i] = m[i] + F2;                   // mass after refreeze
    19001900                                        d[i] = m[i]/dz_0;
    … …  
    21562156        for (int i=0;i<n;i++){
    21572157                if (i==0){
    2158                         dzMax2[i]=dzMin*2;
     2158                        dzMax2[i]=dzMin*2.0;
    21592159                }
    21602160                else{
    21612161                        Zcum[i]=Zcum[i-1]+dz[i];
    21622162                        if (Zcum[i]<=zTop+Dtol){
    2163                                 dzMax2[i]=dzMin*2;
     2163                                dzMax2[i]=dzMin*2.0;
    21642164                                X=i;
    21652165                        }
    21662166                        else{
    2167                                 dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2);
     2167                                dzMax2[i]=max(zY2*dzMin2[i-1],dzMin*2.0);
    21682168                        }
    21692169                }
    … …  
    24162416                                // ERA5 new aIdx=1, swIdx=0
    24172417                                if (aIdx==1 && swIdx==0){
    2418                                         if (fabs(adThresh - 820) < Dtol){
     2418                                        if (fabs(adThresh - 820.0) < Dtol){
    24192419                  // ERA5 v4
    24202420                  M0 = max(1.5131 - (0.1317 * log(C)),0.25);
    … …  
    24612461                                // ERA5 new aIdx=1, swIdx=0
    24622462                                if (aIdx==1 && swIdx==0){
    2463                                         if (fabs(adThresh - 820) < Dtol){
     2463                                        if (fabs(adThresh - 820.0) < Dtol){
    24642464                                                // ERA5 v4
    24652465                                                M0 = max(1.3566 - (0.1350 * log(C)),0.25);
Note: See TracChangeset for help on using the changeset viewer.