Changeset 22628


Ignore:
Timestamp:
03/26/18 11:06:37 (7 years ago)
Author:
tpelle
Message:

CHG: cleaning up

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

Legend:

Unmodified
Added
Removed

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

    r22625 r22628  
    27772777}
    27782778/*}}}*/
    2779 void       Tria::PicoUpdateBoxid(int* pmax_boxid_basin){/*{{{*/
     2779void       Tria::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/
    27802780
    27812781        if(!this->IsIceInElement() || !this->IsFloating()) return;
    27822782
    2783         int i,boxid,basin_id,boxid_max;
    2784         IssmDouble dist_gl,dist_cf,rel_dist_gl,lowbound,highbound;
     2783        int        basin_id;
     2784        IssmDouble dist_gl,dist_cf;
    27852785
    27862786        inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
    2787         boxid_max=pmax_boxid_basin[basin_id];
     2787        IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id]);
    27882788
    27892789        Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
    2790         Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input);
     2790        Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum);  _assert_(dist_cf_input);
    27912791
    27922792        /*Get dist_gl and dist_cf at center of element*/
    … …  
    27942794        dist_gl_input->GetInputValue(&dist_gl,gauss);
    27952795        dist_cf_input->GetInputValue(&dist_cf,gauss);
     2796        delete gauss;
    27962797       
    27972798        /*Ensure values are positive for floating ice*/
    2798         if(dist_gl<0){dist_gl=dist_gl*(-1);}
    2799         if(dist_cf<0){dist_cf=dist_cf*(-1);}
     2799        dist_gl = fabs(dist_gl);
     2800        dist_cf = fabs(dist_cf);
    28002801
    28012802        /*Compute relative distance to grounding line*/
    2802         rel_dist_gl=dist_gl/(dist_gl+dist_cf);
     2803        IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf);
    28032804
    28042805        /*Assign box numbers based on rel_dist_gl*/
    2805         boxid=-1;
    2806         for(i=0;i<boxid_max;i++){
    2807                 lowbound = 1-sqrt((boxid_max-i-0.0)/boxid_max);
    2808                 highbound = 1-sqrt((boxid_max-i-1.0)/boxid_max);
     2806        int boxid = -1;
     2807        for(IssmDouble i=0.;i<boxid_max;i++){
     2808                IssmDouble lowbound  = 1. -sqrt((boxid_max-i   )/boxid_max);
     2809                IssmDouble highbound = 1. -sqrt((boxid_max-i-1.)/boxid_max);
    28092810                if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){
    28102811                        boxid=i;
    … …  
    28152816
    28162817        this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));       
    2817 
    2818         /*Cleanup & return: */
    2819         delete gauss;
    2820 }
    2821 /*}}}*/
     2818}/*}}}*/
    28222819void       Tria::PicoUpdateFirstBox(){/*{{{*/
    28232820
    … …  
    28292826
    28302827        int basinid, maxbox, num_basins, numnodes, M;
    2831         IssmDouble time, rhoi, rhow, earth_grav, rho_star, nu, latentheat, c_p_ocean, lambda, a, b, c;
    2832         IssmDouble alpha, Beta, gamma_T, overturning_coeff;
    2833         IssmDouble thickness, toc_farocean, soc_farocean, area_box1;
    2834         IssmDouble pressure, T_star, g1, s1, p_coeff, q_coeff;
    2835         IssmDouble* boxareas=NULL;
    2836         IssmDouble* t_farocean=NULL;
    2837         IssmDouble* s_farocean=NULL;
     2828        IssmDouble time, gamma_T, overturning_coeff,thickness;
     2829        IssmDouble pressure, T_star,p_coeff, q_coeff,toc_farocean,soc_farocean;
     2830        IssmDouble* boxareas   = NULL;
    28382831
    28392832        /*Get variables*/
    2840         rhoi                    = this->GetMaterialParameter(MaterialsRhoIceEnum);
    2841         rhow                    = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    2842         earth_grav = this->GetMaterialParameter(ConstantsGEnum);
    2843         rho_star                = 1033.;             // kg/m^3
    2844         nu                              = rhoi/rhow;
    2845         latentheat      = this->GetMaterialParameter(MaterialsLatentheatEnum);
    2846         c_p_ocean       = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
    2847         lambda          = latentheat/c_p_ocean;
    2848         a                               = -0.0572;           // K/psu
    2849         b                               = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum);  //K
    2850         c           = 7.77e-4;
    2851         alpha                   = 7.5e-5;            // 1/K
    2852         Beta                    = 7.7e-4;            // K
     2833        IssmDouble rhoi       = this->GetMaterialParameter(MaterialsRhoIceEnum);
     2834        IssmDouble rhow       = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     2835        IssmDouble earth_grav = this->GetMaterialParameter(ConstantsGEnum);
     2836        IssmDouble rho_star   = 1033.;             // kg/m^3
     2837        IssmDouble nu         = rhoi/rhow;
     2838        IssmDouble latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum);
     2839        IssmDouble c_p_ocean  = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
     2840        IssmDouble lambda     = latentheat/c_p_ocean;
     2841        IssmDouble a          = -0.0572;           // K/psu
     2842        IssmDouble b          = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum);  //K
     2843        IssmDouble c          = 7.77e-4;
     2844        IssmDouble alpha      = 7.5e-5;            // 1/K
     2845        IssmDouble Beta       = 7.7e-4;            // K
    28532846
    28542847        /* Get parameters and inputs */
    … …  
    28572850        this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
    28582851        this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
    2859         this->parameters->FindParam(&t_farocean, &num_basins, time, BasalforcingsPicoFarOceantemperatureEnum);
    2860         this->parameters->FindParam(&s_farocean, &num_basins, time, BasalforcingsPicoFarOceansalinityEnum);
     2852        this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
     2853        this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
    28612854        this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
    2862         M = num_basins*maxbox;
    28632855        this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
    28642856        this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    2865        
    28662857        Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
    2867 
    2868         toc_farocean                                                    = t_farocean[basinid];
    2869         soc_farocean                                                    = s_farocean[basinid];
    2870         area_box1                                                               = boxareas[basinid*maxbox+boxid];
    2871         g1                                                                                      = area_box1*gamma_T;
    2872         s1                                                                                      = soc_farocean/(nu*lambda);
    2873 
     2858        IssmDouble area_box1    = boxareas[basinid*maxbox+boxid];
     2859        IssmDouble g1           = area_box1*gamma_T;
     2860        IssmDouble s1           = soc_farocean/(nu*lambda);
     2861
     2862        /*Define new inputs*/
    28742863   IssmDouble basalmeltrates_shelf[NUMVERTICES];
    28752864        IssmDouble potential_pressure_melting_point[NUMVERTICES];   
    … …  
    28802869        /* Start looping on the number of nodes and calculate ocean vars */
    28812870        Gauss* gauss=this->NewGauss();
    2882         for (int i=0;i<NUMVERTICES;i++){
     2871        for(int i=0;i<NUMVERTICES;i++){
    28832872                gauss->GaussVertex(i);
    28842873                thickness_input->GetInputValue(&thickness,gauss);
    28852874                pressure = (rhoi*earth_grav*1e-4)*thickness;
    2886                 T_star = a*soc_farocean+b-c*pressure-toc_farocean;
    2887                 p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
    2888                 q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
     2875                T_star   = a*soc_farocean+b-c*pressure-toc_farocean;
     2876                p_coeff  = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
     2877                q_coeff  = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
    28892878
    28902879                /* To avoid negatives under the square root */
    2891                 if ((0.25*pow(p_coeff,2)-q_coeff)<0){
    2892                         q_coeff = (0.25*pow(p_coeff,2));
     2880                if((0.25*pow(p_coeff,2)-q_coeff)<0){
     2881                        q_coeff = 0.25*p_coeff*p_coeff;
    28932882                }
    28942883
    … …  
    29072896        /*Cleanup and return */
    29082897        delete gauss;
    2909         xDelete<IssmDouble>(t_farocean);
    2910         xDelete<IssmDouble>(s_farocean);
    29112898        xDelete<IssmDouble>(boxareas);
    29122899}
    … …  
    29202907        if(loopboxid!=boxid) return;
    29212908
    2922         int basinid, maxbox, num_basins, numnodes,M;
    2923         IssmDouble rhoi, rhow, earth_grav, rho_star, nu, latentheat, c_p_ocean, lambda, a, b, c;
    2924         IssmDouble alpha, Beta, gamma_T, overturning_coeff;
    2925         IssmDouble thickness, toc_farocean, soc_farocean, area_boxi;
    2926         IssmDouble mean_toc,mean_soc,mean_overturning;
    2927         IssmDouble pressure, T_star, g1, g2, p_coeff, q_coeff;
    2928         IssmDouble* toc_weighted_avg=NULL;
    2929         IssmDouble* soc_weighted_avg=NULL;
    2930         IssmDouble* overturning_weighted_avg=NULL;
    2931         IssmDouble* boxareas=NULL;
     2909        int        basinid, maxbox, num_basins, numnodes,M;
     2910        IssmDouble gamma_T, overturning_coeff;
     2911        IssmDouble thickness, toc_farocean, soc_farocean;
     2912        IssmDouble pressure, T_star,p_coeff, q_coeff;
     2913        IssmDouble* toc_weighted_avg         = NULL;
     2914        IssmDouble* soc_weighted_avg         = NULL;
     2915        IssmDouble* overturning_weighted_avg = NULL;
     2916        IssmDouble* boxareas                 = NULL;
    29322917
    29332918        /*Get variables*/
    2934         rhoi        = this->GetMaterialParameter(MaterialsRhoIceEnum);
    2935         rhow        = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    2936         earth_grav = this->GetMaterialParameter(ConstantsGEnum);
    2937         rho_star    = 1033.;             // kg/m^3
    2938         nu          = rhoi/rhow;
    2939         latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum);
    2940         c_p_ocean   = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
    2941         lambda      = latentheat/c_p_ocean;
    2942         a           = -0.0572;          // K/psu
    2943         b           = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum);  //K
    2944         c           = 7.77e-4;
    2945         alpha       = 7.5e-5;           // 1/K
    2946         Beta        = 7.7e-4;           // K
     2919        IssmDouble rhoi       = this->GetMaterialParameter(MaterialsRhoIceEnum);
     2920        IssmDouble rhow       = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
     2921        IssmDouble earth_grav = this->GetMaterialParameter(ConstantsGEnum);
     2922        IssmDouble rho_star   = 1033.;             // kg/m^3
     2923        IssmDouble nu         = rhoi/rhow;
     2924        IssmDouble latentheat = this->GetMaterialParameter(MaterialsLatentheatEnum);
     2925        IssmDouble c_p_ocean  = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
     2926        IssmDouble lambda     = latentheat/c_p_ocean;
     2927        IssmDouble a          = -0.0572;          // K/psu
     2928        IssmDouble b          = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum);  //K
     2929        IssmDouble c          = 7.77e-4;
     2930        IssmDouble alpha      = 7.5e-5;           // 1/K
     2931        IssmDouble Beta       = 7.7e-4;           // K
    29472932
    29482933        /* Get parameters and inputs */
    … …  
    29542939        this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);
    29552940        this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);
    2956         M=num_basins*maxbox;
    29572941        this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
    29582942        this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    2959 
    29602943   Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
    29612944
    2962         area_boxi                                                               = boxareas[basinid*maxbox+boxid];
    2963    mean_toc                                                                     = toc_weighted_avg[basinid];
    2964         mean_soc                                                                        = soc_weighted_avg[basinid];
    2965         mean_overturning                                                = overturning_weighted_avg[basinid];
    2966         g1                                                                                      = area_boxi*gamma_T;
    2967    g2                                                                                   = g2/(nu*lambda);
     2945        _assert_(basinid<=num_basins);
     2946        IssmDouble area_boxi        = boxareas[basinid*maxbox+boxid];
     2947        IssmDouble mean_toc         = toc_weighted_avg[basinid];
     2948        IssmDouble mean_soc         = soc_weighted_avg[basinid];
     2949        IssmDouble mean_overturning = overturning_weighted_avg[basinid];
     2950        IssmDouble g1               = area_boxi*gamma_T;
     2951        IssmDouble g2               = g2/(nu*lambda);
    29682952
    29692953        IssmDouble basalmeltrates_shelf[NUMVERTICES];
    … …  
    29742958        /* Start looping on the number of nodes and calculate ocean vars */
    29752959        Gauss* gauss=this->NewGauss();
    2976         for (int i=0;i<NUMVERTICES;i++){
     2960        for(int i=0;i<NUMVERTICES;i++){
    29772961                gauss->GaussVertex(i);
    29782962                thickness_input->GetInputValue(&thickness,gauss);
    2979                 pressure = (rhoi*earth_grav*1e-4)*thickness;
    2980                 T_star = a*mean_soc+b-c*pressure-mean_toc;
    2981                 Tocs[i] = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc));
    2982                 Socs[i] = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda));
     2963                pressure = (rhoi*earth_grav*1.e-4)*thickness;
     2964                T_star   = a*mean_soc+b-c*pressure-mean_toc;
     2965                Tocs[i]  = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc));
     2966                Socs[i]  = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda));
    29832967                potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
    29842968                basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);

  • issm/trunk-jpl/src/c/classes/Params/BoolParam.h

    r22625 r22628  
    4545                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4646                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    47                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4847                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4948                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/DataSetParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(FILE** pfile){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a file pointer");}
    5049                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}

  • issm/trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.h

    r22625 r22628  
    4848                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << "cannot return a IssmDouble");}
    4949                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    50                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    5150                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << "cannot return a string");}
    5251                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << "cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/DoubleMatParam.h

    r22625 r22628  
    4848                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4949                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    50                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    5150                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5251                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/DoubleParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){*pIssmDouble=value;}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5049                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.cpp

    r20827 r22628  
    7777
    7878/*DoubleVecParam virtual functions definitions: */
    79 void  DoubleVecParam::GetParameterValue(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/
    80         IssmDouble* output=NULL;
    81         int M;
     79void  DoubleVecParam::GetParameterValue(IssmDouble** poutput,int* pM){/*{{{*/
    8280
    83         M=this->M;
    84         output=xNew<IssmDouble>(M);
     81        IssmDouble* output=xNew<IssmDouble>(M);
    8582        xMemCpy<IssmDouble>(output,values,M);
    8683
    8784        /*Assign output pointers:*/
    8885        if(pM) *pM=M;
    89         *pIssmDoublearray=output;
     86        *poutput=output;
    9087}
    9188/*}}}*/
    92 void  DoubleVecParam::GetParameterValue(IssmDouble** pIssmDoublearray,int* pM,int* pN){/*{{{*/
    93         IssmDouble* output=NULL;
    94         int M;
    95         int N;
     89void  DoubleVecParam::GetParameterValue(IssmDouble** poutput,int* pM,int* pN){/*{{{*/
    9690
    97         N=1;
    98         M=this->M;
    99         output=xNew<IssmDouble>(M);
     91        IssmDouble* output=xNew<IssmDouble>(this->M);
    10092        xMemCpy<IssmDouble>(output,values,M);
    10193
    10294        /*Assign output pointers:*/
    103         if(pM) *pM=M;
    104         if(pN) *pN=N;
    105         *pIssmDoublearray=output;
     95        if(pM) *pM=this->M;
     96        if(pN) *pN=1;
     97        *poutput=output;
    10698}
    10799/*}}}*/

  • issm/trunk-jpl/src/c/classes/Params/DoubleVecParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5049                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/FileParam.h

    r22625 r22628  
    4545                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4646                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    47                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4847                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4948                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/GenericParam.h

    r22625 r22628  
    7171                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble");}
    7272                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDouble for a given time");}
    73                                          void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a IssmDoubleArray for a given time");}
    7473                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string");}
    7574                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/IntMatParam.h

    r22625 r22628  
    4747                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4848                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    49                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    5049                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5150                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/IntParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5049                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/IntVecParam.h

    r22625 r22628  
    4747                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4848                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    49                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    5049                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5150                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/MatrixParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5049                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/Param.h

    r22625 r22628  
    3434                virtual void  GetParameterValue(IssmDouble* pIssmDouble)=0;
    3535                virtual void  GetParameterValue(IssmDouble* pdouble,IssmDouble time)=0;
    36                 virtual void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time)=0;
     36                virtual void  GetParameterValue(IssmDouble* pdouble,int row, IssmDouble time){_error_("not implemented yet");};
    3737                virtual void  GetParameterValue(char** pstring)=0;
    3838                virtual void  GetParameterValue(char*** pstringarray,int* pM)=0;

  • issm/trunk-jpl/src/c/classes/Params/Parameters.cpp

    r22625 r22628  
    300300}
    301301/*}}}*/
    302 void Parameters::FindParam(IssmDouble** pIssmDoublearray, int* pM, IssmDouble time, int param_enum){ _assert_(this);/*{{{*/
    303 
    304         _assert_(param_enum>ParametersSTARTEnum);
    305         _assert_(param_enum<ParametersENDEnum);
    306 
    307         int index = param_enum - ParametersSTARTEnum -1;
    308         if(!this->params[index]) _error_("Parameter " << EnumToStringx(param_enum) <<" not set");
    309         this->params[index]->GetParameterValue(pIssmDoublearray, pM, time);
     302void Parameters::FindParam(IssmDouble* pscalar,int row,IssmDouble time, int param_enum){ _assert_(this);/*{{{*/
     303
     304        _assert_(param_enum>ParametersSTARTEnum);
     305        _assert_(param_enum<ParametersENDEnum);
     306
     307        int index = param_enum - ParametersSTARTEnum -1;
     308        if(!this->params[index]) _error_("Parameter " << EnumToStringx(param_enum) <<" not set");
     309        this->params[index]->GetParameterValue(pscalar,row,time);
    310310}
    311311/*}}}*/

  • issm/trunk-jpl/src/c/classes/Params/Parameters.h

    r22625 r22628  
    4040                void  FindParam(IssmDouble* pscalar, int enum_type);
    4141                void  FindParam(IssmDouble* pscalar, int enum_type,IssmDouble time);
    42                 void    FindParam(IssmDouble** pIssmDoublearray, int* pM, IssmDouble time, int enum_type);
     42                void  FindParam(IssmDouble* pscalar, int row,IssmDouble time,int enum_type);
    4343                void  FindParam(char** pstring,int enum_type);
    4444                void  FindParam(char*** pstringarray,int* pM,int enum_type);

  • issm/trunk-jpl/src/c/classes/Params/StringArrayParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5049                void  GetParameterValue(char*** pstringarray,int* pM);

  • issm/trunk-jpl/src/c/classes/Params/StringParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(char** pstring);
    5049                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.cpp

    r22625 r22628  
    2525
    2626        enum_type=in_enum_type;
    27         N=in_N; //Number of time steps
    28         M=in_M; //Number of rows
     27        M=in_M; //Number of time steps
     28        N=in_N; //Number of rows
    2929        interpolation=interpolation_on;
    3030
    31         values=xNew<IssmDouble>(N*M);
    32         xMemCpy<IssmDouble>(values,in_values,N*M);
     31        values=xNew<IssmDouble>(M*N);
     32        xMemCpy<IssmDouble>(values,in_values,M*N);
    3333
    34         timesteps=xNew<IssmDouble>(N);
    35         xMemCpy<IssmDouble>(timesteps,in_time,N);
     34        timesteps=xNew<IssmDouble>(M);
     35        xMemCpy<IssmDouble>(timesteps,in_time,M);
    3636}
    3737/*}}}*/
    … …  
    3939        xDelete<IssmDouble>(values);
    4040        xDelete<IssmDouble>(timesteps);
    41 }
    42 /*}}}*/
     41}/*}}}*/
    4342
    4443/*Object virtual functions definitions:*/
    4544Param* TransientArrayParam::copy() {/*{{{*/
    4645
    47         return new TransientArrayParam(this->enum_type,this->values,this->timesteps,this->interpolation,this->N,this->M);
     46        return new TransientArrayParam(this->enum_type,this->values,this->timesteps,this->interpolation,this->M,this->N);
    4847
    4948}
    … …  
    6867        _printf_("TransientArrayParam:\n");
    6968        _printf_("   enum: " << this->enum_type << " (" << EnumToStringx(this->enum_type) << ")\n");
    70         _printf_("   size: " << this->N << " by " << this->M << "\n");
     69        _printf_("   size: " << this->M << " by " << this->N << "\n");
    7170
    7271}
    … …  
    8382        MARSHALLING(N);
    8483        if(marshall_direction==MARSHALLING_BACKWARD){
    85                 values=xNew<IssmDouble>(N*M);
    86                 timesteps=xNew<IssmDouble>(N*M);
     84                values    = xNew<IssmDouble>(M*N);
     85                timesteps = xNew<IssmDouble>(N);
    8786        }
    88         MARSHALLING_ARRAY(values,IssmDouble,N*M);
    89         MARSHALLING_ARRAY(timesteps,IssmDouble,N*M);
     87        MARSHALLING_ARRAY(values,IssmDouble,M*N);
     88        MARSHALLING_ARRAY(timesteps,IssmDouble,N);
    9089
    9190}/*}}}*/
    … …  
    9493        return TransientArrayParamEnum;
    9594
    96 }
    97 /*}}}*/
     95}/*}}}*/
    9896
    9997/*TransientArrayParam virtual functions definitions: */
    100 void  TransientArrayParam::GetParameterValue(IssmDouble* pdouble,IssmDouble time, int row){/*{{{*/
     98void  TransientArrayParam::GetParameterValue(IssmDouble* pdouble,int row,IssmDouble time){/*{{{*/
    10199
    102100        IssmDouble output;
    103         bool   found;
     101        bool       found;
    104102
    105103        /*Ok, we have the time, go through the timesteps, and figure out which interval we
    … …  
    142140}
    143141/*}}}*/
    144 void  TransientArrayParam::GetParameterValue(IssmDouble** pIssmDoubleArray, int* pM, IssmDouble time){/*{{{*/
    145 
    146         IssmDouble* output=NULL;
    147         bool   found;
    148         int M=*pM;
    149 
    150         output=xNew<IssmDouble>(M);
    151 
    152         /*Ok, we have the time, go through the timesteps, and figure out which interval we
    153          *fall within. Then interpolate the values on this interval: */
    154         if(time<this->timesteps[0]){
    155                 /*get values for the first time: */
    156                 for(int k=0; k<M;k++){
    157                         output[k]=this->values[k*this->N];
    158                         found=true;
    159                 }
    160         }
    161         else if(time>this->timesteps[this->N-1]){
    162                 /*get values for the last time: */
    163                 for(int k=0; k<M;k++){
    164                         output[k]=this->values[(k+1)*this->N-1];
    165                         found=true;
    166                 }
    167         }
    168         else{
    169                 /*Find which interval we fall within: */
    170                 for(int i=0;i<this->N;i++){
    171                         if(time==this->timesteps[i]){
    172                                 /*We are right on one step time: */
    173                                 for(int k=0; k<M;k++){
    174                                         output[k]=this->values[k*this->N+i];
    175                                 }
    176                                 found=true;
    177                                 break; //we are done with the time interpolation.
    178                         }
    179                         else{
    180                                 if(this->timesteps[i]<time && time<this->timesteps[i+1]){
    181                                         /*ok, we have the interval ]i:i+1[. Interpolate linearly for now: */
    182                                         IssmDouble deltat=this->timesteps[i+1]-this->timesteps[i];
    183                                         IssmDouble alpha=(time-this->timesteps[i])/deltat;
    184                                         for(int k=0; k<M;k++){
    185                                                 if(interpolation==true) output[k]=(1.0-alpha)*this->values[k*this->N+i] + alpha*this->values[k*this->N+i+1];
    186                                                 else output[k]=this->values[k*this->N+i];
    187                                         }
    188                                         found=true;
    189                                         break;
    190                                 }
    191                                 else continue; //keep looking on the next interval
    192                         }
    193                 }
    194         }
    195         if(!found)_error_("did not find time interval on which to interpolate values");
    196         *pIssmDoubleArray=output;
    197 
    198 }/*}}}*/

  • issm/trunk-jpl/src/c/classes/Params/TransientArrayParam.h

    r22625 r22628  
    4848                void  GetParameterValue(int** pintarray,int* pM,int* pN){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a array of integers");}
    4949                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    50                 void  GetParameterValue(IssmDouble* pdouble,IssmDouble time,int row);
    51                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time);
     50                void  GetParameterValue(IssmDouble* pdouble,int row,IssmDouble time);
    5251                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Parameter " <<EnumToStringx(enum_type) << " needs row to be specified");}
    5352                void  GetParameterValue(char** pstring){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a string");}

  • issm/trunk-jpl/src/c/classes/Params/TransientParam.h

    r22625 r22628  
    4848                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4949                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time);
    50                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    5150                void  GetParameterValue(char** pstring){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a string");}
    5251                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/classes/Params/VectorParam.h

    r22625 r22628  
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    4747                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    48                 void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDoubleArray for a given time");}
    4948                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5049                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r22625 r22628  
    209209                                parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.overturning_coeff",BasalforcingsPicoOverturningCoeffEnum));
    210210                                parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_T",BasalforcingsPicoGammaTEnum));
    211                                 iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.farocean_temperature");
    212                                 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,&transparam[0],&transparam[M*(N-1)],interp,M,N));
     211                                iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature");
     212                                parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp,M,N));
    213213                                xDelete<IssmDouble>(transparam);
    214                                 iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.farocean_salinity");
    215                                 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,&transparam[0],&transparam[M*(N-1)],interp,M,N));
     214                                iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_salinity");
     215                                parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp,M,N));
    216216                                xDelete<IssmDouble>(transparam);
    217217                        break;
Note: See TracChangeset for help on using the changeset viewer.