Changeset 27498


Ignore:
Timestamp:
01/04/23 16:20:32 (2 years ago)
Author:
inwoo
Message:

CHG: update SEMIC - enable transient mode.

Location:
issm/trunk-jpl/src
Files:
1 added
13 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk-jpl/src/c/Makefile.am

    r27462 r27498  
    601601if FORTRAN
    602602issm_sources += ./modules/SurfaceMassBalancex/run_semic.f90
     603issm_sources += ./modules/SurfaceMassBalancex/run_semic_transient.f90
    603604endif
    604605endif
  • issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

    r27323 r27498  
    206206                        break;
    207207                case SMBsemicEnum:
    208                         iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
     208                        int ismethod;
     209                        //if(VerboseSolution()) _printf0_("   smb semic: UpdateElements.\n");
     210                        //iomodel->FetchDataToInput(inputs,elements,"md.thermal.spctemperature",ThermalSpctemperatureEnum);
    209211                        iomodel->FetchDataToInput(inputs,elements,"md.smb.s0gcm",SmbS0gcmEnum);
    210212                        iomodel->FetchDataToInput(inputs,elements,"md.smb.dailysnowfall",SmbDailysnowfallEnum);
     
    217219                        iomodel->FetchDataToInput(inputs,elements,"md.smb.dailyairhumidity",SmbDailyairhumidityEnum);
    218220                        iomodel->FetchDataToInput(inputs,elements,"md.smb.dailytemperature",SmbDailytemperatureEnum);
     221                        // assign initial SEMIC temperature from initialization class.
     222                        if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - temperature.\n");
     223                        iomodel->FetchDataToInput(inputs,elements,"md.initialization.temperature",TemperatureSEMICEnum);
     224
     225                        iomodel->FindConstant(&ismethod,"md.smb.ismethod");
     226                        if (ismethod == 1){
     227                                // initial albedo
     228                                if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - albedo.\n");
     229                                iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo",SmbAlbedoEnum);
     230                                iomodel->FetchDataToInput(inputs,elements,"md.smb.albedo_snow",SmbAlbedoSnowEnum);
     231                                if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - Hice/Hsnow.\n");
     232                                iomodel->FetchDataToInput(inputs,elements,"md.smb.hice",SmbHIceEnum);
     233                                iomodel->FetchDataToInput(inputs,elements,"md.smb.hsnow",SmbHSnowEnum);
     234
     235                                // initial Temperature amplitude.
     236                                if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - Tamp.\n");
     237                                iomodel->FetchDataToInput(inputs,elements,"md.smb.Tamp",SmbTampEnum);
     238
     239                                // assign masking
     240                                iomodel->FetchDataToInput(inputs,elements,"md.smb.mask",SmbMaskEnum);
     241                                if(VerboseSolution()) _printf0_("   smb semic: UpdateElements - done.\n");
     242                        }
    219243                        break;
    220244                case SMBdebrisMLEnum:
     
    420444                        break;
    421445                case SMBsemicEnum:
     446                        int ismethod;
     447                        parameters->FindParam(&ismethod,SmbSemicMethodEnum);
     448                        if (ismethod == 1){
     449                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.desfacElevation",SmbDesfacElevEnum));
     450                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.hcrit",SmbSemicHcritEnum));
     451                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.rcrit",SmbSemicRcritEnum));
     452                                /*Define albedo parameters.*/
     453                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.albedo_scheme",SmbAlbedoSchemeEnum));
     454                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.alb_smax",SmbAlbedoSnowMaxEnum));
     455                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.alb_smin",SmbAlbedoSnowMinEnum));
     456                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.albi",SmbAlbedoIceEnum));
     457                                parameters->AddObject(iomodel->CopyConstantObject("md.smb.albl",SmbAlbedoLandEnum));
     458                        }
    422459                        /*Nothing to add to parameters*/
    423460                        break;
     
    514551                case SMBsemicEnum:
    515552                        #ifdef _HAVE_SEMIC_
    516                         if(VerboseSolution())_printf0_("  call smb SEMIC module\n");
    517                         SmbSemicx(femmodel);
     553                        if(VerboseSolution())_printf0_("   call smb SEMIC module\n");
     554                        int ismethod;
     555                        femmodel->parameters->FindParam(&ismethod,SmbSemicMethodEnum);
     556                        SmbSemicx(femmodel,ismethod);
    518557                        #else
    519558                        _error_("SEMIC not installed");
  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r27491 r27498  
    3030extern "C" void run_semic_(IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in,
    3131                        IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *saccu_out, IssmDouble *smelt_out);
     32extern "C" void run_semic_transient_(IssmDouble *sf_in, IssmDouble *rf_in, IssmDouble *swd_in, IssmDouble *lwd_in, IssmDouble *wind_in, IssmDouble *sp_in, IssmDouble *rhoa_in,
     33                        IssmDouble *qq_in, IssmDouble *tt_in, IssmDouble *tsurf_in, IssmDouble *tstic,
     34                        IssmDouble *hcrit, IssmDouble *rcrit,
     35                        IssmDouble *mask, IssmDouble *hice, IssmDouble *hsnow,
     36                        IssmDouble *albedo, IssmDouble *albedo_snow, int *alb_scheme, IssmDouble *alb_smax, IssmDouble *alb_smin, IssmDouble *albi, IssmDouble *albl,
     37                        IssmDouble *Tamp,
     38                        IssmDouble *tsurf_out, IssmDouble *smb_out, IssmDouble *saccu_out, IssmDouble *smelt_out, IssmDouble *albedo_out, IssmDouble *albedo_snow_out, IssmDouble *hsnow_out, IssmDouble *hice_out);
    3239#endif
    3340// _HAVE_SEMIC_
     
    42084215        xDelete<IssmDouble>(saccu_out);
    42094216        xDelete<IssmDouble>(tsurf_out);
     4217        xDelete<IssmDouble>(s);
     4218        xDelete<IssmDouble>(st);
     4219        xDelete<IssmDouble>(s0gcm);
     4220}
     4221/*}}}*/
     4222void       Element::SmbSemicTransient(){/*{{{*/
     4223
     4224        /*only compute SMB at the surface: */
     4225        if (!IsOnSurface()) return;
     4226
     4227        const int NUM_VERTICES                                  = this->GetNumberOfVertices();
     4228        const int NLOOP = 10; // default internal iteration for SEMIC (Ruckamp et al. 2018).
     4229
     4230        IssmDouble* s=xNew<IssmDouble>(NUM_VERTICES);
     4231        IssmDouble* s0gcm=xNew<IssmDouble>(NUM_VERTICES);
     4232        IssmDouble* st=xNew<IssmDouble>(NUM_VERTICES);
     4233
     4234        // daily forcing inputs
     4235        IssmDouble* dailyrainfall   =xNew<IssmDouble>(NUM_VERTICES);
     4236        IssmDouble* dailysnowfall   =xNew<IssmDouble>(NUM_VERTICES);
     4237        IssmDouble* dailydlradiation=xNew<IssmDouble>(NUM_VERTICES);
     4238        IssmDouble* dailydsradiation=xNew<IssmDouble>(NUM_VERTICES);
     4239        IssmDouble* dailywindspeed  =xNew<IssmDouble>(NUM_VERTICES);
     4240        IssmDouble* dailypressure   =xNew<IssmDouble>(NUM_VERTICES);
     4241        IssmDouble* dailyairdensity =xNew<IssmDouble>(NUM_VERTICES);
     4242        IssmDouble* dailyairhumidity=xNew<IssmDouble>(NUM_VERTICES);
     4243        IssmDouble* dailytemperature=xNew<IssmDouble>(NUM_VERTICES);
     4244
     4245        IssmDouble* tsurf_in        =xNew<IssmDouble>(NUM_VERTICES);
     4246        IssmDouble* mask_in         =xNew<IssmDouble>(NUM_VERTICES);
     4247        IssmDouble* Tamp_in         =xNew<IssmDouble>(NUM_VERTICES);
     4248        IssmDouble* albedo_in       =xNew<IssmDouble>(NUM_VERTICES);
     4249        IssmDouble* albedo_snow_in  =xNew<IssmDouble>(NUM_VERTICES);
     4250        IssmDouble* hice_in         =xNew<IssmDouble>(NUM_VERTICES);
     4251        IssmDouble* hsnow_in        =xNew<IssmDouble>(NUM_VERTICES);
     4252
     4253        // daily outputs
     4254        IssmDouble* tsurf_out  =xNew<IssmDouble>(NUM_VERTICES); memset(tsurf_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4255        IssmDouble* smb_out    =xNew<IssmDouble>(NUM_VERTICES); memset(smb_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4256        IssmDouble* saccu_out  =xNew<IssmDouble>(NUM_VERTICES); memset(saccu_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4257        IssmDouble* smelt_out  =xNew<IssmDouble>(NUM_VERTICES); memset(smelt_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4258        IssmDouble* albedo_out =xNew<IssmDouble>(NUM_VERTICES); memset(albedo_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4259        IssmDouble* albedo_snow_out =xNew<IssmDouble>(NUM_VERTICES); memset(albedo_snow_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4260        IssmDouble* hsnow_out    =xNew<IssmDouble>(NUM_VERTICES); memset(hsnow_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4261        IssmDouble* hice_out    =xNew<IssmDouble>(NUM_VERTICES); memset(hice_out, 0., NUM_VERTICES*sizeof(IssmDouble));
     4262
     4263        IssmDouble rho_water,rho_ice,desfac,desfacElev,rlaps,rdl;
     4264        IssmDouble alb_smax, alb_smin, albi, albl;
     4265        int alb_scheme;
     4266        IssmDouble hcrit, rcrit;
     4267
     4268        IssmDouble time,yts,time_yr,dt;
     4269
     4270        /* Get time: */
     4271        this->parameters->FindParam(&time,TimeEnum);
     4272        this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
     4273        this->parameters->FindParam(&yts,ConstantsYtsEnum);
     4274        time_yr=floor(time/yts)*yts;
     4275        //dt = dt * yts;
     4276
     4277        /*Get material parameters :*/
     4278        rho_water=this->FindParam(MaterialsRhoSeawaterEnum);
     4279        rho_ice=this->FindParam(MaterialsRhoIceEnum);
     4280        desfac=this->FindParam(SmbDesfacEnum);
     4281        desfacElev=this->FindParam(SmbDesfacElevEnum);
     4282        rlaps=this->FindParam(SmbRlapsEnum);
     4283        rdl=this->FindParam(SmbRdlEnum);
     4284
     4285        this->FindParam(&alb_scheme,SmbAlbedoSchemeEnum);
     4286        this->FindParam(&hcrit,SmbSemicHcritEnum);
     4287        this->FindParam(&rcrit,SmbSemicRcritEnum);
     4288        alb_smax=this->FindParam(SmbAlbedoSnowMaxEnum);
     4289        alb_smin=this->FindParam(SmbAlbedoSnowMinEnum);
     4290        albi=this->FindParam(SmbAlbedoIceEnum);
     4291        albl=this->FindParam(SmbAlbedoLandEnum);
     4292
     4293        /* Recover info at the vertices: */
     4294        GetInputListOnVertices(&s[0],SurfaceEnum);
     4295        GetInputListOnVertices(&s0gcm[0],SmbS0gcmEnum);
     4296
     4297        if(VerboseSmb() && this->Sid()==0){
     4298                _printf0_("smb core: allocate inputs.\n");
     4299                _printf0_("smb core: time_yr  : " << time_yr/yts <<"\n");
     4300                _printf0_("smb core: time     : " << time <<"\n");
     4301                _printf0_("smb core: dt       : " << dt <<"\n");
     4302        }
     4303        /* loop over vertices and days */ //FIXME account for leap years (365 -> 366)
     4304        Gauss* gauss=this->NewGauss();
     4305        /* Retrieve inputs: */
     4306        Input* dailysnowfall_input    = this->GetInput(SmbDailysnowfallEnum,time); _assert_(dailysnowfall_input);
     4307        Input* dailyrainfall_input    = this->GetInput(SmbDailyrainfallEnum,time); _assert_(dailyrainfall_input);
     4308        Input* dailydlradiation_input = this->GetInput(SmbDailydlradiationEnum,time); _assert_(dailydlradiation_input);
     4309        Input* dailydsradiation_input = this->GetInput(SmbDailydsradiationEnum,time); _assert_(dailydsradiation_input);
     4310        Input* dailywindspeed_input   = this->GetInput(SmbDailywindspeedEnum,time); _assert_(dailywindspeed_input);
     4311        Input* dailypressure_input    = this->GetInput(SmbDailypressureEnum,time); _assert_(dailypressure_input);
     4312        Input* dailyairdensity_input  = this->GetInput(SmbDailyairdensityEnum,time); _assert_(dailyairdensity_input);
     4313        Input* dailyairhumidity_input = this->GetInput(SmbDailyairhumidityEnum,time); _assert_(dailyairhumidity_input);
     4314        Input* dailytemperature_input = this->GetInput(SmbDailytemperatureEnum,time); _assert_(dailytemperature_input);
     4315
     4316        Input* tsurf_input            = this->GetInput(TemperatureSEMICEnum); _assert_(tsurf_in);
     4317        Input* mask_input        = this->GetInput(SmbMaskEnum); _assert_(mask_input);
     4318        Input* Tamp_input        = this->GetInput(SmbTampEnum); _assert_(Tamp_input);
     4319        Input* albedo_input      = this->GetInput(SmbAlbedoEnum); _assert_(albedo_input);
     4320        Input* albedo_snow_input = this->GetInput(SmbAlbedoSnowEnum); _assert_(albedo_snow_input);
     4321        Input* hice_input        = this->GetInput(SmbHIceEnum); _assert_(hice_input);
     4322        Input* hsnow_input       = this->GetInput(SmbHSnowEnum); _assert_(hsnow_input);
     4323
     4324        if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: assign inputs.\n");
     4325        for(int iv=0;iv<NUM_VERTICES;iv++){
     4326                gauss->GaussVertex(iv);
     4327                /* get forcing */
     4328                dailyrainfall_input->GetInputValue(&dailyrainfall[iv],gauss);
     4329                dailysnowfall_input->GetInputValue(&dailysnowfall[iv],gauss);
     4330                dailydlradiation_input->GetInputValue(&dailydlradiation[iv],gauss);
     4331                dailydsradiation_input->GetInputValue(&dailydsradiation[iv],gauss);
     4332                dailywindspeed_input->GetInputValue(&dailywindspeed[iv],gauss);
     4333                dailypressure_input->GetInputValue(&dailypressure[iv],gauss);
     4334                dailyairdensity_input->GetInputValue(&dailyairdensity[iv],gauss);
     4335                dailyairhumidity_input->GetInputValue(&dailyairhumidity[iv],gauss);
     4336                dailytemperature_input->GetInputValue(&dailytemperature[iv],gauss);
     4337                tsurf_input->GetInputValue(&tsurf_in[iv],gauss);
     4338
     4339                /* Get Albedo information */
     4340                mask_input->GetInputValue(&mask_in[iv],gauss);
     4341                Tamp_input->GetInputValue(&Tamp_in[iv],gauss);
     4342                albedo_input->GetInputValue(&albedo_in[iv],gauss);
     4343                albedo_snow_input->GetInputValue(&albedo_snow_in[iv],gauss);
     4344
     4345                hsnow_input->GetInputValue(&hsnow_in[iv],gauss);
     4346                hice_input->GetInputValue(&hice_in[iv],gauss);
     4347
     4348                /* Surface temperature correction */
     4349                st[iv]=(s[iv]-s0gcm[iv])/1000.;
     4350                dailytemperature[iv]=dailytemperature[iv]-rlaps *st[iv];
     4351
     4352                /* Precipitation correction (Vizcaino et al. 2010) */
     4353                if (s0gcm[iv] < desfacElev) {
     4354                        dailysnowfall[iv] = dailysnowfall[iv]*exp(desfac*(max(s[iv],desfacElev)-desfacElev));
     4355                        dailyrainfall[iv] = dailyrainfall[iv]*exp(desfac*(max(s[iv],desfacElev)-desfacElev));
     4356                }else{
     4357                        dailysnowfall[iv] = dailysnowfall[iv]*exp(desfac*(max(s[iv],desfacElev)-s0gcm[iv]));
     4358                        dailyrainfall[iv] = dailyrainfall[iv]*exp(desfac*(max(s[iv],desfacElev)-s0gcm[iv]));
     4359                }
     4360
     4361                /* downward longwave radiation correction (Marty et al. 2002) */
     4362                st[iv]=(s[iv]-s0gcm[iv])/1000.;
     4363                dailydlradiation[iv]=dailydlradiation[iv]+rdl*st[iv];
     4364        }
     4365        if(VerboseSmb() && this->Sid()==0){
     4366                _printf0_("smb core: assign tsurf_in        :" << tsurf_in[0] << "\n");
     4367                _printf0_("smb core: assign dailytemperature:" << dailytemperature[0] << "\n");
     4368                _printf0_("smb core: assign hsnow           :" << hsnow_in[0] << "\n");
     4369                _printf0_("smb core: assign hice            :" << hice_in[0] << "\n");
     4370                _printf0_("smb core: assign mask            :" << mask_in[0] << "\n");
     4371                _printf0_("smb core: assign albedo          :" << albedo_in[0] << "\n");
     4372                _printf0_("smb core: assign albedo_snow     :" << albedo_snow_in[0] << "\n");
     4373        }
     4374
     4375        if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: call semic module.\n");
     4376        for (int iv = 0; iv<NUM_VERTICES; iv++){
     4377                /* call semic */
     4378                run_semic_transient_(&dailysnowfall[iv], &dailyrainfall[iv], &dailydsradiation[iv], &dailydlradiation[iv],
     4379                                &dailywindspeed[iv], &dailypressure[iv], &dailyairdensity[iv], &dailyairhumidity[iv], &dailytemperature[iv], &tsurf_in[iv], &dt,
     4380                                &hcrit, &rcrit,
     4381                                &mask_in[iv], &hice_in[iv], &hsnow_in[iv],
     4382                                &albedo_in[iv], &albedo_snow_in[iv], &alb_scheme, &alb_smax, &alb_smin, &albi, &albl,
     4383                                &Tamp_in[iv],
     4384                                &tsurf_out[iv], &smb_out[iv], &saccu_out[iv], &smelt_out[iv],&albedo_out[iv], &albedo_snow_out[iv], &hsnow_out[iv], &hice_out[iv]);
     4385        }
     4386
     4387        if(VerboseSmb() && this->Sid()==0){
     4388                _printf0_("smb core: tsurf_out " << tsurf_out[0] << "\n");
     4389        }
     4390
     4391        switch(this->ObjectEnum()){
     4392                case TriaEnum:
     4393                        //this->AddInput(TemperatureEnum,&tsurf_out[0],P1Enum); // TODO add TemperatureSEMICEnum to EnumDefinitions
     4394                        this->AddInput(TemperatureSEMICEnum,&tsurf_out[0],P1Enum); // TODO add TemperatureSEMICEnum to EnumDefinitions
     4395                        this->AddInput(SmbMassBalanceEnum,&smb_out[0],P1Enum);
     4396                        this->AddInput(SmbAccumulationEnum,&saccu_out[0],P1Enum);
     4397                        this->AddInput(SmbMeltEnum,&smelt_out[0],P1Enum);
     4398                        this->AddInput(SmbAlbedoEnum,&albedo_out[0],P1Enum);
     4399                        this->AddInput(SmbAlbedoSnowEnum,&albedo_snow_out[0],P1Enum);
     4400                        this->AddInput(SmbHSnowEnum,&hsnow_out[0],P1Enum);
     4401                        this->AddInput(SmbHIceEnum,&hice_out[0],P1Enum);
     4402                        break;
     4403                case PentaEnum:
     4404                        // TODO
     4405                        break;
     4406                case TetraEnum:
     4407                        // TODO
     4408                        break;
     4409                default: _error_("Not implemented yet");
     4410        }
     4411
     4412        /*clean-up*/
     4413        delete gauss;
     4414        xDelete<IssmDouble>(dailysnowfall);
     4415        xDelete<IssmDouble>(dailyrainfall);
     4416        xDelete<IssmDouble>(dailydlradiation);
     4417        xDelete<IssmDouble>(dailydsradiation);
     4418        xDelete<IssmDouble>(dailywindspeed);
     4419        xDelete<IssmDouble>(dailypressure);
     4420        xDelete<IssmDouble>(dailyairdensity);
     4421        xDelete<IssmDouble>(dailyairhumidity);
     4422        xDelete<IssmDouble>(dailypressure);
     4423        xDelete<IssmDouble>(dailytemperature);
     4424
     4425        xDelete<IssmDouble>(tsurf_out);
     4426        xDelete<IssmDouble>(smb_out);
     4427        xDelete<IssmDouble>(saccu_out);
     4428        xDelete<IssmDouble>(smelt_out);
     4429        xDelete<IssmDouble>(albedo_out);
     4430        xDelete<IssmDouble>(albedo_snow_out);
     4431
     4432        xDelete<IssmDouble>(hsnow_in);
     4433        xDelete<IssmDouble>(hice_in);
     4434        xDelete<IssmDouble>(mask_in);
     4435        xDelete<IssmDouble>(Tamp_in);
     4436        xDelete<IssmDouble>(albedo_snow_in);
     4437        xDelete<IssmDouble>(albedo_in);
     4438        xDelete<IssmDouble>(tsurf_in);
    42104439        xDelete<IssmDouble>(s);
    42114440        xDelete<IssmDouble>(st);
  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r27491 r27498  
    187187                void               SetIntInput(Inputs* inputs,int enum_in,int value);
    188188                void               SmbSemic();
     189                void               SmbSemicTransient();
    189190                int                Sid();
    190191                void               SmbGemb(IssmDouble timeinputs, int count);
  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r27462 r27498  
    487487                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.rlaps",SmbRlapsEnum));
    488488                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.rdl",SmbRdlEnum));
     489                        parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismethod",SmbSemicMethodEnum));
    489490                        break;
    490491                case SMBdebrisMLEnum:
  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp

    r27491 r27498  
    655655}/*}}}*/
    656656#ifdef _HAVE_SEMIC_
    657 void SmbSemicx(FemModel* femmodel){/*{{{*/
    658 
    659         for(Object* & object : femmodel->elements->objects){
    660                 Element* element=xDynamicCast<Element*>(object);
    661                 element->SmbSemic();
     657void SmbSemicx(FemModel* femmodel,int ismethod){/*{{{*/
     658
     659        for(Object* & object : femmodel->elements->objects){
     660                Element* element=xDynamicCast<Element*>(object);
     661                if (ismethod == 1) element->SmbSemicTransient(); // Inwoo's version.
     662                else element->SmbSemic(); // original SmbSEMIC
    662663        }
    663664
  • issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h

    r27297 r27498  
    2525void SmbDebrisMLx(FemModel* femmodel);
    2626/* SEMIC: */
    27 void SmbSemicx(FemModel* femmodel);
     27void SmbSemicx(FemModel* femmodel, int ismethod);
    2828/*GEMB: */
    2929void       Gembx(FemModel* femmodel);
  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r27491 r27498  
    524524syn keyword cConstant SmbAccurefEnum
    525525syn keyword cConstant SmbAdThreshEnum
     526syn keyword cConstant SmbAlbedoSchemeEnum
     527syn keyword cConstant SmbAlbedoSnowMaxEnum
     528syn keyword cConstant SmbAlbedoSnowMinEnum
     529syn keyword cConstant SmbAlbedoIceEnum
     530syn keyword cConstant SmbAlbedoLandEnum
    526531syn keyword cConstant SmbARMATimestepEnum
    527532syn keyword cConstant SmbARMAarOrderEnum
     
    529534syn keyword cConstant SmbAveragingEnum
    530535syn keyword cConstant SmbDesfacEnum
     536syn keyword cConstant SmbDesfacElevEnum
    531537syn keyword cConstant SmbDpermilEnum
    532538syn keyword cConstant SmbDsnowIdxEnum
     
    579585syn keyword cConstant SmbRunoffrefEnum
    580586syn keyword cConstant SmbSealevEnum
     587syn keyword cConstant SmbSemicMethodEnum
     588syn keyword cConstant SmbSemicHcritEnum
     589syn keyword cConstant SmbSemicRcritEnum
    581590syn keyword cConstant SmbStepsPerStepEnum
    582591syn keyword cConstant SmbSwIdxEnum
     
    10101019syn keyword cConstant SmbAccumulatedRefreezeEnum
    10111020syn keyword cConstant SmbAccumulatedRunoffEnum
     1021syn keyword cConstant SmbAlbedoEnum
     1022syn keyword cConstant SmbAlbedoSnowEnum
    10121023syn keyword cConstant SmbAEnum
    10131024syn keyword cConstant SmbAdiffEnum
     
    10571068syn keyword cConstant SmbGspEnum
    10581069syn keyword cConstant SmbGspiniEnum
     1070syn keyword cConstant SmbHIceEnum
     1071syn keyword cConstant SmbHSnowEnum
    10591072syn keyword cConstant SmbHrefEnum
    10601073syn keyword cConstant SmbIsInitializedEnum
     
    10631076syn keyword cConstant SmbMassBalanceSubstepEnum
    10641077syn keyword cConstant SmbMassBalanceTransientEnum
     1078syn keyword cConstant SmbMaskEnum
    10651079syn keyword cConstant SmbMeanLHFEnum
    10661080syn keyword cConstant SmbMeanSHFEnum
     
    10961110syn keyword cConstant SmbTEnum
    10971111syn keyword cConstant SmbTaEnum
     1112syn keyword cConstant SmbTampEnum
    10981113syn keyword cConstant SmbTeValueEnum
    10991114syn keyword cConstant SmbTemperaturesAnomalyEnum
     
    17261741syn keyword cType Cfsurfacesquare
    17271742syn keyword cType Channel
    1728 syn keyword cType classes
    17291743syn keyword cType Constraint
    17301744syn keyword cType Constraints
     
    17331747syn keyword cType ControlInput
    17341748syn keyword cType Covertree
     1749syn keyword cType DataSetParam
    17351750syn keyword cType DatasetInput
    1736 syn keyword cType DataSetParam
    17371751syn keyword cType Definition
    17381752syn keyword cType DependentObject
     
    17471761syn keyword cType ElementInput
    17481762syn keyword cType ElementMatrix
     1763syn keyword cType ElementVector
    17491764syn keyword cType Elements
    1750 syn keyword cType ElementVector
    17511765syn keyword cType ExponentialVariogram
    17521766syn keyword cType ExternalResult
     
    17551769syn keyword cType Friction
    17561770syn keyword cType Gauss
    1757 syn keyword cType GaussianVariogram
    1758 syn keyword cType gaussobjects
    17591771syn keyword cType GaussPenta
    17601772syn keyword cType GaussSeg
    17611773syn keyword cType GaussTetra
    17621774syn keyword cType GaussTria
     1775syn keyword cType GaussianVariogram
    17631776syn keyword cType GenericExternalResult
    17641777syn keyword cType GenericOption
     
    17771790syn keyword cType IssmDirectApplicInterface
    17781791syn keyword cType IssmParallelDirectApplicInterface
    1779 syn keyword cType krigingobjects
    17801792syn keyword cType Load
    17811793syn keyword cType Loads
     
    17881800syn keyword cType Matice
    17891801syn keyword cType Matlitho
    1790 syn keyword cType matrixobjects
    17911802syn keyword cType MatrixParam
    17921803syn keyword cType Misfit
     
    18011812syn keyword cType Observations
    18021813syn keyword cType Option
     1814syn keyword cType OptionUtilities
    18031815syn keyword cType Options
    1804 syn keyword cType OptionUtilities
    18051816syn keyword cType Param
    18061817syn keyword cType Parameters
     
    18161827syn keyword cType Regionaloutput
    18171828syn keyword cType Results
     1829syn keyword cType RiftStruct
    18181830syn keyword cType Riftfront
    1819 syn keyword cType RiftStruct
    18201831syn keyword cType SealevelGeometry
    18211832syn keyword cType Seg
    18221833syn keyword cType SegInput
     1834syn keyword cType SegRef
    18231835syn keyword cType Segment
    1824 syn keyword cType SegRef
    18251836syn keyword cType SpcDynamic
    18261837syn keyword cType SpcStatic
     
    18411852syn keyword cType Vertex
    18421853syn keyword cType Vertices
     1854syn keyword cType classes
     1855syn keyword cType gaussobjects
     1856syn keyword cType krigingobjects
     1857syn keyword cType matrixobjects
    18431858syn keyword cType AdjointBalancethickness2Analysis
    18441859syn keyword cType AdjointBalancethicknessAnalysis
  • issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

    r27491 r27498  
    518518        SmbAccurefEnum,
    519519        SmbAdThreshEnum,
     520        SmbAlbedoSchemeEnum,
     521        SmbAlbedoSnowMaxEnum,
     522        SmbAlbedoSnowMinEnum,
     523        SmbAlbedoIceEnum,
     524        SmbAlbedoLandEnum,
    520525        SmbARMATimestepEnum,
    521526   SmbARMAarOrderEnum,
     
    523528        SmbAveragingEnum,
    524529        SmbDesfacEnum,
     530        SmbDesfacElevEnum,
    525531        SmbDpermilEnum,
    526532        SmbDsnowIdxEnum,
     
    573579        SmbRunoffrefEnum,
    574580        SmbSealevEnum,
     581        SmbSemicMethodEnum,
     582        SmbSemicHcritEnum,
     583        SmbSemicRcritEnum,
    575584        SmbStepsPerStepEnum,
    576585        SmbSwIdxEnum,
     
    10061015        SmbAccumulatedRefreezeEnum,
    10071016        SmbAccumulatedRunoffEnum,
     1017        SmbAlbedoEnum,
     1018        SmbAlbedoSnowEnum,
    10081019        SmbAEnum,
    10091020        SmbAdiffEnum,
     
    10531064        SmbGspEnum,
    10541065        SmbGspiniEnum,
     1066        SmbHIceEnum,
     1067        SmbHSnowEnum,
    10551068        SmbHrefEnum,
    10561069        SmbIsInitializedEnum,
     
    10591072   SmbMassBalanceSubstepEnum,
    10601073   SmbMassBalanceTransientEnum,
     1074        SmbMaskEnum,
    10611075        SmbMeanLHFEnum,
    10621076        SmbMeanSHFEnum,
     
    10931107        SmbTEnum,
    10941108        SmbTaEnum,
     1109        SmbTampEnum,
    10951110        SmbTeValueEnum,
    10961111        SmbTemperaturesAnomalyEnum,
  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r27491 r27498  
    526526                case SmbAccurefEnum : return "SmbAccuref";
    527527                case SmbAdThreshEnum : return "SmbAdThresh";
     528                case SmbAlbedoSchemeEnum : return "SmbAlbedoScheme";
     529                case SmbAlbedoSnowMaxEnum : return "SmbAlbedoSnowMax";
     530                case SmbAlbedoSnowMinEnum : return "SmbAlbedoSnowMin";
     531                case SmbAlbedoIceEnum : return "SmbAlbedoIce";
     532                case SmbAlbedoLandEnum : return "SmbAlbedoLand";
    528533                case SmbARMATimestepEnum : return "SmbARMATimestep";
    529534                case SmbARMAarOrderEnum : return "SmbARMAarOrder";
     
    531536                case SmbAveragingEnum : return "SmbAveraging";
    532537                case SmbDesfacEnum : return "SmbDesfac";
     538                case SmbDesfacElevEnum : return "SmbDesfacElev";
    533539                case SmbDpermilEnum : return "SmbDpermil";
    534540                case SmbDsnowIdxEnum : return "SmbDsnowIdx";
     
    581587                case SmbRunoffrefEnum : return "SmbRunoffref";
    582588                case SmbSealevEnum : return "SmbSealev";
     589                case SmbSemicMethodEnum : return "SmbSemicMethod";
     590                case SmbSemicHcritEnum : return "SmbSemicHcrit";
     591                case SmbSemicRcritEnum : return "SmbSemicRcrit";
    583592                case SmbStepsPerStepEnum : return "SmbStepsPerStep";
    584593                case SmbSwIdxEnum : return "SmbSwIdx";
     
    10121021                case SmbAccumulatedRefreezeEnum : return "SmbAccumulatedRefreeze";
    10131022                case SmbAccumulatedRunoffEnum : return "SmbAccumulatedRunoff";
     1023                case SmbAlbedoEnum : return "SmbAlbedo";
     1024                case SmbAlbedoSnowEnum : return "SmbAlbedoSnow";
    10141025                case SmbAEnum : return "SmbA";
    10151026                case SmbAdiffEnum : return "SmbAdiff";
     
    10591070                case SmbGspEnum : return "SmbGsp";
    10601071                case SmbGspiniEnum : return "SmbGspini";
     1072                case SmbHIceEnum : return "SmbHIce";
     1073                case SmbHSnowEnum : return "SmbHSnow";
    10611074                case SmbHrefEnum : return "SmbHref";
    10621075                case SmbIsInitializedEnum : return "SmbIsInitialized";
     
    10651078                case SmbMassBalanceSubstepEnum : return "SmbMassBalanceSubstep";
    10661079                case SmbMassBalanceTransientEnum : return "SmbMassBalanceTransient";
     1080                case SmbMaskEnum : return "SmbMask";
    10671081                case SmbMeanLHFEnum : return "SmbMeanLHF";
    10681082                case SmbMeanSHFEnum : return "SmbMeanSHF";
     
    10981112                case SmbTEnum : return "SmbT";
    10991113                case SmbTaEnum : return "SmbTa";
     1114                case SmbTampEnum : return "SmbTamp";
    11001115                case SmbTeValueEnum : return "SmbTeValue";
    11011116                case SmbTemperaturesAnomalyEnum : return "SmbTemperaturesAnomaly";
  • issm/trunk-jpl/src/c/shared/Enum/Enumjl.vim

    r27491 r27498  
    517517syn keyword juliaConstC SmbAccurefEnum
    518518syn keyword juliaConstC SmbAdThreshEnum
     519syn keyword juliaConstC SmbAlbedoSchemeEnum
     520syn keyword juliaConstC SmbAlbedoSnowMaxEnum
     521syn keyword juliaConstC SmbAlbedoSnowMinEnum
     522syn keyword juliaConstC SmbAlbedoIceEnum
     523syn keyword juliaConstC SmbAlbedoLandEnum
    519524syn keyword juliaConstC SmbARMATimestepEnum
    520525syn keyword juliaConstC SmbARMAarOrderEnum
     
    522527syn keyword juliaConstC SmbAveragingEnum
    523528syn keyword juliaConstC SmbDesfacEnum
     529syn keyword juliaConstC SmbDesfacElevEnum
    524530syn keyword juliaConstC SmbDpermilEnum
    525531syn keyword juliaConstC SmbDsnowIdxEnum
     
    572578syn keyword juliaConstC SmbRunoffrefEnum
    573579syn keyword juliaConstC SmbSealevEnum
     580syn keyword juliaConstC SmbSemicMethodEnum
     581syn keyword juliaConstC SmbSemicHcritEnum
     582syn keyword juliaConstC SmbSemicRcritEnum
    574583syn keyword juliaConstC SmbStepsPerStepEnum
    575584syn keyword juliaConstC SmbSwIdxEnum
     
    10031012syn keyword juliaConstC SmbAccumulatedRefreezeEnum
    10041013syn keyword juliaConstC SmbAccumulatedRunoffEnum
     1014syn keyword juliaConstC SmbAlbedoEnum
     1015syn keyword juliaConstC SmbAlbedoSnowEnum
    10051016syn keyword juliaConstC SmbAEnum
    10061017syn keyword juliaConstC SmbAdiffEnum
     
    10501061syn keyword juliaConstC SmbGspEnum
    10511062syn keyword juliaConstC SmbGspiniEnum
     1063syn keyword juliaConstC SmbHIceEnum
     1064syn keyword juliaConstC SmbHSnowEnum
    10521065syn keyword juliaConstC SmbHrefEnum
    10531066syn keyword juliaConstC SmbIsInitializedEnum
     
    10561069syn keyword juliaConstC SmbMassBalanceSubstepEnum
    10571070syn keyword juliaConstC SmbMassBalanceTransientEnum
     1071syn keyword juliaConstC SmbMaskEnum
    10581072syn keyword juliaConstC SmbMeanLHFEnum
    10591073syn keyword juliaConstC SmbMeanSHFEnum
     
    10891103syn keyword juliaConstC SmbTEnum
    10901104syn keyword juliaConstC SmbTaEnum
     1105syn keyword juliaConstC SmbTampEnum
    10911106syn keyword juliaConstC SmbTeValueEnum
    10921107syn keyword juliaConstC SmbTemperaturesAnomalyEnum
  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r27491 r27498  
    538538              else if (strcmp(name,"SmbAccuref")==0) return SmbAccurefEnum;
    539539              else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
     540              else if (strcmp(name,"SmbAlbedoScheme")==0) return SmbAlbedoSchemeEnum;
     541              else if (strcmp(name,"SmbAlbedoSnowMax")==0) return SmbAlbedoSnowMaxEnum;
     542              else if (strcmp(name,"SmbAlbedoSnowMin")==0) return SmbAlbedoSnowMinEnum;
     543              else if (strcmp(name,"SmbAlbedoIce")==0) return SmbAlbedoIceEnum;
     544              else if (strcmp(name,"SmbAlbedoLand")==0) return SmbAlbedoLandEnum;
    540545              else if (strcmp(name,"SmbARMATimestep")==0) return SmbARMATimestepEnum;
    541546              else if (strcmp(name,"SmbARMAarOrder")==0) return SmbARMAarOrderEnum;
     
    543548              else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum;
    544549              else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum;
     550              else if (strcmp(name,"SmbDesfacElev")==0) return SmbDesfacElevEnum;
    545551              else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum;
    546552              else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum;
     
    593599              else if (strcmp(name,"SmbRunoffref")==0) return SmbRunoffrefEnum;
    594600              else if (strcmp(name,"SmbSealev")==0) return SmbSealevEnum;
     601              else if (strcmp(name,"SmbSemicMethod")==0) return SmbSemicMethodEnum;
     602              else if (strcmp(name,"SmbSemicHcrit")==0) return SmbSemicHcritEnum;
     603              else if (strcmp(name,"SmbSemicRcrit")==0) return SmbSemicRcritEnum;
    595604              else if (strcmp(name,"SmbStepsPerStep")==0) return SmbStepsPerStepEnum;
    596605              else if (strcmp(name,"SmbSwIdx")==0) return SmbSwIdxEnum;
     
    620629              else if (strcmp(name,"StressbalanceRequestedOutputs")==0) return StressbalanceRequestedOutputsEnum;
    621630              else if (strcmp(name,"StressbalanceRestol")==0) return StressbalanceRestolEnum;
    622               else if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"StressbalanceRiftPenaltyThreshold")==0) return StressbalanceRiftPenaltyThresholdEnum;
    623635              else if (strcmp(name,"StressbalanceShelfDampening")==0) return StressbalanceShelfDampeningEnum;
    624636              else if (strcmp(name,"ThermalForcingMonthlyEffects")==0) return ThermalForcingMonthlyEffectsEnum;
     
    629641              else if (strcmp(name,"ThermalNumRequestedOutputs")==0) return ThermalNumRequestedOutputsEnum;
    630642              else if (strcmp(name,"ThermalPenaltyFactor")==0) return ThermalPenaltyFactorEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
     643              else if (strcmp(name,"ThermalPenaltyLock")==0) return ThermalPenaltyLockEnum;
    635644              else if (strcmp(name,"ThermalPenaltyThreshold")==0) return ThermalPenaltyThresholdEnum;
    636645              else if (strcmp(name,"ThermalReltol")==0) return ThermalReltolEnum;
     
    743752              else if (strcmp(name,"Bed")==0) return BedEnum;
    744753              else if (strcmp(name,"BedGRD")==0) return BedGRDEnum;
    745               else if (strcmp(name,"BedEast")==0) return BedEastEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"BedEast")==0) return BedEastEnum;
    746758              else if (strcmp(name,"BedEastGRD")==0) return BedEastGRDEnum;
    747759              else if (strcmp(name,"BedNorth")==0) return BedNorthEnum;
     
    752764              else if (strcmp(name,"BottomPressureOld")==0) return BottomPressureOldEnum;
    753765              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
     766              else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    758767              else if (strcmp(name,"CalvingAblationrate")==0) return CalvingAblationrateEnum;
    759768              else if (strcmp(name,"CalvingMeltingrate")==0) return CalvingMeltingrateEnum;
     
    866875              else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum;
    867876              else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum;
    868               else if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"HydrologydcSedimentTransmitivity")==0) return HydrologydcSedimentTransmitivityEnum;
    869881              else if (strcmp(name,"HydrologyDrainageRate")==0) return HydrologyDrainageRateEnum;
    870882              else if (strcmp(name,"HydrologyEnglacialInput")==0) return HydrologyEnglacialInputEnum;
     
    875887              else if (strcmp(name,"HydrologyGapHeightYY")==0) return HydrologyGapHeightYYEnum;
    876888              else if (strcmp(name,"HydrologyHead")==0) return HydrologyHeadEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
     889              else if (strcmp(name,"HydrologyHeadOld")==0) return HydrologyHeadOldEnum;
    881890              else if (strcmp(name,"HydrologyMoulinInput")==0) return HydrologyMoulinInputEnum;
    882891              else if (strcmp(name,"HydrologyNeumannflux")==0) return HydrologyNeumannfluxEnum;
     
    989998              else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum;
    990999              else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum;
    991               else if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"BslcRate")==0) return BslcRateEnum;
    9921004              else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum;
    9931005              else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum;
     
    9981010              else if (strcmp(name,"SealevelUNorthEsa")==0) return SealevelUNorthEsaEnum;
    9991011              else if (strcmp(name,"SealevelchangeIndices")==0) return SealevelchangeIndicesEnum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum;
     1012              else if (strcmp(name,"SealevelchangeConvolutionVertices")==0) return SealevelchangeConvolutionVerticesEnum;
    10041013              else if (strcmp(name,"SealevelchangeAlphaIndex")==0) return SealevelchangeAlphaIndexEnum;
    10051014              else if (strcmp(name,"SealevelchangeAzimuthIndex")==0) return SealevelchangeAzimuthIndexEnum;
     
    10361045              else if (strcmp(name,"SmbAccumulatedRefreeze")==0) return SmbAccumulatedRefreezeEnum;
    10371046              else if (strcmp(name,"SmbAccumulatedRunoff")==0) return SmbAccumulatedRunoffEnum;
     1047              else if (strcmp(name,"SmbAlbedo")==0) return SmbAlbedoEnum;
     1048              else if (strcmp(name,"SmbAlbedoSnow")==0) return SmbAlbedoSnowEnum;
    10381049              else if (strcmp(name,"SmbA")==0) return SmbAEnum;
    10391050              else if (strcmp(name,"SmbAdiff")==0) return SmbAdiffEnum;
     
    10831094              else if (strcmp(name,"SmbGsp")==0) return SmbGspEnum;
    10841095              else if (strcmp(name,"SmbGspini")==0) return SmbGspiniEnum;
     1096              else if (strcmp(name,"SmbHIce")==0) return SmbHIceEnum;
     1097              else if (strcmp(name,"SmbHSnow")==0) return SmbHSnowEnum;
    10851098              else if (strcmp(name,"SmbHref")==0) return SmbHrefEnum;
    10861099              else if (strcmp(name,"SmbIsInitialized")==0) return SmbIsInitializedEnum;
     
    10891102              else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
    10901103              else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
     1104              else if (strcmp(name,"SmbMask")==0) return SmbMaskEnum;
    10911105              else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
    10921106              else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
     
    11071121              else if (strcmp(name,"SmbPrecipitationsReconstructed")==0) return SmbPrecipitationsReconstructedEnum;
    11081122              else if (strcmp(name,"SmbRain")==0) return SmbRainEnum;
    1109               else if (strcmp(name,"SmbRe")==0) return SmbReEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"SmbRe")==0) return SmbReEnum;
    11101127              else if (strcmp(name,"SmbRefreeze")==0) return SmbRefreezeEnum;
    11111128              else if (strcmp(name,"SmbReini")==0) return SmbReiniEnum;
     
    11211138              else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum;
    11221139              else if (strcmp(name,"SmbT")==0) return SmbTEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
     1140              else if (strcmp(name,"SmbTa")==0) return SmbTaEnum;
     1141              else if (strcmp(name,"SmbTamp")==0) return SmbTampEnum;
    11271142              else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum;
    11281143              else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum;
     
    12291244              else if (strcmp(name,"WaterfractionDrainageIntegrated")==0) return WaterfractionDrainageIntegratedEnum;
    12301245              else if (strcmp(name,"Waterfraction")==0) return WaterfractionEnum;
    1231               else if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
     1246         else stage=11;
     1247   }
     1248   if(stage==11){
     1249              if (strcmp(name,"Waterheight")==0) return WaterheightEnum;
    12321250              else if (strcmp(name,"WaterPressureArmaPerturbation")==0) return WaterPressureArmaPerturbationEnum;
    12331251              else if (strcmp(name,"WaterPressureValuesAutoregression")==0) return WaterPressureValuesAutoregressionEnum;
     
    12441262              else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    12451263              else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    1246          else stage=11;
    1247    }
    1248    if(stage==11){
    1249               if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
     1264              else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    12501265              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    12511266              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
     
    13521367              else if (strcmp(name,"AmrNeopz")==0) return AmrNeopzEnum;
    13531368              else if (strcmp(name,"AndroidFrictionCoefficient")==0) return AndroidFrictionCoefficientEnum;
    1354               else if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
     1369         else stage=12;
     1370   }
     1371   if(stage==12){
     1372              if (strcmp(name,"Arrhenius")==0) return ArrheniusEnum;
    13551373              else if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum;
    13561374              else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
     
    13671385              else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum;
    13681386              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    1369          else stage=12;
    1370    }
    1371    if(stage==12){
    1372               if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
     1387              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    13731388              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    13741389              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
     
    14751490              else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
    14761491              else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum;
    1477               else if (strcmp(name,"Gset")==0) return GsetEnum;
     1492         else stage=13;
     1493   }
     1494   if(stage==13){
     1495              if (strcmp(name,"Gset")==0) return GsetEnum;
    14781496              else if (strcmp(name,"Gsl")==0) return GslEnum;
    14791497              else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
     
    14901508              else if (strcmp(name,"HydrologyShreveAnalysis")==0) return HydrologyShreveAnalysisEnum;
    14911509              else if (strcmp(name,"HydrologySolution")==0) return HydrologySolutionEnum;
    1492          else stage=13;
    1493    }
    1494    if(stage==13){
    1495               if (strcmp(name,"HydrologySubsteps")==0) return HydrologySubstepsEnum;
     1510              else if (strcmp(name,"HydrologySubsteps")==0) return HydrologySubstepsEnum;
    14961511              else if (strcmp(name,"HydrologySubTime")==0) return HydrologySubTimeEnum;
    14971512              else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
     
    15981613              else if (strcmp(name,"None")==0) return NoneEnum;
    15991614              else if (strcmp(name,"Numberedcostfunction")==0) return NumberedcostfunctionEnum;
    1600               else if (strcmp(name,"NyeCO2")==0) return NyeCO2Enum;
     1615         else stage=14;
     1616   }
     1617   if(stage==14){
     1618              if (strcmp(name,"NyeCO2")==0) return NyeCO2Enum;
    16011619              else if (strcmp(name,"NyeH2O")==0) return NyeH2OEnum;
    16021620              else if (strcmp(name,"Numericalflux")==0) return NumericalfluxEnum;
     
    16131631              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    16141632              else if (strcmp(name,"P1P1")==0) return P1P1Enum;
    1615          else stage=14;
    1616    }
    1617    if(stage==14){
    1618               if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
     1633              else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum;
    16191634              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    16201635              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
     
    17211736              else if (strcmp(name,"TransientParam")==0) return TransientParamEnum;
    17221737              else if (strcmp(name,"TransientSolution")==0) return TransientSolutionEnum;
    1723               else if (strcmp(name,"Tria")==0) return TriaEnum;
     1738         else stage=15;
     1739   }
     1740   if(stage==15){
     1741              if (strcmp(name,"Tria")==0) return TriaEnum;
    17241742              else if (strcmp(name,"TriaInput")==0) return TriaInputEnum;
    17251743              else if (strcmp(name,"UzawaPressureAnalysis")==0) return UzawaPressureAnalysisEnum;
     
    17361754              else if (strcmp(name,"XYZ")==0) return XYZEnum;
    17371755              else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
    1738          else stage=15;
    1739    }
    1740    if(stage==15){
    1741               if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
     1756              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
    17421757              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    17431758              else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
  • issm/trunk-jpl/src/m/classes/SMBsemic.m

    r27454 r27498  
    1515                dailyairhumidity        = NaN;
    1616                dailytemperature        = NaN;
     17                Tamp              = NaN;
     18                mask              = NaN;
     19                hice              = NaN;
     20                hsnow             = NaN;
    1721                desfac                          = 0;
     22                desfacElevation   = 0;
    1823                rlaps                           = 0;
    1924                rdl                                     = 0;
     
    2227                averaging                       = 0;
    2328                requested_outputs       = {};
     29
     30                hcrit             = 0;
     31                rcrit             = 0;
     32
     33                % albedo
     34                albedo            = 0;
     35                albedo_snow       = 0;
     36                albedo_scheme     = 0;
     37                alb_smax = NaN;
     38                alb_smin = NaN;
     39                albi = NaN;
     40                albl = NaN;
     41
     42                % method
     43                ismethod          = 0;
    2444        end
    2545        methods
     
    5474                                disp('      no SMBsemic.s0gcm specified: values set as zero');
    5575                        end
    56 
     76                        self.Tamp       = 3*ones(md.mesh.numberofvertices,1);
     77                        self.albedo     = 0.8*ones(md.mesh.numberofvertices,1);
     78                        self.albedo_snow= 0.5*ones(md.mesh.numberofvertices,1);
     79                        self.hice       = zeros(md.mesh.numberofvertices,1);
     80                        self.hsnow      = 5*ones(md.mesh.numberofvertices,1);
    5781                end % }}}
    5882                function self = setdefaultparameters(self) % {{{
    5983
    60                         self.desfac             = -log(2.0)/1000;
    61                         self.rlaps              = 7.4;
    62                         self.rdl                        = 0.29;
     84                        self.albedo_scheme   = 0;
     85                        self.alb_smax = 0.79;
     86                        self.alb_smin = 0.6;
     87                        self.albi = 0.41;
     88                        self.albl = 0.07;
     89
     90                        self.hcrit = 0.028;% from Krapp et al. (2017)
     91                        self.rcrit = 0.85; % from Krapp et al. (2017)
     92               
     93                        self.desfac                   = -log(2.0)/1000;
     94                        self.desfacElevation = 2000;
     95                        self.rlaps                    = 7.4;
     96                        self.rdl                              = 0.29;
     97
     98                        self.ismethod        = 0;
    6399                        self.requested_outputs={'default'};
    64 
    65100                end % }}}zo
    66101                function md = checkconsistency(self,md,solution,analyses) % {{{
     
    80115                                md = checkfield(md,'fieldname','smb.dailyairhumidity','timeseries',1,'NaN',1,'Inf',1);
    81116                                md = checkfield(md,'fieldname','smb.dailytemperature','timeseries',1,'NaN',1,'Inf',1);
     117
     118                                % TODO: transient model should be merged with SEMIC developed by Ruckamp et al. (2018)
     119
     120                                md = checkfield(md,'fieldname','smb.ismethod','numel',1,'values',[0,1]);
     121                                if self.ismethod == 1 % transient mode.
     122                                        md = checkfield(md,'fieldname','smb.desfacElevation','>=',0,'numel',1);
     123
     124                                        md = checkfield(md,'fieldname','smb.albedo_scheme','NaN',1,'Inf',1,'numel',1,'values',[0,1,2,3,4]);
     125                                        md = checkfield(md,'fieldname','smb.alb_smax','>=',0,'NaN',1,'Inf',1,'numel',1);
     126                                        md = checkfield(md,'fieldname','smb.mask','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1],'values',[0,1,2]);
     127
     128                                        % initial values
     129                                        md = checkfield(md,'fieldname','smb.albedo','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
     130                                        md = checkfield(md,'fieldname','smb.albedo_snow','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
     131                                        md = checkfield(md,'fieldname','smb.alb_smax','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
     132                                        md = checkfield(md,'fieldname','smb.alb_smin','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
     133                                        md = checkfield(md,'fieldname','smb.albi','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
     134                                        md = checkfield(md,'fieldname','smb.albl','>=',0,'<=',1,'NaN',1,'Inf',1,'numel',1);
     135                                        md = checkfield(md,'fieldname','smb.hice','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
     136                                        md = checkfield(md,'fieldname','smb.hsnow','NaN',1,'Inf',1,'size',[md.mesh.numberofvertices, 1]);
     137                                end
    82138                        end
    83139                        md = checkfield(md,'fieldname','smb.steps_per_step','>=',1,'numel',[1]);
     
    106162                        fielddisplay(self,'rdl','longwave downward radiation decrease (default is 0.29 [W/m^2/km]; Marty et al. 2002)');
    107163                        fielddisplay(self,'s0gcm','GCM reference elevation; (default is 0) [m]');
     164
     165                        fielddisplay(self,'ismethod','method for calculating SMB with SEMIC. Default version of SEMIC is really slow. 0: steady, 1: transient (default: 0)');
     166                        if self.ismethod == 1 % tranisent mode
     167                                fielddisplay(self,'desfacElevation','desertification elevation (default is 2000 m; Vizcaino et al. 2010)');
     168                                fielddisplay(self,'Tamp','amplitude of diurnal cycle [K]');
     169                                fielddisplay(self,'albedo','initial albedo [no unit]');
     170                                fielddisplay(self,'albedo_snow','initial albedo for snow [no unit]');
     171                                fielddisplay(self,'hice','initial thickness of ice [unit: m]');
     172                                fielddisplay(self,'hsnow','initial thickness of snow [unit: m]');
     173                                fielddisplay(self,'mask','masking for albedo. 0: ocean, 1: land, 2: ice (default: 2)');
     174                                fielddisplay(self,'hcrit','critical snow height for albedo [unit: m]');
     175                                fielddisplay(self,'rcrit','critical refreezing height for albedo [no unit]');
     176
     177                                disp(sprintf('\nSEMIC albedo parameters.'));
     178                                fielddisplay(self,'albedo_scheme','albedo scheme for SEMIC. 0: none, 1: slater, 2: isba, 3: denby, 4: alex (default is 0)');
     179                                fielddisplay(self,'abl_smax','maximum snow albedo (default: 0.79)');
     180                                fielddisplay(self,'abl_smin','minimum snow albedo (default: 0.6)');
     181                                fielddisplay(self,'albi','background albeod for bare ice (default: 0.41)');
     182                                fielddisplay(self,'albl','background albeod for bare land (default: 0.07)');
     183                        end
     184
    108185                        fielddisplay(self, 'steps_per_step', 'number of smb steps per time step');
    109186                        fielddisplay(self, 'averaging', 'averaging methods from short to long steps');
     
    119196                        WriteData(fid,prefix,'name','md.smb.model','data',12,'format','Integer');
    120197
     198                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','ismethod','format','Integer','values',[0, 1]);
    121199                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfac','format','Double');
     200                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','desfacElevation','format','Double');
    122201                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','s0gcm','format','DoubleMat','mattype',1);
    123202                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','rlaps','format','Double');
     
    132211                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailyairhumidity','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
    133212                        WriteData(fid,prefix,'object',self,'class','smb','fieldname','dailytemperature','format','DoubleMat','mattype',1,'timeserieslength',md.mesh.numberofvertices+1,'yts',yts);
     213                        % TODO: transient mode should be merged with SEMIC developed by Ruckamp et al. (2018).
     214                        if self.ismethod == 1% transient mode.
     215                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','Tamp','format','DoubleMat','mattype',1);
     216                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','mask','format','DoubleMat','mattype',1);
     217                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','hice','format','DoubleMat','mattype',1);
     218                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','hsnow','format','DoubleMat','mattype',1);
     219
     220                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','hcrit','format','Double');
     221                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','rcrit','format','Double');
     222
     223                                %albedo...
     224                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo','format','DoubleMat','mattype',1);
     225                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_snow','format','DoubleMat','mattype',1);
     226                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','albedo_scheme','format','Integer');
     227                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','albi','format','Double');
     228                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','albl','format','Double');
     229                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smin','format','Double');
     230                                WriteData(fid,prefix,'object',self,'class','smb','fieldname','alb_smax','format','Double');
     231                        end
     232
    134233                        WriteData(fid,prefix,'object',self,'fieldname','steps_per_step','format','Integer');
    135234                        WriteData(fid,prefix,'object',self,'fieldname','averaging','format','Integer');
     235
    136236                        %process requested outputs
    137237                        outputs = self.requested_outputs;
Note: See TracChangeset for help on using the changeset viewer.