Changeset 23317


Ignore:
Timestamp:
09/19/18 11:50:04 (7 years ago)
Author:
rueckamp
Message:

NEW: SICOPOLIS PDD scheme

Location:
issm/trunk-jpl
Files:
4 added
13 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/Makefile.am

    r23252 r23317  
    175175                                        ./shared/Elements/PrintArrays.cpp\
    176176                                        ./shared/Elements/PddSurfaceMassBalance.cpp\
     177                                        ./shared/Elements/PddSurfaceMassBalanceSicopolis.cpp\
    177178                                        ./shared/Elements/ComputeDelta18oTemperaturePrecipitation.cpp\
    178179                                        ./shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp\

  • issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

    r23066 r23317  
    4343                int smb_model;
    4444                iomodel->FindConstant(&smb_model,"md.smb.model");
    45                 if(smb_model==SMBpddEnum)     isdynamic=true;
    46                 if(smb_model==SMBd18opddEnum) isdynamic=true;
     45                if(smb_model==SMBpddEnum)                               isdynamic=true;
     46                if(smb_model==SMBd18opddEnum)                   isdynamic=true;
     47                if(smb_model==SMBpddSicopolisEnum)      isdynamic=true;
    4748        }
    4849
    … …  
    12001201                _assert_((Hc+Ht)>0.);
    12011202                lambda = Hc/(Hc+Ht);
     1203                _assert_(lambda>=0.);
     1204                _assert_(lambda<=1.);
    12021205                kappa  = kappa_c*kappa_t/(lambda*kappa_t+(1.-lambda)*kappa_c); // ==(lambda/kappa_c + (1.-lambda)/kappa_t)^-1
    12031206        }       

  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r23066 r23317  
    55#include "../modules/modules.h"
    66
     7// FIX
     8#include "./shared/io/Print/Print.h"
     9
    710/*Model processing*/
    811void SmbAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    … …  
    2124
    2225        int    smb_model;
    23         bool   isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled;
     26        bool   isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled,isfirnwarming;
    2427
    2528        /*Update elements: */
    … …  
    8891                        }
    8992                        break;
     93                case SMBpddSicopolisEnum:
     94                        iomodel->FetchDataToInput(elements,"md.smb.s0p",SmbS0pEnum);
     95                        iomodel->FetchDataToInput(elements,"md.smb.s0t",SmbS0tEnum);
     96                        iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
     97                        iomodel->FetchDataToInput(elements,"md.smb.smb_corr",SmbSmbCorrEnum);
     98                        iomodel->FetchDataToInput(elements,"md.smb.precipitation",SmbPrecipitationEnum);
     99                        iomodel->FetchDataToInput(elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
     100                        iomodel->FetchDataToInput(elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
     101                        iomodel->FetchDataToInput(elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
     102                        break;
    90103                case SMBd18opddEnum:
    91104                        iomodel->FindConstant(&istemperaturescaled,"md.smb.istemperaturescaled");
    … …  
    148161        int     numoutputs;
    149162        char**  requestedoutputs = NULL;
    150         bool    isdelta18o,ismungsm,isd18opd,issetpddfac,interp;
     163        bool    isdelta18o,ismungsm,isd18opd,issetpddfac,interp,isfirnwarming;
    151164        int     smb_model;
    152165        IssmDouble *temp = NULL;
    … …  
    215228                        }
    216229                        break;
     230                case SMBpddSicopolisEnum:
     231                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.isfirnwarming",SmbIssetpddfacEnum));
     232                        break;
    217233                case SMBd18opddEnum:
    218234                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismungsm",SmbIsmungsmEnum));
    … …  
    288304                        PositiveDegreeDayx(femmodel);
    289305                        break;
     306                case SMBpddSicopolisEnum:
     307                        if(VerboseSolution()) _printf0_("   call SICOPOLIS positive degree day module\n");
     308                        PositiveDegreeDaySicopolisx(femmodel);
     309                        break;
    290310                case SMBd18opddEnum:
    291311                        bool isd18opd;
    … …  
    296316                                if(VerboseSolution()) _printf0_("   call positive degree day module\n");
    297317                                PositiveDegreeDayx(femmodel);
    298                         } 
     318                        }
    299319                        break;
    300320                case SMBgradientsEnum:

  • issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp

    r22874 r23317  
    3333                int smb_model;
    3434                iomodel->FindConstant(&smb_model,"md.smb.model");
    35                 if(smb_model==SMBpddEnum) isdynamic=true;
    36                 if(smb_model==SMBd18opddEnum) isdynamic=true;
     35                if(smb_model==SMBpddEnum)                               isdynamic=true;
     36                if(smb_model==SMBd18opddEnum)                   isdynamic=true;
     37                if(smb_model==SMBpddSicopolisEnum)  isdynamic=true;
    3738        }
    3839        else{

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r23242 r23317  
    25592559        xDelete<IssmDouble>(tmp);
    25602560
     2561}
     2562/*}}}*/
     2563void       Element::PositiveDegreeDaySicopolis(bool isfirnwarming){/*{{{*/
     2564
     2565        /* General FIXMEs: get Tmelting point, pddicefactor, pddsnowfactor, sigma from parameters/user input */
     2566
     2567        int  numvertices = this->GetNumberOfVertices();
     2568
     2569        int        i;
     2570        IssmDouble* smb=xNew<IssmDouble>(numvertices);          // surface mass balance
     2571        IssmDouble* melt=xNew<IssmDouble>(numvertices);         // melting comp. of surface mass balance
     2572        IssmDouble* accu=xNew<IssmDouble>(numvertices);         // accuumulation comp. of surface mass balance
     2573        IssmDouble* melt_star=xNew<IssmDouble>(numvertices);
     2574        IssmDouble* monthlytemperatures=xNew<IssmDouble>(12*numvertices);
     2575        IssmDouble* monthlyprec=xNew<IssmDouble>(12*numvertices);
     2576        IssmDouble* yearlytemperatures=xNew<IssmDouble>(numvertices); memset(yearlytemperatures, 0., numvertices*sizeof(IssmDouble));
     2577        IssmDouble* s=xNew<IssmDouble>(numvertices);                    // actual surface height
     2578        IssmDouble* s0p=xNew<IssmDouble>(numvertices);          // reference elevation for precip.
     2579        IssmDouble* s0t=xNew<IssmDouble>(numvertices);          // reference elevation for temperature
     2580        IssmDouble* smbcorr=xNew<IssmDouble>(numvertices); // surface mass balance correction; will be added after pdd call
     2581        IssmDouble* p_ampl=xNew<IssmDouble>(numvertices);       // precip anomaly
     2582        IssmDouble* t_ampl=xNew<IssmDouble>(numvertices);       // remperature anomaly
     2583        IssmDouble rho_water,rho_ice,desfac,rlaps;
     2584        IssmDouble inv_twelve=1./12.;                                                           //factor for monthly average
     2585        IssmDouble time,yts,time_yr;
     2586
     2587        /*Get material parameters :*/
     2588        rho_water=this->matpar->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     2589        rho_ice=this->matpar->GetMaterialParameter(MaterialsRhoIceEnum);
     2590       
     2591        /*Get parameters for height corrections*/
     2592        desfac=this->matpar->GetMaterialParameter(SmbDesfacEnum);
     2593        rlaps=this->matpar->GetMaterialParameter(SmbRlapsEnum);
     2594       
     2595        /*Recover monthly temperatures and precipitation*/
     2596        Input*     input=this->inputs->GetInput(SmbMonthlytemperaturesEnum); _assert_(input);
     2597        Input*     input2=this->inputs->GetInput(SmbPrecipitationEnum); _assert_(input2);
     2598        /*Recover smb correction term */
     2599        Input*     input3=this->inputs->GetInput(SmbSmbCorrEnum); _assert_(input3);
     2600       
     2601        /* Get time */
     2602        this->parameters->FindParam(&time,TimeEnum);
     2603        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     2604        time_yr=floor(time/yts)*yts;
     2605
     2606        /* Set parameters for finrnwarming */
     2607        IssmDouble MU_0         = 9.7155; //Firn-warming correction, in (d*deg C)/(mm WE)
     2608        IssmDouble mu           = MU_0*(1000.0*86400.0)*(rho_ice/rho_water);   // (d*deg C)/(mm WE) --> (s*deg C)/(m IE)
     2609       
     2610        /*loop over vertices: */
     2611        Gauss* gauss=this->NewGauss();
     2612        for(int month=0;month<12;month++){
     2613                for(int iv=0;iv<numvertices;iv++){
     2614                        gauss->GaussVertex(iv);
     2615                        input->GetInputValue(&monthlytemperatures[iv*12+month],gauss,(month+1)/12.*yts);
     2616                        monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module
     2617                        input2->GetInputValue(&monthlyprec[iv*12+month],gauss,(month+1)/12.*yts);
     2618                        monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts;
     2619                }
     2620        }
     2621
     2622        /*Recover info at the vertices: */
     2623        GetInputListOnVertices(&s[0],SurfaceEnum);
     2624        GetInputListOnVertices(&s0p[0],SmbS0pEnum);
     2625        GetInputListOnVertices(&s0t[0],SmbS0tEnum);
     2626        GetInputListOnVertices(&smbcorr[0],SmbSmbCorrEnum);
     2627        GetInputListOnVertices(&t_ampl[0],SmbTemperaturesAnomalyEnum);
     2628        GetInputListOnVertices(&p_ampl[0],SmbPrecipitationsAnomalyEnum);
     2629
     2630        /*measure the surface mass balance*/
     2631        for (int iv = 0; iv<numvertices; iv++){
     2632                smb[iv]=PddSurfaceMassBalanceSicopolis(&monthlytemperatures[iv*12], &monthlyprec[iv*12],
     2633                                        &melt[iv], &accu[iv], &melt_star[iv], &t_ampl[iv], &p_ampl[iv], yts, s[iv],
     2634                                        desfac, s0t[iv], s0p[iv],rlaps,rho_water,rho_ice);
     2635
     2636                /* make correction */
     2637                smb[iv] = smb[iv]+smbcorr[iv];
     2638                /*Get yearlytemperatures */
     2639                for(int month=0;month<12;month++) yearlytemperatures[iv]=yearlytemperatures[iv]+((monthlytemperatures[iv*12+month]+273.15)+t_ampl[iv])*inv_twelve; // Has to be in Kelvin
     2640       
     2641                if(isfirnwarming){
     2642                        if(melt_star[iv]>=melt[iv]){
     2643                                yearlytemperatures[iv]= yearlytemperatures[iv]+mu*(melt_star[iv]-melt[iv]);
     2644                        }
     2645                        else{
     2646                                yearlytemperatures[iv]= yearlytemperatures[iv];
     2647                        }
     2648                }
     2649                if (yearlytemperatures[iv]>273.15) yearlytemperatures[iv]=273.15;
     2650        }
     2651
     2652        switch(this->ObjectEnum()){
     2653                case TriaEnum:
     2654                        // this->inputs->AddInput(new TriaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
     2655                        this->inputs->AddInput(new TriaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
     2656                        this->inputs->AddInput(new TriaInput(SmbMassBalanceEnum,&smb[0],P1Enum));
     2657                        this->inputs->AddInput(new TriaInput(SmbAccumulationEnum,&accu[0],P1Enum));
     2658                        this->inputs->AddInput(new TriaInput(SmbMeltEnum,&melt[0],P1Enum));
     2659                        break;
     2660                case PentaEnum:
     2661                        bool isthermal;
     2662                        this->parameters->FindParam(&isthermal,TransientIsthermalEnum);
     2663                        if(isthermal){
     2664                                bool isenthalpy;
     2665                                this->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
     2666                                if(IsOnSurface()){
     2667                                        /*Here, we want to change the BC of the thermal model, keep
     2668                                         * the temperatures as they are for the base of the penta and
     2669                                         * use yearlytemperatures for the top*/
     2670                                        GetInputListOnVertices(&s[0],TemperatureEnum);
     2671                                        yearlytemperatures[0] = s[0];
     2672                                        yearlytemperatures[1] = s[1];
     2673                                        yearlytemperatures[2] = s[2];
     2674                                        this->inputs->AddInput(new PentaInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
     2675                                        if(isenthalpy){
     2676                                                /*Convert that to enthalpy for the enthalpy model*/
     2677                                                IssmDouble enthalpy[6];
     2678                                                GetInputListOnVertices(&enthalpy[0],EnthalpyEnum);
     2679                                                ThermalToEnthalpy(&enthalpy[3],yearlytemperatures[3],0.,0.);
     2680                                                ThermalToEnthalpy(&enthalpy[4],yearlytemperatures[4],0.,0.);
     2681                                                ThermalToEnthalpy(&enthalpy[5],yearlytemperatures[5],0.,0.);
     2682                                                this->inputs->AddInput(new PentaInput(EnthalpyEnum,&enthalpy[0],P1Enum));
     2683                                        }
     2684                                }
     2685                        }
     2686                        this->inputs->AddInput(new PentaInput(SmbMassBalanceEnum,&smb[0],P1Enum));
     2687                        this->inputs->AddInput(new PentaInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
     2688                        this->inputs->AddInput(new PentaInput(SmbAccumulationEnum,&accu[0],P1Enum));
     2689                        this->inputs->AddInput(new PentaInput(SmbMeltEnum,&melt[0],P1Enum));
     2690                        this->InputExtrude(TemperaturePDDEnum,-1);
     2691                        this->InputExtrude(SmbMassBalanceEnum,-1);
     2692                        this->InputExtrude(SmbAccumulationEnum,-1);
     2693                        this->InputExtrude(SmbMeltEnum,-1);
     2694                        break;
     2695                case TetraEnum:
     2696                        if(IsOnSurface()){
     2697                                GetInputListOnVertices(&s[0],TemperatureEnum);
     2698                                yearlytemperatures[0] = s[0];
     2699                                yearlytemperatures[1] = s[1];
     2700                                yearlytemperatures[2] = s[2];
     2701                                this->inputs->AddInput(new TetraInput(TemperatureEnum,&yearlytemperatures[0],P1Enum));
     2702                        }
     2703                        this->inputs->AddInput(new TetraInput(SmbMassBalanceEnum,&smb[0],P1Enum));
     2704                        this->inputs->AddInput(new TetraInput(TemperaturePDDEnum,&yearlytemperatures[0],P1Enum));
     2705                        this->InputExtrude(TemperaturePDDEnum,-1);
     2706                        this->InputExtrude(SmbMassBalanceEnum,-1);
     2707                        break;
     2708                default: _error_("Not implemented yet");
     2709        }
     2710
     2711        /*clean-up*/
     2712        delete gauss;
     2713        xDelete<IssmDouble>(monthlytemperatures);
     2714        xDelete<IssmDouble>(monthlyprec);
     2715        xDelete<IssmDouble>(smb);
     2716        xDelete<IssmDouble>(melt);
     2717        xDelete<IssmDouble>(accu);
     2718        xDelete<IssmDouble>(yearlytemperatures);
     2719        xDelete<IssmDouble>(s);
     2720        xDelete<IssmDouble>(s0t);
     2721        xDelete<IssmDouble>(s0p);
     2722        xDelete<IssmDouble>(t_ampl);
     2723        xDelete<IssmDouble>(p_ampl);
     2724        xDelete<IssmDouble>(smbcorr);
     2725        xDelete<IssmDouble>(melt_star);
    25612726}
    25622727/*}}}*/

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r23058 r23317  
    145145                ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
    146146                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
     147                void               PositiveDegreeDaySicopolis(bool isfirnwarming);
    147148                IssmDouble         PureIceEnthalpy(IssmDouble pressure);
    148149                void               ResultInterpolation(int* pinterpolation,int*nodesperelement,int* parray_size, int output_enum);

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp

    r23093 r23317  
    106106                                        iomodel->FindConstant(&this->rlaps,"md.smb.rlaps");
    107107                                        iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm");
     108                                        break;
     109                                case SMBpddSicopolisEnum:
     110                                        iomodel->FindConstant(&this->desfac,"md.smb.desfac");
     111                                        iomodel->FindConstant(&this->rlaps,"md.smb.rlaps");
    108112                                        break;
    109113                                case SMBd18opddEnum:
    … …  
    210214                                                                iomodel->FindConstant(&this->rlaps,"md.smb.rlaps");
    211215                                                                iomodel->FindConstant(&this->rlapslgm,"md.smb.rlapslgm");
     216                                                                break;
     217                                                        case SMBpddSicopolisEnum:
     218                                                                iomodel->FindConstant(&this->desfac,"md.smb.desfac");
     219                                                                iomodel->FindConstant(&this->rlaps,"md.smb.rlaps");
    212220                                                                break;
    213221                                                        case SMBd18opddEnum:

  • issm/trunk-jpl/src/c/cores/transient_core.cpp

    r23252 r23317  
    371371                                femmodel->parameters->FindParam(&smb_model,SmbEnum);
    372372                                if(isenthalpy){
    373                                         if(smb_model==SMBpddEnum)     ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
    374                                         if(smb_model==SMBd18opddEnum) ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
     373                                        if(smb_model==SMBpddEnum)                               ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
     374                                        if(smb_model==SMBd18opddEnum)                   ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
     375                                        if(smb_model==SMBpddSicopolisEnum)      ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
    375376                                }
    376377                                else{
    377                                         if(smb_model==SMBpddEnum)     ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
    378                                         if(smb_model==SMBd18opddEnum) ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
     378                                        if(smb_model==SMBpddEnum)                               ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
     379                                        if(smb_model==SMBd18opddEnum)                   ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
     380                                        if(smb_model==SMBpddSicopolisEnum)      ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
    379381                                }
    380382                        }

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

    r23066 r23317  
    262262        xDelete<IssmDouble>(pds);
    263263}/*}}}*/
     264void PositiveDegreeDaySicopolisx(FemModel* femmodel){/*{{{*/
     265       
     266        bool isfirnwarming;
     267        femmodel->parameters->FindParam(&isfirnwarming,SmbIssetpddfacEnum);
     268       
     269        for(int i=0;i<femmodel->elements->Size();i++){
     270                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     271                element->PositiveDegreeDaySicopolis(isfirnwarming);
     272        }
     273
     274}/*}}}*/
    264275void SmbHenningx(FemModel* femmodel){/*{{{*/
    265276

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

    r22482 r23317  
    1616void Delta18opdParameterizationx(FemModel* femmodel);
    1717void PositiveDegreeDayx(FemModel* femmodel);
     18void PositiveDegreeDaySicopolisx(FemModel* femmodel);
    1819void SmbHenningx(FemModel* femmodel);
    1920void SmbComponentsx(FemModel* femmodel);

  • issm/trunk-jpl/src/c/shared/Elements/elements.h

    r23255 r23317  
    2424                                 IssmDouble TdiffTime,IssmDouble sealevTime,IssmDouble pddsnowfac,IssmDouble pddicefac,
    2525                                 IssmDouble rho_water, IssmDouble rho_ice);
     26IssmDouble PddSurfaceMassBalanceSicopolis(IssmDouble* monthlytemperatures,  IssmDouble* monthlyprec,
     27                                 IssmDouble* melt, IssmDouble* accu, IssmDouble* melt_star, IssmDouble* t_ampl, IssmDouble* p_ampl,
     28                                 IssmDouble yts, IssmDouble s, IssmDouble desfac,IssmDouble s0t,
     29                                 IssmDouble s0p, IssmDouble rlaps, IssmDouble rho_water, IssmDouble rho_ice);
    2630void ComputeDelta18oTemperaturePrecipitation(IssmDouble Delta18oSurfacePresent, IssmDouble Delta18oSurfaceLgm, IssmDouble Delta18oSurfaceTime,
    2731                                             IssmDouble Delta18oPresent, IssmDouble Delta18oLgm, IssmDouble Delta18oTime,

  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r23242 r23317  
    11571157        XYEnum,
    11581158        XYZEnum,
     1159        SMBpddSicopolisEnum,
     1160        SmbIsfirnwarmingEnum,
     1161        SmbSmbCorrEnum,
     1162        SmbTemperaturesAnomalyEnum,
     1163        SmbPrecipitationsAnomalyEnum,
    11591164        /*}}}*/
    11601165        /*Unused?{{{*/

  • issm/trunk-jpl/src/c/shared/io/Marshalling/IoCodeConversions.cpp

    r23093 r23317  
    173173                case 8: return SMBgembEnum;
    174174                case 9: return SMBgradientselaEnum;
     175                case 10: return SMBpddSicopolisEnum;
    175176                default: _error_("Marshalled SMB code \""<<enum_in<<"\" not supported yet");
    176177        }
Note: See TracChangeset for help on using the changeset viewer.