Ignore:
Timestamp:
11/15/23 12:14:04 (16 months ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 28011

Location:
issm/trunk
Files:
5 edited
1 copied

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

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

    r27347 r28013  
    240240
    241241                // change in pure snow albedo due to soot loading
    242                 IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S1,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
     242                IssmDouble dac2 = max(0.04 - as2, pow(-c2,0.55)/(0.16 + 0.6*pow(S2,0.5) + (1.8*pow(c2,0.6))*pow(S2,-0.25)));
    243243
    244244                // determine the effective change due to finite depth and soot loading

  • issm/trunk/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

    r27347 r28013  
    6767                        }
    6868
    69                         smb[v]=smb[v]/1000*rho_water/rho_ice;      // SMB in m/y ice
    7069                }  //end of the loop over the vertices
    7170
    … …  
    187186        femmodel->parameters->FindParam(&arlagcoefs,&M,&N,SmbARMAarlagcoefsEnum);             _assert_(M==numbasins); _assert_(N==arorder);
    188187   femmodel->parameters->FindParam(&malagcoefs,&M,&N,SmbARMAmalagcoefsEnum);             _assert_(M==numbasins); _assert_(N==maorder);
    189    femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum);                 _assert_(M==numbasins); _assert_(N==numelevbins);
    190    femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum);                _assert_(M==numbasins); _assert_(N==numelevbins-1);
     188   femmodel->parameters->FindParam(&lapserates,&M,&N,SmbLapseRatesEnum);                 _assert_(M==numbasins); _assert_(N==numelevbins*12);
     189   femmodel->parameters->FindParam(&elevbins,&M,&N,SmbElevationBinsEnum);                _assert_(M==numbasins); _assert_(N==(numelevbins-1)*12);
    191190   femmodel->parameters->FindParam(&refelevation,&M,SmbRefElevationEnum);                _assert_(M==numbasins);
    192191
    … …  
    504503
    505504}/*}}}*/
    506 void SmbDebrisMLx(FemModel* femmodel){/*{{{*/
    507 
    508         //      The function is based on:
    509         //      Evatt GW, Abrahams ID, Heil M, Mayer C, Kingslake J, Mitchell SL, et al. Glacial melt under a porous debris layer. Journal of Glaciology 61 (2015) 825–836, doi:10.3189/2
    510         //      Constants/Values are taken from Mayer, Licciulli (2021): https://www.frontiersin.org/articles/10.3389/feart.2021.710276/full#B7
    511         //      function taken from https://github.com/carlolic/DebrisExp/blob/main/USFs/USF_DebrisCoverage.f90
    512 
    513         /*Intermediaries*/
    514         // altitude gradients of the crucial parameters (radiation from Marty et al., TaAClimat; 2002)
    515         IssmDouble LW=2.9;          // W/m^2 /100m                       2.9
    516         IssmDouble SW=1.3;          // W/m^2 /100m                       1.3
    517         IssmDouble HumidityG=0;     // % /100m         rough estimate
    518         IssmDouble AirTemp=0.7;     // C /100m
    519         IssmDouble WindSpeed=0.02;  // m/s /100m       rough estimate    0.2
    520 
    521         // accumulation follows a linear increase above the ELA up to a plateau
    522         IssmDouble AccG=0.1;                    // m w.e. /100m
    523         IssmDouble AccMax=1.;                    // m w.e.
    524         IssmDouble ReferenceElevation=2200.;     // m M&L
    525         IssmDouble AblationDays=120.;            //
    526 
    527         IssmDouble In=100.;                 // Wm^-2        incoming long wave
    528         IssmDouble Q=500.;                  // Wm^-2        incoming short wave
    529         IssmDouble K=0.585;                // Wm^-1K^-1    thermal conductivity          0.585
    530         IssmDouble Qm=0.0012;              // kg m^-3      measured humiditiy level
    531         IssmDouble Qh=0.006 ;              // kg m^-3      saturated humidity level
    532         IssmDouble Tm=2.;                   // C            air temperature
    533         IssmDouble Rhoaa=1.22;             // kgm^-3       air densitiy
    534         IssmDouble Um=1.5;                 // ms^-1        measured wind speed
    535         IssmDouble Xm=1.5;                 // ms^-1        measurement height
    536         IssmDouble Xr=0.01;                // ms^-1        surface roughness             0.01
    537         IssmDouble Alphad=0.07;            //              debris albedo                 0.07
    538         IssmDouble Alphai=0.4;             //              ice ablbedo
    539         IssmDouble Ustar=0.16;             // ms^-1        friction velocity             0.16
    540         IssmDouble Ca=1000.;                // jkg^-1K^-1   specific heat capacity of air
    541         IssmDouble Lm;//=3.34E+05;            // jkg^-1K^-1   latent heat of ice melt
    542         IssmDouble Lv=2.50E+06;            // jkg^-1K^-1   latent heat of evaporation
    543         IssmDouble Tf=273.;                 // K            water freeezing temperature
    544         IssmDouble Eps=0.95;               //              thermal emissivity
    545         IssmDouble Rhoi=900.;               // kgm^-3       ice density
    546         IssmDouble Sigma=5.67E-08;         // Wm^-2K^-4    Stefan Boltzmann constant
    547         IssmDouble Kstar=0.4;              //              von kármán constant
    548         IssmDouble Gamma=180.;              // m^-1         wind speed attenuation        234
    549         IssmDouble PhiD;//=0.005;              //              debris packing fraction       0.01
    550         IssmDouble Humidity=0.2;           //              relative humidity
    551 
    552         IssmDouble smb,yts,z,debris;
    553         IssmDouble MassBalanceCmDayDebris,MassBalanceMYearDebris;
    554 
    555         /*Get material parameters and constants */
    556         //femmodel->parameters->FindParam(&Rhoi,MaterialsRhoIceEnum); // Note Carlo's model used as  benchmark was run with different densities for debris and FS
    557         femmodel->parameters->FindParam(&Lm,MaterialsLatentheatEnum);
    558         femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
    559         PhiD=0.;
    560         //if(isdebris) femmodel->parameters->FindParam(&PhiD,DebrisPackingFractionEnum);
    561 
    562         /* Loop over all the elements of this partition */
    563         for(Object* & object : femmodel->elements->objects){
    564                 Element* element=xDynamicCast<Element*>(object);
    565 
    566                 /* Allocate all arrays */
    567                 int         numvertices=element->GetNumberOfVertices();
    568                 IssmDouble* surfacelist=xNew<IssmDouble>(numvertices);
    569                 IssmDouble* smb=xNew<IssmDouble>(numvertices);
    570                 IssmDouble* debriscover=xNew<IssmDouble>(numvertices);
    571                 element->GetInputListOnVertices(surfacelist,SurfaceEnum);
    572 
    573                 /* Get inputs */
    574                 element->GetInputListOnVertices(debriscover,DebrisThicknessEnum);
    575 
    576                 /*Loop over all vertices of element and calculate SMB as function of Debris Cover and z */
    577                 for(int v=0;v<numvertices;v++){
    578 
    579                         /*Get vertex elevation */
    580                         z=surfacelist[v];
    581 
    582                         /* Get debris cover */
    583                         debris=debriscover[v];
    584 
    585                         /*IssmDouble dk=1e-5; // TODO make Alphad and Alphai a user input
    586                         IssmDouble n=debris/dk;
    587                         IssmDouble nmax=1000;
    588                         IssmDouble Alphaeff;
    589                         if(n>nmax){
    590                                 Alphaeff=Alphad;
    591                         } else {
    592                                 Alphaeff=Alphai+n*(Alphad-Alphai)/nmax;
    593                         }*/
    594 
    595                         // M&L
    596                        IssmDouble Alphaeff=Alphad;
    597 
    598                         /* compute smb */
    599                         for (int ismb=0;ismb<2;ismb++){
    600                                 if(ismb==0){
    601                                         // calc a reference smb to identify accum and melt region; debris only develops in ablation area
    602                                         debris=0.;
    603                                 }else{
    604                                         // only in the meltregime debris develops
    605                                         if(-MassBalanceCmDayDebris<0) debris=debriscover[v];
    606                                 }
    607                                 MassBalanceCmDayDebris=(((In-(z-ReferenceElevation)*LW/100.)-(Eps*Sigma*(Tf*Tf*Tf*Tf))+
    608                                     (Q+(z-ReferenceElevation)*SW/100.)*(1.-Alphaeff)+
    609                                     (Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
    610                                     WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))*(Tm-(z-
    611                                     ReferenceElevation)*AirTemp/100.))/((1-PhiD)*Rhoi*Lm)/(1.+
    612                                     ((Rhoaa*Ca*Ustar*Ustar)/((Um-(z-ReferenceElevation)*
    613                                     WindSpeed/100.)-Ustar*(2.-(exp(Gamma*Xr))))+4.*Eps*Sigma*(Tf*Tf*Tf))/
    614                                     K*debris)-(Lv*Ustar*Ustar*(Qh-(Qh*(Humidity-(z-
    615                                     ReferenceElevation)*HumidityG/100.)))*(exp(-Gamma*Xr)))/((1.-PhiD)*
    616                                     Rhoi*Lm*Ustar)/((((Um-(z-ReferenceElevation)*WindSpeed/100.)
    617                                     -2.*Ustar)*exp(-Gamma*Xr))/Ustar+exp(Gamma*debris)))*100.*24.*60.*60.;
    618                         }
    619 
    620                         /* account form ablation days, and convert to m/s */
    621                         MassBalanceMYearDebris=-MassBalanceCmDayDebris/100.*AblationDays/yts;
    622 
    623                         /*Update array accordingly*/
    624                         smb[v]=MassBalanceMYearDebris;
    625                 }
    626 
    627                 /*Add input to element and Free memory*/
    628                 element->AddInput(SmbMassBalanceEnum,smb,P1Enum);
    629                 xDelete<IssmDouble>(surfacelist);
    630                 xDelete<IssmDouble>(smb);
    631                 xDelete<IssmDouble>(debriscover);
    632         }
     505void SmbDebrisEvattx(FemModel* femmodel){/*{{{*/
     506        for(Object* & object : femmodel->elements->objects){
     507                Element* element=xDynamicCast<Element*>(object);
     508                element->SmbDebrisEvatt();
     509        }
    633510}/*}}}*/
    634511void SmbGradientsComponentsx(FemModel* femmodel){/*{{{*/
    … …  
    641518}/*}}}*/
    642519#ifdef _HAVE_SEMIC_
    643 void SmbSemicx(FemModel* femmodel){/*{{{*/
    644 
    645         for(Object* & object : femmodel->elements->objects){
    646                 Element* element=xDynamicCast<Element*>(object);
    647                 element->SmbSemic();
     520void SmbSemicx(FemModel* femmodel,int ismethod){/*{{{*/
     521
     522        for(Object* & object : femmodel->elements->objects){
     523                Element* element=xDynamicCast<Element*>(object);
     524                if (ismethod == 1) element->SmbSemicTransient(); // Inwoo's version.
     525                else element->SmbSemic(); // original SmbSEMIC
    648526        }
    649527

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

    r27347 r28013  
    2323void SmbMeltComponentsx(FemModel* femmodel);
    2424void SmbGradientsComponentsx(FemModel* femmodel);
    25 void SmbDebrisMLx(FemModel* femmodel);
     25void SmbDebrisEvattx(FemModel* femmodel);
    2626/* SEMIC: */
    27 void SmbSemicx(FemModel* femmodel);
     27void SmbSemicx(FemModel* femmodel, int ismethod);
    2828/*GEMB: */
    2929void       Gembx(FemModel* femmodel);
Note: See TracChangeset for help on using the changeset viewer.