Ignore:
Timestamp:
11/19/19 16:40:16 (5 years ago)
Author:
Mathieu Morlighem
Message:

CHG: almost done with Gemb

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r24359 r24360  
    20672067
    20682068        /*update input*/
    2069         this->SetBoolInput(this->inputs2,name,constant);
     2069        this->SetIntInput(this->inputs2,name,constant);
    20702070}
    20712071/*}}}*/
    … …  
    33013301                        *parray_size      = 1;
    33023302                        break;
     3303                case ArrayInput2Enum:{
     3304                        int M;
     3305                        this->inputs2->GetArray(output_enum,this->lid,NULL,&M);
     3306                        *pinterpolation   = P0ArrayEnum;
     3307                        *pnodesperelement = 1;
     3308                        *parray_size      = M;
     3309                        }
     3310                        break;
    33033311                default:
    33043312                        _error_("Input type \""<<EnumToStringx(this->inputs2->GetInputObjectEnum(output_enum))<<"\" not supported yet");
    … …  
    33213329void       Element::ResultToMatrix(IssmDouble* values,int ncols,int output_enum){/*{{{*/
    33223330
    3323         Input2* input=this->GetInput2(output_enum);
    3324         if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
    3325 
    3326         _error_("not implemented yet");
    3327         //input->ResultToMatrix(values,ncols,this->Sid());
     3331        IssmDouble* array = NULL;
     3332        int         m;
     3333        this->inputs2->GetArray(output_enum,this->lid,&array,&m);
     3334        for(int i=0;i<m;i++) values[this->Sid()*ncols + i] = array[i];
     3335        xDelete<IssmDouble>(array);
    33283336
    33293337} /*}}}*/
    … …  
    36173625
    36183626        /*Intermediary variables: {{{*/
    3619         IssmDouble isinitialized=0.0;
     3627        bool       isinitialized;
    36203628        IssmDouble zTop=0.0;
    36213629        IssmDouble dzTop=0.0;
    … …  
    37353743        /*}}}*/
    37363744        /*Retrieve inputs: {{{*/
    3737         _error_("fix....");
    37383745        Input2 *zTop_input          = this->GetInput2(SmbZTopEnum);         _assert_(zTop_input);
    37393746        Input2 *dzTop_input         = this->GetInput2(SmbDzTopEnum);        _assert_(dzTop_input);
    … …  
    37493756        Input2 *teValue_input       = this->GetInput2(SmbTeValueEnum);      _assert_(teValue_input);
    37503757        Input2 *aValue_input        = this->GetInput2(SmbAValueEnum);       _assert_(aValue_input);
    3751 
    3752         TransientInput2 *Ta_input   = this->inputs2->GetTransientInput(SmbTaEnum);    _assert_(Ta_input);
    3753         TransientInput2 *V_input    = this->inputs2->GetTransientInput(SmbVEnum);     _assert_(V_input);
    3754         TransientInput2 *Dlwr_input = this->inputs2->GetTransientInput(SmbDlwrfEnum); _assert_(Dlwr_input);
    3755         TransientInput2 *Dswr_input = this->inputs2->GetTransientInput(SmbDswrfEnum); _assert_(Dswr_input);
    3756         TransientInput2 *P_input    = this->inputs2->GetTransientInput(SmbPEnum);     _assert_(P_input);
    3757         TransientInput2 *eAir_input = this->inputs2->GetTransientInput(SmbEAirEnum);  _assert_(eAir_input);
    3758         TransientInput2 *pAir_input = this->inputs2->GetTransientInput(SmbPAirEnum);  _assert_(pAir_input);
     3758        Input2 *EC_input            = this->GetInput2(SmbECiniEnum);        _assert_(EC_input);
    37593759
    37603760        /*Retrieve input values:*/
    … …  
    37783778
    37793779        /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
    3780         if(isinitialized==0.0){
     3780        if(!isinitialized){
    37813781                if(VerboseSmb() && this->Sid()==0)_printf0_("smb core: Initializing grid\n");
    37823782                //if(this->Sid()==1) for(int i=0;i<m;i++)_printf_("z[" << i << "]=" <<
    37833783                //dz[i] << "\n");
    37843784
    3785                 DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDziniEnum)); _assert_(dz_input);
    3786                 DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDiniEnum));_assert_(d_input);
    3787                 DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReiniEnum));_assert_(re_input);
    3788                 DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdniniEnum));_assert_(gdn_input);
    3789                 DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspiniEnum));_assert_(gsp_input);
    3790                 DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECiniEnum));_assert_(EC_input);
    3791                 DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWiniEnum));_assert_(W_input);
    3792                 DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAiniEnum));_assert_(a_input);
    3793                 DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTiniEnum));_assert_(T_input);
    3794 
    3795                 dz_input->GetValues(&dzini,&m);
    3796                 d_input->GetValues(&dini,&m);
    3797                 re_input->GetValues(&reini,&m);
    3798                 gdn_input->GetValues(&gdnini,&m);
    3799                 gsp_input->GetValues(&gspini,&m);
    3800                 EC_input->GetInputValue(&EC);
    3801                 W_input->GetValues(&Wini,&m);
    3802                 a_input->GetValues(&aini,&m);
    3803                 T_input->GetValues(&Tini,&m);
     3785                this->inputs2->GetArray(SmbDziniEnum,this->lid,&dzini,&m);
     3786                this->inputs2->GetArray(SmbDiniEnum,this->lid,&dini,&m);
     3787                this->inputs2->GetArray(SmbReiniEnum,this->lid,&reini,&m);
     3788                this->inputs2->GetArray(SmbGdniniEnum,this->lid,&gdnini,&m);
     3789                this->inputs2->GetArray(SmbGspiniEnum,this->lid,&gspini,&m);
     3790                this->inputs2->GetArray(SmbWiniEnum,this->lid,&Wini,&m);
     3791                this->inputs2->GetArray(SmbAiniEnum,this->lid,&aini,&m);
     3792                this->inputs2->GetArray(SmbTiniEnum,this->lid,&Tini,&m);
     3793                EC_input->GetInputAverage(&EC);
    38043794
    38053795                /*Retrive the correct value of m (without the zeroes at the end)*/
    3806                 Input* Size_input=this->GetInput(SmbSizeiniEnum); _assert_(Size_input);
    3807                 Size_input->GetInputValue(&m);
     3796                this->GetInput2Value(&m,SmbSizeiniEnum);
    38083797
    38093798                if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
    … …  
    38133802                        GembgridInitialize(&dz, &m, zTop, dzTop, zMax, zY);
    38143803
    3815                         d = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3]
    3816                         re = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0];         //set grain size to old snow [mm]
    3817                         gdn = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0];         //set grain dentricity to old snow
    3818                         gsp = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0];         //set grain sphericity to old snow
    3819                         W = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
    3820                         a = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
    3821                         T = xNewZeroInit<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
     3804                        d = xNew<IssmDouble>(m); for(int i=0;i<m;i++)d[i]=dini[0]; //ice density [kg m-3]
     3805                        re = xNew<IssmDouble>(m); for(int i=0;i<m;i++)re[i]=reini[0];         //set grain size to old snow [mm]
     3806                        gdn = xNew<IssmDouble>(m); for(int i=0;i<m;i++)gdn[i]=gdnini[0];         //set grain dentricity to old snow
     3807                        gsp = xNew<IssmDouble>(m); for(int i=0;i<m;i++)gsp[i]=gspini[0];         //set grain sphericity to old snow
     3808                        W = xNew<IssmDouble>(m); for(int i=0;i<m;i++)W[i]=Wini[0];             //set water content to zero [kg m-2]
     3809                        a = xNew<IssmDouble>(m); for(int i=0;i<m;i++)a[i]=aini[0];         //set albedo equal to fresh snow [fraction]
     3810                        T = xNew<IssmDouble>(m); for(int i=0;i<m;i++)T[i]=Tmean;         //set initial grid cell temperature to the annual mean temperature [K]
    38223811                        /*/!\ Default value of T can not be retrived from SMBgemb.m (like other snow properties)
    38233812                         *    because don't know Tmean yet when set default values.
    … …  
    38303819                        //            if(VerboseSmb() && this->Sid()==0)_printf0_("Snow properties initialized w RESTART values\n");
    38313820
    3832                         dz = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i];
    3833                         d = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i];
    3834                         re = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i];
    3835                         gdn = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i];
    3836                         gsp = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i];
    3837                         W = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
    3838                         a = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
    3839                         T = xNewZeroInit<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];
     3821                        dz = xNew<IssmDouble>(m);for(int i=0;i<m;i++)dz[i]=dzini[i];
     3822                        d = xNew<IssmDouble>(m);for(int i=0;i<m;i++)d[i]=dini[i];
     3823                        re = xNew<IssmDouble>(m);for(int i=0;i<m;i++)re[i]=reini[i];
     3824                        gdn = xNew<IssmDouble>(m);for(int i=0;i<m;i++)gdn[i]=gdnini[i];
     3825                        gsp = xNew<IssmDouble>(m);for(int i=0;i<m;i++)gsp[i]=gspini[i];
     3826                        W = xNew<IssmDouble>(m);for(int i=0;i<m;i++)W[i]=Wini[i];
     3827                        a = xNew<IssmDouble>(m);for(int i=0;i<m;i++)a[i]=aini[i];
     3828                        T = xNew<IssmDouble>(m);for(int i=0;i<m;i++)T[i]=Tini[i];
    38403829
    38413830                        //fixed lower temperature bounday condition - T is fixed
     3831                        _assert_(m>0);
    38423832                        T_bottom=T[m-1];
    38433833                }
    38443834
    38453835                /*Flag the initialization:*/
    3846                 this->AddInput(new DoubleInput(SmbIsInitializedEnum,1.0));
     3836                this->SetBoolInput(this->inputs2,SmbIsInitializedEnum,true);
    38473837        }
    38483838        else{
    38493839                /*Recover inputs: */
    3850                 DoubleArrayInput* dz_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDzEnum)); _assert_(dz_input);
    3851                 DoubleArrayInput* d_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbDEnum));_assert_(d_input);
    3852                 DoubleArrayInput* re_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbReEnum));_assert_(re_input);
    3853                 DoubleArrayInput* gdn_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGdnEnum));_assert_(gdn_input);
    3854                 DoubleArrayInput* gsp_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbGspEnum));_assert_(gsp_input);
    3855                 DoubleInput* EC_input= dynamic_cast<DoubleInput*>(this->GetInput(SmbECEnum));_assert_(EC_input);
    3856                 DoubleArrayInput* W_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbWEnum));_assert_(W_input);
    3857                 DoubleArrayInput* a_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbAEnum));_assert_(a_input);
    3858                 DoubleArrayInput* T_input= dynamic_cast<DoubleArrayInput*>(this->GetInput(SmbTEnum));_assert_(T_input);
    3859 
    3860                 /*Recover arrays: */
    3861                 dz_input->GetValues(&dz,&m);
    3862                 d_input->GetValues(&d,&m);
    3863                 re_input->GetValues(&re,&m);
    3864                 gdn_input->GetValues(&gdn,&m);
    3865                 gsp_input->GetValues(&gsp,&m);
    3866                 EC_input->GetInputValue(&EC);
    3867                 W_input->GetValues(&W,&m);
    3868                 a_input->GetValues(&a,&m);
    3869                 T_input->GetValues(&T,&m);
     3840                this->inputs2->GetArray(SmbDzEnum,this->lid,&dzini,&m);
     3841                this->inputs2->GetArray(SmbDEnum,this->lid,&dini,&m);
     3842                this->inputs2->GetArray(SmbReEnum,this->lid,&reini,&m);
     3843                this->inputs2->GetArray(SmbGdnEnum,this->lid,&gdnini,&m);
     3844                this->inputs2->GetArray(SmbGspEnum,this->lid,&gspini,&m);
     3845                this->inputs2->GetArray(SmbWEnum,this->lid,&Wini,&m);
     3846                this->inputs2->GetArray(SmbAEnum,this->lid,&aini,&m);
     3847                this->inputs2->GetArray(SmbTEnum,this->lid,&Tini,&m);
     3848                EC_input->GetInputAverage(&EC);
    38703849
    38713850                //fixed lower temperature bounday condition - T is fixed
     3851                _assert_(m>0);
    38723852                T_bottom=T[m-1];
    38733853
    … …  
    38883868        if (isclimatology){
    38893869                //If this is a climatology, we need to repeat the forcing after the final time
    3890                 offsetend = Ta_input->GetTimeInputOffset(finaltime);
    3891                 time0     = Ta_input->GetTimeByOffset(-1);
    3892                 timeend   = Ta_input->GetTimeByOffset(offsetend);
     3870                TransientInput2* Ta_input_tr  = this->inputs2->GetTransientInput(SmbTaEnum);    _assert_(Ta_input_tr);
     3871                offsetend = Ta_input_tr->GetTimeInputOffset(finaltime);
     3872                time0     = Ta_input_tr->GetTimeByOffset(-1);
     3873                timeend   = Ta_input_tr->GetTimeByOffset(offsetend);
    38933874                if (time>time0 & timeend>time0){
    38943875                        delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
    … …  
    39033884                if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");
    39043885
    3905                 Input2* Ta_input2  = this->GetInput2(SmbTaEnum,t-time+timeclim);    _assert_(Ta_input2);
    3906                 Input2 *V_input2   = this->GetInput2(SmbVEnum,t-time+timeclim);     _assert_(V_input);
    3907                 Input2 *Dlwr_input2= this->GetInput2(SmbDlwrfEnum,t-time+timeclim); _assert_(Dlwr_input);
    3908                 Input2 *Dswr_input2= this->GetInput2(SmbDswrfEnum,t-time+timeclim); _assert_(Dswr_input);
    3909                 Input2 *P_input2   = this->GetInput2(SmbPEnum,t-time+timeclim);     _assert_(P_input);
    3910                 Input2 *eAir_input2= this->GetInput2(SmbEAirEnum,t-time+timeclim);  _assert_(eAir_input);
    3911                 Input2 *pAir_input2= this->GetInput2(SmbPAirEnum,t-time+timeclim);  _assert_(pAir_input);
     3886                Input2* Ta_input  = this->GetInput2(SmbTaEnum,t-time+timeclim);    _assert_(Ta_input);
     3887                Input2 *V_input   = this->GetInput2(SmbVEnum,t-time+timeclim);     _assert_(V_input);
     3888                Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,t-time+timeclim); _assert_(Dlwr_input);
     3889                Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,t-time+timeclim); _assert_(Dswr_input);
     3890                Input2 *P_input   = this->GetInput2(SmbPEnum,t-time+timeclim);     _assert_(P_input);
     3891                Input2 *eAir_input= this->GetInput2(SmbEAirEnum,t-time+timeclim);  _assert_(eAir_input);
     3892                Input2 *pAir_input= this->GetInput2(SmbPAirEnum,t-time+timeclim);  _assert_(pAir_input);
    39123893
    39133894                /*extract daily data:{{{*/
    3914                 Ta_input2->GetInputValue(&Ta,gauss);//screen level air temperature [K]
     3895                Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
    39153896                V_input->GetInputValue(&V,gauss);  //wind speed [m s-1]
    39163897                Dlwr_input->GetInputValue(&dlw,gauss);   //downward longwave radiation flux [W m-2]
Note: See TracChangeset for help on using the changeset viewer.