Changeset 24683


Ignore:
Timestamp:
04/01/20 18:53:08 (5 years ago)
Author:
schlegel
Message:

CHG: speed up GEMB, move time loop into Gembx outside of element loop

Location:
issm/trunk-jpl/src/c
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r24652 r24683  
    562562        this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum);
    563563
    564                 /*Recover present day temperature and precipitation*/
     564        /*Recover present day temperature and precipitation*/
    565565        DatasetInput2 *dinput3 = NULL;
    566566        DatasetInput2 *dinput4 = NULL;
    … …  
    35883588}
    35893589/*}}}*/
    3590 void       Element::SmbGemb(){/*{{{*/
     3590void       Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/
    35913591
    35923592        /*only compute SMB at the surface: */
    … …  
    36053605        IssmDouble C=0.0;
    36063606        IssmDouble Tz,Vz=0.0;
    3607         IssmDouble time,dt,starttime,finaltime;
    3608         IssmDouble timeclim=0.0;
    3609         IssmDouble t,smb_dt;
    36103607        IssmDouble yts;
    36113608        IssmDouble Ta=0.0;
    … …  
    36183615        IssmDouble teValue=1.0;
    36193616        IssmDouble aValue=0.0;
     3617        IssmDouble dt,time,smb_dt;
    36203618        int        aIdx=0;
    36213619        int        denIdx=0;
    … …  
    36273625        IssmDouble dayEC=0.0;
    36283626        IssmDouble initMass=0.0;
    3629         IssmDouble sumR=0.0;
    3630         IssmDouble sumM=0.0;
    3631         IssmDouble sumEC=0.0;
    3632         IssmDouble sumP=0.0;
    3633         IssmDouble sumW=0.0;
    3634         IssmDouble sumMassAdd=0.0;
    3635         IssmDouble sumdz_add=0.0;
     3627   IssmDouble sumR=0.0;
     3628   IssmDouble sumM=0.0;
     3629   IssmDouble sumEC=0.0;
     3630   IssmDouble sumP=0.0;
     3631   IssmDouble sumW=0.0;
     3632   IssmDouble sumMassAdd=0.0;
    36363633        IssmDouble fac=0.0;
    36373634        IssmDouble sumMass=0.0;
    … …  
    36423639        IssmDouble thermo_scaling=1.0;
    36433640        IssmDouble adThresh=1023.0;
    3644         int offsetend=-1;
    3645         IssmDouble time0, timeend, delta;
    36463641        /*}}}*/
    36473642        /*Output variables:{{{ */
    … …  
    36763671        IssmDouble* Tini = NULL;
    36773672        int         m=0;
    3678         int         count=0;
    36793673        /*}}}*/
    36803674
    … …  
    36873681        parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
    36883682        parameters->FindParam(&yts,ConstantsYtsEnum);
    3689         parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
    3690         parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    36913683        parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
    36923684        parameters->FindParam(&aIdx,SmbAIdxEnum);
    … …  
    36983690        parameters->FindParam(&t0dry,SmbT0dryEnum);
    36993691        parameters->FindParam(&K,SmbKEnum);
    3700         parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
    37013692        parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum);
    37023693        parameters->FindParam(&isalbedo,SmbIsalbedoEnum);
    … …  
    37233714        Input2 *Tz_input            = this->GetInput2(SmbTzEnum);           _assert_(Tz_input);
    37243715        Input2 *Vz_input            = this->GetInput2(SmbVzEnum);           _assert_(Vz_input);
    3725         Input2 *teValue_input       = this->GetInput2(SmbTeValueEnum);      _assert_(teValue_input);
    3726         Input2 *aValue_input        = this->GetInput2(SmbAValueEnum);       _assert_(aValue_input);
    37273716        Input2 *EC_input            = NULL;
    37283717
    … …  
    37423731        Tz_input->GetInputValue(&Tz,gauss);
    37433732        Vz_input->GetInputValue(&Vz,gauss);
    3744         teValue_input->GetInputValue(&teValue,gauss);
    3745         aValue_input->GetInputValue(&aValue,gauss);
    37463733        /*}}}*/
    37473734
    … …  
    37633750                EC_input->GetInputAverage(&EC);
    37643751
    3765                 /*Retrive the correct value of m (without the zeroes at the end)*/
     3752                /*Retrieve the correct value of m (without the zeroes at the end)*/
    37663753                this->GetInput2Value(&m,SmbSizeiniEnum);
    37673754
    … …  
    38163803                this->inputs2->GetArray(SmbAEnum,this->lid,&a,&m);
    38173804                this->inputs2->GetArray(SmbTEnum,this->lid,&T,&m);
    3818                 EC_input = this->GetInput2(SmbECEnum);  _assert_(EC_input);
     3805                EC_input = this->GetInput2(SmbECDtEnum);  _assert_(EC_input);
    38193806                EC_input->GetInputAverage(&EC);
    3820                 EC=EC*dt*rho_ice;
    38213807
    38223808                //fixed lower temperature bounday condition - T is fixed
    … …  
    38303816        // initialize cumulative variables
    38313817        sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
    3832         sumdz_add=0;
    38333818
    38343819        //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
    3835         //go back to time - deltaT:
    3836         time-=dt;
    3837 
    3838         timeclim=time;
    3839         if (isclimatology){
    3840                 //If this is a climatology, we need to repeat the forcing after the final time
    3841                 TransientInput2* Ta_input_tr  = this->inputs2->GetTransientInput(SmbTaEnum);    _assert_(Ta_input_tr);
    3842                 offsetend = Ta_input_tr->GetTimeInputOffset(finaltime);
    3843                 time0     = Ta_input_tr->GetTimeByOffset(-1);
    3844                 timeend   = Ta_input_tr->GetTimeByOffset(offsetend);
    3845                 if (time>time0 & timeend>time0){
    3846                         delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
    3847                         timeclim=time0+delta;
    3848                 }
    3849         }
    3850 
    3851         /*Start loop: */
    3852         count=1;
    3853         for (t=time;t<time+dt;t=t+smb_dt){
    3854 
    3855                 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");
    3856 
    3857                 Input2 *Ta_input  = this->GetInput2(SmbTaEnum,t-time+timeclim);    _assert_(Ta_input);
    3858                 Input2 *V_input   = this->GetInput2(SmbVEnum,t-time+timeclim);     _assert_(V_input);
    3859                 Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,t-time+timeclim); _assert_(Dlwr_input);
    3860                 Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,t-time+timeclim); _assert_(Dswr_input);
    3861                 Input2 *P_input   = this->GetInput2(SmbPEnum,t-time+timeclim);     _assert_(P_input);
    3862                 Input2 *eAir_input= this->GetInput2(SmbEAirEnum,t-time+timeclim);  _assert_(eAir_input);
    3863                 Input2 *pAir_input= this->GetInput2(SmbPAirEnum,t-time+timeclim);  _assert_(pAir_input);
    3864 
    3865                 /*extract daily data:{{{*/
    3866                 Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
    3867                 V_input->GetInputValue(&V,gauss);  //wind speed [m s-1]
    3868                 Dlwr_input->GetInputValue(&dlw,gauss);   //downward longwave radiation flux [W m-2]
    3869                 Dswr_input->GetInputValue(&dsw,gauss);   //downward shortwave radiation flux [W m-2]
    3870                 P_input->GetInputValue(&P,gauss);        //precipitation [kg m-2]
    3871                 eAir_input->GetInputValue(&eAir,gauss);  //screen level vapor pressure [Pa]
    3872                 pAir_input->GetInputValue(&pAir,gauss);  // screen level air pressure [Pa]
    3873                 teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
    3874                 aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
    3875                 //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
    3876                 /*}}}*/
    3877 
    3878                 /*Snow grain metamorphism:*/
    3879                 if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
    3880 
    3881                 /*Snow, firn and ice albedo:*/
    3882                 if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
    3883 
    3884                 /*Distribution of absorbed short wave radation with depth:*/
    3885                 if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
    3886 
    3887                 /*Calculate net shortwave [W m-2]*/
    3888                 netSW = netSW + cellsum(swf,m)*smb_dt/dt;
    3889 
    3890                 /*Thermal profile computation:*/
    3891                 if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
    3892 
    3893                 /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
    3894                  * need to fix this in case all or more of cell evaporates */
    3895                 dz[0] = dz[0] + EC / d[0];
    3896 
    3897                 /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
    3898                 if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
    3899 
    3900                 /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
    3901                  * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
    3902                 if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
    3903 
    3904                 /*Allow non-melt densification and determine compaction [m]*/
    3905                 if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
    3906 
    3907                 /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
    3908                  * sub-time step in thermo equations*/
    3909                 //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
    3910 
    3911                 /*Calculate net longwave [W m-2]*/
    3912                 meanULW = meanULW + ulw*smb_dt/dt;
    3913                 netLW = netLW + (dlw - ulw)*smb_dt/dt;
    3914 
    3915                 /*Calculate turbulent heat fluxes [W m-2]*/
    3916                 if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
    3917 
    3918                 /*Verbose some results in debug mode: {{{*/
    3919                 if(VerboseSmb() && 0){
    3920                         _printf_("smb log: count[" << count << "] m[" << m << "] "
    3921                                                 << setprecision(16)   << "T[" << cellsum(T,m)  << "] "
    3922                                                 << "d[" << cellsum(d,m)  << "] "
    3923                                                 << "dz[" << cellsum(dz,m)  << "] "
    3924                                                 << "a[" << cellsum(a,m)  << "] "
    3925                                                 << "W[" << cellsum(W,m)  << "] "
    3926                                                 << "re[" << cellsum(re,m)  << "] "
    3927                                                 << "gdn[" << cellsum(gdn,m)  << "] "
    3928                                                 << "gsp[" << cellsum(gsp,m)  << "] "
    3929                                                 << "swf[" << netSW << "] "
    3930                                                 << "lwf[" << netLW << "] "
    3931                                                 << "a[" << a << "] "
    3932                                                 << "te[" << teValue << "] "
    3933                                                 << "\n");
    3934                 } /*}}}*/
    3935 
    3936                 meanLHF = meanLHF + lhf*smb_dt/dt;
    3937                 meanSHF = meanSHF + shf*smb_dt/dt;
    3938 
    3939                 /*Sum component mass changes [kg m-2]*/
    3940                 sumMassAdd = mAdd + sumMassAdd;
    3941                 sumM = M + sumM;
    3942                 sumR = R + sumR;
    3943                 sumW = cellsum(W,m);
    3944                 sumP = P +  sumP;
    3945                 sumEC = sumEC + EC;  // evap (-)/cond(+)
    3946 
    3947                 /*Calculate total system mass:*/
    3948                 sumMass=0;
    3949                 fac=0;
    3950                 for(int i=0;i<m;i++){
    3951                         sumMass += dz[i]*d[i];
    3952                         fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
    3953                 }
    3954 
    3955                 #if defined(_HAVE_AD_)
    3956                 /*we want to avoid the round operation at all cost. Not differentiable.*/
    3957                 _error_("not implemented yet");
    3958                 #else
    3959                 dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
    3960                 dMass = round(dMass * 100.0)/100.0;
    3961 
    3962                 /*Check mass conservation:*/
    3963                 if (dMass != 0.0) _printf_("total system mass not conserved in MB function \n");
    3964                 #endif
    3965 
    3966                 /*Check bottom grid cell T is unchanged:*/
    3967                 if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
    3968                         if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
    3969                 }
    3970 
    3971                 /*Free ressources: */
    3972                 xDelete<IssmDouble>(swf);
    3973 
    3974                 /*increase counter:*/
    3975                 count++;
    3976         } //for (t=time;t<time+dt;t=t+smb_dt)
     3820   //go back to time - deltaT:
     3821   time-=dt;
     3822
     3823        if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << timeinputs/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");
     3824
     3825        /*Get daily accumulated inputs {{{*/
     3826        if (count>1){
     3827                Input2 *sumEC_input         = this->GetInput2(SmbECEnum);  _assert_(sumEC_input);
     3828                Input2 *sumM_input          = this->GetInput2(SmbMeltEnum);  _assert_(sumM_input);
     3829                Input2 *sumR_input          = this->GetInput2(SmbRunoffEnum);  _assert_(sumR_input);
     3830                Input2 *sumP_input          = this->GetInput2(SmbPrecipitationEnum);  _assert_(sumP_input);
     3831                Input2 *ULW_input           = this->GetInput2(SmbMeanULWEnum);  _assert_(ULW_input);
     3832                Input2 *LW_input            = this->GetInput2(SmbNetLWEnum);  _assert_(LW_input);
     3833                Input2 *SW_input            = this->GetInput2(SmbNetSWEnum);  _assert_(SW_input);
     3834                Input2 *LHF_input           = this->GetInput2(SmbMeanLHFEnum);  _assert_(LHF_input);
     3835                Input2 *SHF_input           = this->GetInput2(SmbMeanSHFEnum);  _assert_(SHF_input);
     3836                Input2 *DzAdd_input         = this->GetInput2(SmbDzAddEnum);  _assert_(DzAdd_input);
     3837                Input2 *MassAdd_input       = this->GetInput2(SmbMAddEnum);  _assert_(MassAdd_input);
     3838                Input2 *InitMass_input      = this->GetInput2(SmbMInitnum);  _assert_(InitMass_input);
     3839
     3840                ULW_input->GetInputAverage(&meanULW);
     3841                LW_input->GetInputAverage(&netLW);
     3842                SW_input->GetInputAverage(&netSW);
     3843                LHF_input->GetInputAverage(&meanLHF);
     3844                SHF_input->GetInputAverage(&meanSHF);
     3845                DzAdd_input->GetInputAverage(&dz_add);
     3846                MassAdd_input->GetInputAverage(&sumMassAdd);
     3847                sumMassAdd=sumMassAdd*dt;
     3848                InitMass_input->GetInputAverage(&initMass);
     3849                sumEC_input->GetInputAverage(&sumEC);
     3850                sumEC=sumEC*dt*rho_ice;
     3851                sumM_input->GetInputAverage(&sumM);
     3852                sumM=sumM*dt*rho_ice;
     3853                sumR_input->GetInputAverage(&sumR);
     3854                sumR=sumR*dt*rho_ice;
     3855                sumP_input->GetInputAverage(&sumP);
     3856                sumP=sumP*dt*rho_ice;
     3857        }
     3858        /*}}}*/
     3859
     3860        // Get time forcing inputs
     3861        Input2 *Ta_input  = this->GetInput2(SmbTaEnum,timeinputs);    _assert_(Ta_input);
     3862        Input2 *V_input   = this->GetInput2(SmbVEnum,timeinputs);     _assert_(V_input);
     3863        Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input);
     3864        Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,timeinputs); _assert_(Dswr_input);
     3865        Input2 *P_input   = this->GetInput2(SmbPEnum,timeinputs);     _assert_(P_input);
     3866        Input2 *eAir_input= this->GetInput2(SmbEAirEnum,timeinputs);  _assert_(eAir_input);
     3867        Input2 *pAir_input= this->GetInput2(SmbPAirEnum,timeinputs);  _assert_(pAir_input);
     3868        Input2 *teValue_input= this->GetInput2(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
     3869        Input2 *aValue_input= this->GetInput2(SmbAValueEnum,timeinputs); _assert_(aValue_input);
     3870
     3871        /*extract daily data:{{{*/
     3872        Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
     3873        V_input->GetInputValue(&V,gauss);  //wind speed [m s-1]
     3874        Dlwr_input->GetInputValue(&dlw,gauss);   //downward longwave radiation flux [W m-2]
     3875        Dswr_input->GetInputValue(&dsw,gauss);   //downward shortwave radiation flux [W m-2]
     3876        P_input->GetInputValue(&P,gauss);        //precipitation [kg m-2]
     3877        eAir_input->GetInputValue(&eAir,gauss);  //screen level vapor pressure [Pa]
     3878        pAir_input->GetInputValue(&pAir,gauss);  // screen level air pressure [Pa]
     3879        teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
     3880        aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
     3881        //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
     3882        /*}}}*/
     3883
     3884        /*Snow grain metamorphism:*/
     3885        if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
     3886
     3887        /*Snow, firn and ice albedo:*/
     3888        if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
     3889
     3890        /*Distribution of absorbed short wave radation with depth:*/
     3891        if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
     3892
     3893        /*Calculate net shortwave [W m-2]*/
     3894        netSW = netSW + cellsum(swf,m)*smb_dt/dt;
     3895
     3896        /*Thermal profile computation:*/
     3897        if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
     3898
     3899        /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
     3900         * need to fix this in case all or more of cell evaporates */
     3901        dz[0] = dz[0] + EC / d[0];
     3902
     3903        /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
     3904        if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
     3905
     3906        /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
     3907         * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
     3908        if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
     3909
     3910        /*Allow non-melt densification and determine compaction [m]*/
     3911        if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
     3912
     3913        /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
     3914         * sub-time step in thermo equations*/
     3915        //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
     3916
     3917        /*Calculate net longwave [W m-2]*/
     3918        meanULW = meanULW + ulw*smb_dt/dt;
     3919        netLW = netLW + (dlw - ulw)*smb_dt/dt;
     3920
     3921        /*Calculate turbulent heat fluxes [W m-2]*/
     3922        if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
     3923
     3924        /*Verbose some results in debug mode: {{{*/
     3925        if(VerboseSmb() && 0){
     3926                _printf_("smb log: count[" << count << "] m[" << m << "] "
     3927                                        << setprecision(16)   << "T[" << cellsum(T,m)  << "] "
     3928                                        << "d[" << cellsum(d,m)  << "] "
     3929                                        << "dz[" << cellsum(dz,m)  << "] "
     3930                                        << "a[" << cellsum(a,m)  << "] "
     3931                                        << "W[" << cellsum(W,m)  << "] "
     3932                                        << "re[" << cellsum(re,m)  << "] "
     3933                                        << "gdn[" << cellsum(gdn,m)  << "] "
     3934                                        << "gsp[" << cellsum(gsp,m)  << "] "
     3935                                        << "swf[" << netSW << "] "
     3936                                        << "lwf[" << netLW << "] "
     3937                                        << "a[" << a << "] "
     3938                                        << "te[" << teValue << "] "
     3939                                        << "\n");
     3940        } /*}}}*/
     3941
     3942        meanLHF = meanLHF + lhf*smb_dt/dt;
     3943        meanSHF = meanSHF + shf*smb_dt/dt;
     3944
     3945        /*Sum component mass changes [kg m-2]*/
     3946        sumMassAdd = mAdd + sumMassAdd;
     3947        sumM = M + sumM;
     3948        sumR = R + sumR;
     3949        sumW = cellsum(W,m);
     3950        sumP = P +  sumP;
     3951        sumEC = sumEC + EC;  // evap (-)/cond(+)
     3952
     3953        /*Calculate total system mass:*/
     3954        sumMass=0;
     3955        fac=0;
     3956        for(int i=0;i<m;i++){
     3957                sumMass += dz[i]*d[i];
     3958                fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
     3959        }
     3960
     3961        #if defined(_HAVE_AD_)
     3962        /*we want to avoid the round operation at all cost. Not differentiable.*/
     3963        _error_("not implemented yet");
     3964        #else
     3965        dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
     3966        dMass = round(dMass * 100.0)/100.0;
     3967
     3968        /*Check mass conservation:*/
     3969        if (dMass != 0.0){
     3970                _printf_("total system mass not conserved in MB function \n");
     3971        }
     3972        #endif
     3973
     3974        /*Check bottom grid cell T is unchanged:*/
     3975        if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
     3976                if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
     3977        }
    39773978
    39783979        /*Save generated inputs: */
    … …  
    39953996        this->SetElementInput(SmbMeanLHFEnum,meanLHF);
    39963997        this->SetElementInput(SmbMeanSHFEnum,meanSHF);
    3997         this->SetElementInput(SmbDzAddEnum,sumdz_add);
     3998        this->SetElementInput(SmbDzAddEnum,dz_add);
     3999        this->SetElementInput(SmbMInitnum,initMass);
    39984000        this->SetElementInput(SmbMAddEnum,sumMassAdd/dt);
     4001        this->SetElementInput(SmbWAddEnum,sumW/dt);
    39994002        this->SetElementInput(SmbFACEnum,fac/1000.); // output in meters
     4003        this->SetElementInput(SmbECDtEnum,EC);
    40004004
    40014005        /*Free allocations:{{{*/
    … …  
    40164020        if(aini) xDelete<IssmDouble>(aini);
    40174021        if(Tini) xDelete<IssmDouble>(Tini);
     4022        if(swf) xDelete<IssmDouble>(swf);
    40184023
    40194024        delete gauss;

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

    r24565 r24683  
    172172                void               SmbSemic();
    173173                int                Sid();
    174                 void               SmbGemb();
     174                void               SmbGemb(IssmDouble timeinputs, int count);
    175175                void               StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input);
    176176                void               StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input,Input2* vz_input);

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

    r24509 r24683  
    66#include "../../shared/shared.h"
    77#include "../../toolkits/toolkits.h"
     8#include "../modules.h"
     9#include "../../classes/Inputs2/TransientInput2.h"
    810
    911const double Pi = 3.141592653589793;
    … …  
    3234void Gembx(FemModel* femmodel){  /*{{{*/
    3335
    34         for(int i=0;i<femmodel->elements->Size();i++){
    35         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    36         element->SmbGemb();
     36        int        count=0;
     37        IssmDouble time,dt,finaltime,starttime;
     38        IssmDouble timeclim=0.0;
     39        IssmDouble t,smb_dt;
     40   IssmDouble delta;
     41        bool       isclimatology=false;
     42
     43        femmodel->parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
     44   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
     45        femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
     46        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     47   femmodel->parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
     48        femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
     49
     50        //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
     51        //go back to time - deltaT:
     52        time-=dt;
     53
     54        IssmDouble timeinputs = time;
     55
     56        /*Start loop: */
     57        count=1;
     58        for (t=time;t<time+dt;t=t+smb_dt){
     59
     60                for(int i=0;i<femmodel->elements->Size();i++){
     61                        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     62
     63                        timeclim=time;
     64                        if (isclimatology){
     65                                //If this is a climatology, we need to repeat the forcing after the final time
     66                                TransientInput2* Ta_input_tr  = element->inputs2->GetTransientInput(SmbTaEnum);    _assert_(Ta_input_tr);
     67
     68                                /*Get temperature climatology value*/
     69                                int offsetend = Ta_input_tr->GetTimeInputOffset(finaltime);
     70                                IssmDouble time0     = Ta_input_tr->GetTimeByOffset(-1);
     71                                IssmDouble timeend   = Ta_input_tr->GetTimeByOffset(offsetend);
     72                                if (time>time0 & timeend>time0){
     73                                        delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
     74                                        timeclim=time0+delta;
     75                                }
     76                        }
     77                        timeinputs = t-time+timeclim;
     78                        element->SmbGemb(timeinputs,count);
     79                }
     80                count=count+1;
    3781        }
    3882
    … …  
    10451089                        // SWs and SWss coefficients need to be better constranted. Greuell
    10461090                        // and Konzelmann 1994 used SWs = 0.36 and SWss = 0.64 as this the
    1047                         // the // of SW radiation with wavelengths > and < 800 nm
     1091                        // the % of SW radiation with wavelengths > and < 800 nm
    10481092                        // respectively.  This, however, may not account for the fact that
    10491093                        // the albedo of wavelengths > 800 nm has a much lower albedo.
    … …  
    21072151
    21082152        // calculate the Bulk Richardson Number (Ri)
    2109         Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
    2110 
    2111         // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
     2153        Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
     2154
     2155        // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
    21122156
    21132157        // do not allow Ri to exceed 0.19

  • issm/trunk-jpl/src/c/shared/Enum/Enum.vim

    r24626 r24683  
    149149syn keyword cConstant FrictionCouplingEnum
    150150syn keyword cConstant FrictionDeltaEnum
     151syn keyword cConstant FrictionEffectivePressureLimitEnum
    151152syn keyword cConstant FrictionFEnum
    152153syn keyword cConstant FrictionGammaEnum
    … …  
    717718syn keyword cConstant SmbEAirEnum
    718719syn keyword cConstant SmbECEnum
     720syn keyword cConstant SmbECDtEnum
    719721syn keyword cConstant SmbECiniEnum
    720722syn keyword cConstant SmbElaEnum
    … …  
    774776syn keyword cConstant SmbVzEnum
    775777syn keyword cConstant SmbWEnum
     778syn keyword cConstant SmbWAddEnum
    776779syn keyword cConstant SmbWiniEnum
    777780syn keyword cConstant SmbZMaxEnum
    … …  
    13351338syn keyword cType Cfsurfacesquare
    13361339syn keyword cType Channel
    1337 syn keyword cType classes
    13381340syn keyword cType Constraint
    13391341syn keyword cType Constraints
    … …  
    13421344syn keyword cType ControlInput2
    13431345syn keyword cType Covertree
     1346syn keyword cType DataSetParam
    13441347syn keyword cType DatasetInput2
    1345 syn keyword cType DataSetParam
    13461348syn keyword cType Definition
    13471349syn keyword cType DependentObject
    … …  
    13551357syn keyword cType ElementInput2
    13561358syn keyword cType ElementMatrix
     1359syn keyword cType ElementVector
    13571360syn keyword cType Elements
    1358 syn keyword cType ElementVector
    13591361syn keyword cType ExponentialVariogram
    13601362syn keyword cType ExternalResult
    … …  
    13631365syn keyword cType Friction
    13641366syn keyword cType Gauss
    1365 syn keyword cType GaussianVariogram
    1366 syn keyword cType gaussobjects
    13671367syn keyword cType GaussPenta
    13681368syn keyword cType GaussSeg
    13691369syn keyword cType GaussTetra
    13701370syn keyword cType GaussTria
     1371syn keyword cType GaussianVariogram
    13711372syn keyword cType GenericExternalResult
    13721373syn keyword cType GenericOption
    … …  
    13831384syn keyword cType IssmDirectApplicInterface
    13841385syn keyword cType IssmParallelDirectApplicInterface
    1385 syn keyword cType krigingobjects
    13861386syn keyword cType Load
    13871387syn keyword cType Loads
    … …  
    13941394syn keyword cType Matice
    13951395syn keyword cType Matlitho
    1396 syn keyword cType matrixobjects
    13971396syn keyword cType MatrixParam
    13981397syn keyword cType Misfit
    … …  
    14071406syn keyword cType Observations
    14081407syn keyword cType Option
     1408syn keyword cType OptionUtilities
    14091409syn keyword cType Options
    1410 syn keyword cType OptionUtilities
    14111410syn keyword cType Param
    14121411syn keyword cType Parameters
    … …  
    14221421syn keyword cType Regionaloutput
    14231422syn keyword cType Results
     1423syn keyword cType RiftStruct
    14241424syn keyword cType Riftfront
    1425 syn keyword cType RiftStruct
    14261425syn keyword cType Seg
    14271426syn keyword cType SegInput2
     1427syn keyword cType SegRef
    14281428syn keyword cType Segment
    1429 syn keyword cType SegRef
    14301429syn keyword cType SpcDynamic
    14311430syn keyword cType SpcStatic
    … …  
    14461445syn keyword cType Vertex
    14471446syn keyword cType Vertices
     1447syn keyword cType classes
     1448syn keyword cType gaussobjects
     1449syn keyword cType krigingobjects
     1450syn keyword cType matrixobjects
    14481451syn keyword cType AdjointBalancethickness2Analysis
    14491452syn keyword cType AdjointBalancethicknessAnalysis
    … …  
    14641467syn keyword cType FreeSurfaceBaseAnalysis
    14651468syn keyword cType FreeSurfaceTopAnalysis
     1469syn keyword cType GLheightadvectionAnalysis
    14661470syn keyword cType GiaIvinsAnalysis
    1467 syn keyword cType GLheightadvectionAnalysis
    14681471syn keyword cType HydrologyDCEfficientAnalysis
    14691472syn keyword cType HydrologyDCInefficientAnalysis

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

    r24666 r24683  
    715715        SmbEAirEnum,
    716716        SmbECEnum,
     717        SmbECDtEnum,
    717718        SmbECiniEnum,
    718719        SmbElaEnum,
    … …  
    734735        SmbMeanULWEnum,
    735736        SmbMeltEnum,
     737        SmbMInitnum,
    736738        SmbMonthlytemperaturesEnum,
    737739        SmbNetLWEnum,
    … …  
    772774        SmbVzEnum,
    773775        SmbWEnum,
     776        SmbWAddEnum,
    774777        SmbWiniEnum,
    775778        SmbZMaxEnum,

  • issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

    r24666 r24683  
    720720                case SmbEAirEnum : return "SmbEAir";
    721721                case SmbECEnum : return "SmbEC";
     722                case SmbECDtEnum : return "SmbECDt";
    722723                case SmbECiniEnum : return "SmbECini";
    723724                case SmbElaEnum : return "SmbEla";
    … …  
    777778                case SmbVzEnum : return "SmbVz";
    778779                case SmbWEnum : return "SmbW";
     780                case SmbWAddEnum : return "SmbWAdd";
    779781                case SmbWiniEnum : return "SmbWini";
    780782                case SmbZMaxEnum : return "SmbZMax";

  • issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

    r24666 r24683  
    735735              else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
    736736              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
     737              else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    737738              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    738739              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    … …  
    751752              else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
    752753              else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
    753               else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
     757              if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
     758              else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
    758759              else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    759760              else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
    … …  
    795796              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    796797              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
     798              else if (strcmp(name,"SmbWAdd")==0) return SmbWAddEnum;
    797799              else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
    798800              else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
    … …  
    873875              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    874876              else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    875               else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    876               else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
     880              if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
     881              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
     882              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    881883              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    882884              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    … …  
    996998              else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum;
    997999              else if (strcmp(name,"IntInput2")==0) return IntInput2Enum;
    998               else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    999               else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
     1003              if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
     1004              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
     1005              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    10041006              else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    10051007              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    … …  
    11191121              else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
    11201122              else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
    1121               else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    1122               else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
     1126              if (strcmp(name,"Incremental")==0) return IncrementalEnum;
     1127              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
     1128              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    11271129              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
    11281130              else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum;
    … …  
    12421244              else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
    12431245              else if (strcmp(name,"Regular")==0) return RegularEnum;
    1244               else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
    1245               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
     1249              if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
     1250              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
     1251              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    12501252              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    12511253              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
Note: See TracChangeset for help on using the changeset viewer.