Changeset 22625


Ignore:
Timestamp:
03/26/18 09:55:52 (7 years ago)
Author:
tpelle
Message:

NEW: added pico ocean melt parameterization and its nightly run test

Location:
issm/trunk-jpl
Files:
5 added
37 edited

Legend:

Unmodified
Added
Removed

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

    r22612 r22625  
    223223                                        ./modules/InputUpdateFromVectorx/InputUpdateFromVectorx.cpp\
    224224                                        ./modules/FloatingiceMeltingRatex/FloatingiceMeltingRatex.cpp\
     225                                        ./modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp\
    225226                                        ./modules/ConfigureObjectsx/ConfigureObjectsx.cpp\
    226227                                        ./modules/SpcNodesx/SpcNodesx.cpp\

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

    r22576 r22625  
    23412341}
    23422342/*}}}*/
    2343 void       Element::PicoFloatingiceMeltingRate(){/*{{{*/
    2344 
    2345         int numvertices      = this->GetNumberOfVertices();
    2346         IssmDouble  deepwaterel,upperwaterel,deepwatermelt;
    2347         IssmDouble* base     = xNew<IssmDouble>(numvertices);
    2348         IssmDouble* values   = xNew<IssmDouble>(numvertices);
    2349         IssmDouble time;
    2350 
    2351         parameters->FindParam(&time,TimeEnum);
    2352         parameters->FindParam(&deepwaterel,BasalforcingsDeepwaterElevationEnum,time);
    2353         parameters->FindParam(&upperwaterel,BasalforcingsUpperwaterElevationEnum,time);
    2354         parameters->FindParam(&deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum,time);
    2355 
    2356         this->GetInputListOnVertices(base,BaseEnum);
    2357         for(int i=0;i<numvertices;i++){
    2358                 if(base[i]>upperwaterel)      values[i]=0;
    2359                 else if (base[i]<deepwaterel) values[i]=deepwatermelt;
    2360                 else values[i]=deepwatermelt*(base[i]-upperwaterel)/(deepwaterel-upperwaterel);
    2361         }
    2362 
    2363         this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
    2364         xDelete<IssmDouble>(base);
    2365         xDelete<IssmDouble>(values);
    2366 
    2367 }/*}}}*/
    23682343void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
    23692344

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

    r22515 r22625  
    142142                ElementMatrix*     NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
    143143                ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
    144                 void                                     PicoFloatingiceMeltingRate();
    145144                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
    146145                IssmDouble         PureIceEnthalpy(IssmDouble pressure);
    … …  
    213212                virtual Element*   GetBasalElement(void)=0;
    214213                virtual int        GetElementType(void)=0;
     214                virtual IssmDouble GetHorizontalSurfaceArea(void){_error_("not implemented");};
    215215                virtual void       GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
    216216                virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
    … …  
    278278                virtual int        NumberofNodesVelocity(void)=0;
    279279                virtual void       PicoUpdateBoxid(int* pmax_boxid_basin)=0;
    280                 virtual void       PicoUpdateFirstBox(void){_error_("not implemented");};
     280                virtual void       PicoUpdateFirstBox()=0;
     281                virtual void       PicoUpdateNextBox(int loopboxid)=0;
    281282                virtual void       PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
    282283                virtual int        PressureInterpolation()=0;

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

    r22471 r22625  
    140140                int            NumberofNodesVelocity(void);
    141141                void                            PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};
     142                void                            PicoUpdateFirstBox(){_error_("not implemented yet");};
     143                void                            PicoUpdateNextBox(int loopboxid){_error_("not implemented yet");};
    142144                void           PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    143145                int            PressureInterpolation();

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

    r22471 r22625  
    129129                int         NumberofNodesVelocity(void){_error_("not implemented yet");};
    130130                void        PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};
     131                void                    PicoUpdateFirstBox(){_error_("not implemented yet");};
     132                void                    PicoUpdateNextBox(int loopboxid){_error_("not implemented yet");};
    131133                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
    132134                int         PressureInterpolation(void){_error_("not implemented yet");};

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

    r22471 r22625  
    137137                int         NumberofNodesVelocity(void);
    138138                void                    PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};       
     139                void                    PicoUpdateFirstBox(){_error_("not implemented yet");};
     140                void                    PicoUpdateNextBox(int loopboxid){_error_("not implemented yet");};
    139141                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
    140142                int         PressureInterpolation(void);

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

    r22515 r22625  
    11071107        _assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
    11081108        return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
     1109}
     1110/*}}}*/
     1111IssmDouble Tria::GetHorizontalSurfaceArea(void){/*{{{*/
     1112
     1113        return this->GetArea();
    11091114}
    11101115/*}}}*/
    … …  
    27762781        if(!this->IsIceInElement() || !this->IsFloating()) return;
    27772782
    2778         //Load inputs and params
    2779         int i, boxid;
    2780         int basin_id, boxid_max;
    2781         IssmDouble dist_gl;
    2782         IssmDouble dist_cf;
    2783         IssmDouble rel_dist_gl;
     2783        int i,boxid,basin_id,boxid_max;
     2784        IssmDouble dist_gl,dist_cf,rel_dist_gl,lowbound,highbound;
    27842785
    27852786        inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
    2786         boxid_max=pmax_boxid_basin[basin_id]; //0-based
    2787 
    2788         Input* dist_gl_input=NULL;
    2789         Input* dist_cf_input=NULL;
    2790         dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
    2791         dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input);
    2792 
     2787        boxid_max=pmax_boxid_basin[basin_id];
     2788
     2789        Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
     2790        Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum); _assert_(dist_cf_input);
     2791
     2792        /*Get dist_gl and dist_cf at center of element*/
    27932793        Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
    27942794        dist_gl_input->GetInputValue(&dist_gl,gauss);
    27952795        dist_cf_input->GetInputValue(&dist_cf,gauss);
    2796 
     2796       
     2797        /*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);}
     2800
     2801        /*Compute relative distance to grounding line*/
    27972802        rel_dist_gl=dist_gl/(dist_gl+dist_cf);
    27982803
     2804        /*Assign box numbers based on rel_dist_gl*/
    27992805        boxid=-1;
    28002806        for(i=0;i<boxid_max;i++){
    2801                 if(rel_dist_gl>=(1-sqrt((boxid_max-i)/boxid_max)) && rel_dist_gl<=(1-sqrt((boxid_max-i-1)/boxid_max))){
    2802                         boxid=i; break;
    2803                 }
    2804         }
    2805         if(boxid==-1){_error_("no boxid found for element" << this->Sid());}
     2807                lowbound = 1-sqrt((boxid_max-i-0.0)/boxid_max);
     2808                highbound = 1-sqrt((boxid_max-i-1.0)/boxid_max);
     2809                if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){
     2810                        boxid=i;
     2811                        break;
     2812                }
     2813        }
     2814        if(boxid==-1){_error_("No boxid found for element " << this->Sid() << "!");}
    28062815
    28072816        this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));       
    … …  
    28132822void       Tria::PicoUpdateFirstBox(){/*{{{*/
    28142823
    2815         IssmDouble rhoi, rhow, earth_grav, rho_star, nu, latentheat, c_p_ocean, lambda, a, b, c;
    2816         IssmDouble alpha, Beta, gamma_T, overturning_coeff, t_farocean, s_farocean;
    2817         IssmDouble basinid, boxid, thickness, toc_farocean, soc_farocean, area_box1;
    2818         IssmDouble pressure, T_star, g1, s1, p_coeff, q_coeff, Toc, Soc, potential_pressure_melting_point;
    2819         IssmDouble basalmeltrate_shelf, overturning, maxbox;
    2820 
    2821         //Get variables
     2824        if(!this->IsIceInElement() || !this->IsFloating()) return;
     2825
     2826        int boxid;
     2827        this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
     2828        if(!boxid==0) return;
     2829
     2830        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;
     2838
     2839        /*Get variables*/
    28222840        rhoi                    = this->GetMaterialParameter(MaterialsRhoIceEnum);
    28232841        rhow                    = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    28242842        earth_grav  = this->GetMaterialParameter(ConstantsGEnum);
    2825         rho_star                = 1033.;             //kg/m^3
     2843        rho_star                = 1033.;             // kg/m^3
    28262844        nu                              = rhoi/rhow;
    28272845        latentheat      = this->GetMaterialParameter(MaterialsLatentheatEnum);
    2828         c_p_ocean       = this->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum);
     2846        c_p_ocean       = this->GetMaterialParameter(MaterialsMixedLayerCapacityEnum);
    28292847        lambda          = latentheat/c_p_ocean;
    2830         a                               = -0.0572;          //K/psu
     2848        a                               = -0.0572;           // K/psu
    28312849        b                               = 0.0788 + this->GetMaterialParameter(MaterialsMeltingpointEnum);  //K
    2832         c                               = this->GetMaterialParameter(MaterialsBetaEnum);
    2833         alpha                   = 0.000075;         //1/K
    2834         Beta                    = 0.00077;          //K
    2835 
     2850        c           = 7.77e-4;
     2851        alpha                   = 7.5e-5;            // 1/K
     2852        Beta                    = 7.7e-4;            // K
     2853
     2854        /* Get parameters and inputs */
     2855        this->parameters->FindParam(&time,TimeEnum);
     2856        this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);
    28362857        this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
    28372858        this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
    2838         this->parameters->FindParam(&t_farocean,BasalforcingsPicoFarOceansalinityEnum);
    2839         this->parameters->FindParam(&s_farocean,BasalforcingsPicoFarOceansalinityEnum);
    2840 
     2859        this->parameters->FindParam(&t_farocean, &num_basins, time, BasalforcingsPicoFarOceantemperatureEnum);
     2860        this->parameters->FindParam(&s_farocean, &num_basins, time, BasalforcingsPicoFarOceansalinityEnum);
     2861        this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
     2862        M = num_basins*maxbox;
     2863        this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
    28412864        this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
     2865       
     2866        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
     2874   IssmDouble basalmeltrates_shelf[NUMVERTICES];
     2875        IssmDouble potential_pressure_melting_point[NUMVERTICES];   
     2876        IssmDouble Tocs[NUMVERTICES];                                                   
     2877   IssmDouble Socs[NUMVERTICES];
     2878   IssmDouble overturnings[NUMVERTICES];
     2879
     2880        /* Start looping on the number of nodes and calculate ocean vars */
     2881        Gauss* gauss=this->NewGauss();
     2882        for (int i=0;i<NUMVERTICES;i++){
     2883                gauss->GaussVertex(i);
     2884                thickness_input->GetInputValue(&thickness,gauss);
     2885                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)));
     2889
     2890                /* 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));
     2893                }
     2894
     2895                Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff));
     2896                Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]);
     2897                potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
     2898                basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
     2899                overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i]));
     2900        }
     2901
     2902        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
     2903        this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
     2904        this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
     2905        this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum);
     2906
     2907        /*Cleanup and return */
     2908        delete gauss;
     2909        xDelete<IssmDouble>(t_farocean);
     2910        xDelete<IssmDouble>(s_farocean);
     2911        xDelete<IssmDouble>(boxareas);
     2912}
     2913/*}}}*/
     2914void       Tria::PicoUpdateNextBox(int loopboxid){/*{{{*/       
     2915
     2916        if(!this->IsIceInElement() || !this->IsFloating()) return;
     2917       
     2918        int boxid;
    28422919        this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
    2843         this->inputs->GetInputValue(&maxbox,BasalforcingsPicoMaxboxcountEnum);
    2844         this->inputs->GetInputValue(&thickness,ThicknessEnum);
    2845 
    2846         _error_("to be continued");
    2847         /*
    2848 
    2849         toc_farocean = t_farocean[basinid];
    2850         soc_farocean = s_farocean[basinid];
    2851         area_box1 = areas[basinid*maxbox+boxid];
    2852         pressure = (rhoi*earth_grav)*thickness;
    2853         T_star = a*soc_farocean+b-c*pressure-toc_farocean;
    2854         g1 = area_box1*gamma_T;
    2855         s1 = soc_farocean/(nu*lambda);
    2856         p_coeff = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
    2857         q_coeff = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
    2858         if ((0.25*(p_coeff)^(2) - q_coeff) < 0)
    2859           {q_coeff = (0.25*(p_coeff)^(2))}
    2860         Toc = toc_farocean-(-0.5*p_coeff+sqrt(0.25*p_coeff^(2)-q_coeff));
    2861         Soc = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Toc);
    2862         potential_pressure_melting_point = a*Soc+b-c*pressure;
    2863         basalmeltrate_shelf = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point-Toc);
    2864         overturning = overturning_coeff*rho_star*(Beta*(soc_farocean-Soc)-alpha*(toc_farocean-Toc));
    2865 
    2866         this->AddInput(BasalforcingsGroundediceMeltingRateEnum,basalmeltingrates,this->GetElementType());
    2867         */
    2868 
     2920        if(loopboxid!=boxid) return;
     2921
     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;
     2932
     2933        /*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
     2947
     2948        /* Get parameters and inputs */
     2949        this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);
     2950        this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
     2951        this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
     2952        this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
     2953        this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum);
     2954        this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);
     2955        this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);
     2956        M=num_basins*maxbox;
     2957        this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
     2958        this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
     2959
     2960   Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
     2961
     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);
     2968
     2969        IssmDouble basalmeltrates_shelf[NUMVERTICES];
     2970        IssmDouble potential_pressure_melting_point[NUMVERTICES];   
     2971        IssmDouble Tocs[NUMVERTICES];                                                   
     2972   IssmDouble Socs[NUMVERTICES];
     2973
     2974        /* Start looping on the number of nodes and calculate ocean vars */
     2975        Gauss* gauss=this->NewGauss();
     2976        for (int i=0;i<NUMVERTICES;i++){
     2977                gauss->GaussVertex(i);
     2978                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));
     2983                potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
     2984                basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
     2985        }
     2986
     2987        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
     2988        this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
     2989        this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
     2990
     2991        /*Cleanup and return */
     2992        delete gauss;
     2993        xDelete<IssmDouble>(boxareas);
     2994        xDelete<IssmDouble>(toc_weighted_avg);
     2995        xDelete<IssmDouble>(soc_weighted_avg);
     2996        xDelete<IssmDouble>(overturning_weighted_avg);
    28692997}
    28702998/*}}}*/
    … …  
    37203848        if(minvalue>fieldvalue || maxvalue<fieldvalue) return;
    37213849
     3850        /* check #2: If only one vertex is on fieldvalue, there is no segment here */
     3851        IssmDouble lsf[NUMVERTICES];
     3852        this->GetInputListOnVertices(&lsf[0],fieldenum);
     3853        int nrice=0;       
     3854        for(int i=0;i<NUMVERTICES;i++) if(lsf[i]==fieldvalue) nrice++;
     3855        if(nrice==1) return;
    37223856
    37233857        /*2. Write segments*/

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

    r22515 r22625  
    111111                int         NumberofNodesVelocity(void);
    112112                void        PicoUpdateBoxid(int* pmax_boxid_basin);
    113                 void        PicoUpdateFirstBox(void);
     113                void        PicoUpdateFirstBox();
     114                void        PicoUpdateNextBox(int loopboxid);
    114115                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    115116                int         PressureInterpolation();
    … …  
    165166                void           DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,IoModel* iomodel,int input_enum);
    166167                IssmDouble     GetArea(void);
     168                IssmDouble     GetHorizontalSurfaceArea(void);
    167169                IssmDouble         GetArea3D(void);
    168170                IssmDouble         GetAreaIce(void);

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

    r20827 r22625  
    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");}
    4748                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4849                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

    r22588 r22625  
    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");}
    4849                void  GetParameterValue(FILE** pfile){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a file pointer");}
    4950                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}

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

    r20827 r22625  
    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");}
    5051                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << "cannot return a string");}
    5152                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

    r20690 r22625  
    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");}
    5051                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5152                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

    r20827 r22625  
    4545                void  GetParameterValue(int** pintarray,int* pM,int* pN);
    4646                void  GetParameterValue(IssmDouble* pIssmDouble){*pIssmDouble=value;}
     47                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");}
    4749                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    48                 void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a IssmDouble for a given time");}
    4950                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}
    5051                void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM);

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

    r20827 r22625  
    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");}
    4849                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4950                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

    r20827 r22625  
    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");}
    4748                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4849                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

    r20827 r22625  
    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");}
    7374                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string");}
    7475                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(myEnumVal) << " cannot return a string array");}

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

    r20827 r22625  
    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");}
    4950                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5051                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

    r20827 r22625  
    3939                int   ObjectEnum();
    4040                /*}}}*/
    41                 /*Param vritual function definitions: {{{*/
     41                /*Param virtual function definitions: {{{*/
    4242                void  GetParameterValue(bool* pbool){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a bool");}
    4343                void  GetParameterValue(int* pinteger){*pinteger=value;}
    … …  
    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");}
    4849                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4950                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

    r20827 r22625  
    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");}
    4950                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    5051                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

    r20827 r22625  
    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");}
    4849                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4950                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

    r22617 r22625  
    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;
    3637                virtual void  GetParameterValue(char** pstring)=0;
    3738                virtual void  GetParameterValue(char*** pstringarray,int* pM)=0;

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

    r22617 r22625  
    300300}
    301301/*}}}*/
     302void 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);
     310}
     311/*}}}*/
    302312void Parameters::FindParam(char** pstring,int param_enum){ _assert_(this);/*{{{*/
    303313
    … …  
    349359
    350360        int index = param_enum - ParametersSTARTEnum -1;
     361
    351362        if(!this->params[index]) _error_("Parameter " << EnumToStringx(param_enum) <<" not set");
    352363        this->params[index]->GetParameterValue(pIssmDoublearray,pM);
     364
    353365}
    354366/*}}}*/

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

    r22617 r22625  
    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);
    4243                void  FindParam(char** pstring,int enum_type);
    4344                void  FindParam(char*** pstringarray,int* pM,int enum_type);

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

    r20827 r22625  
    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");}
    4849                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4950                void  GetParameterValue(char*** pstringarray,int* pM);

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

    r20827 r22625  
    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");}
    4849                void  GetParameterValue(char** pstring);
    4950                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

    r22471 r22625  
    142142}
    143143/*}}}*/
     144void  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

    r22471 r22625  
    4949                void  GetParameterValue(IssmDouble* pIssmDouble){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a IssmDouble");}
    5050                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time,int row);
     51                void  GetParameterValue(IssmDouble** pIssmDoublearray,int* pM, IssmDouble time);
    5152                void  GetParameterValue(IssmDouble* pdouble,IssmDouble time){_error_("Parameter " <<EnumToStringx(enum_type) << " needs row to be specified");}
    5253                void  GetParameterValue(char** pstring){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a string");}

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

    r20827 r22625  
    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");}
    5051                void  GetParameterValue(char** pstring){_error_("Parameter " <<EnumToStringx(enum_type) << " cannot return a string");}
    5152                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

    r20827 r22625  
    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");}
    4849                void  GetParameterValue(char** pstring){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string");}
    4950                void  GetParameterValue(char*** pstringarray,int* pM){_error_("Param "<< EnumToStringx(enum_type) << " cannot return a string array");}

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

    r22477 r22625  
    3131                case BasalforcingsPicoEnum:
    3232                        if(VerboseSolution())_printf0_(" call Pico Floating melting rate module\n");
    33                         PicoFloatingiceMeltingRatex(femmodel);
     33                        FloatingiceMeltingRatePicox(femmodel);
    3434                        break;
    3535                default:
    … …  
    5656/*}}}*/
    5757
    58 void PicoFloatingiceMeltingRatex(FemModel* femmodel){/*{{{*/
    59 
    60         int i,k,p,numvertices;
    61         int num_basins,maxbox,basinid,boxid;
    62         IssmDouble dist_max,area,earth_grav,rhoi,rhow,rho_star,nu,latentheat;
    63         IssmDouble c_p_ocean,lambda,a,b,c,alpha,Beta,gamma_T,overturning_coeff;
    64         IssmDouble g1,g2,s1,p_coeff,tavg,savg,overturning_avg;
    65         IssmDouble mean_toc,mean_soc,mean_overturning,pressure,T_star;
    66         IssmDouble q_coeff,potential_pressure_melting_point,T_pressure_melting,overturning;
    67         IssmDouble thickness,toc_farocean,soc_farocean,area_box1,Toc,Soc,basalmeltrate_shelf;
    68         int* nd=NULL;
    69         IssmDouble* dmax_basin=NULL;
    70         IssmDouble* t_farocean=NULL;
    71         IssmDouble* s_farocean=NULL;
    72         IssmDouble* distances=NULL;
    73         IssmDouble* boxareas=NULL;
    74         Element* element=NULL;
    75 
    76         femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
    77         femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
    78         dmax_basin=xNew<IssmDouble>(num_basins);
    79         femmodel->DistanceToFieldValue(MaskGroundediceLevelsetEnum,0.,DistanceToGroundinglineEnum);
    80         femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0.,DistanceToCalvingfrontEnum);
    81         // find maximum distance to grounding line
    82         dist_max=-1.;
    83         for(i=0;i<num_basins;i++){
    84                 dmax_basin[i]=-1;
    85         }
    86         for(i=0;i<femmodel->elements->Size();i++){
    87                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    88                 if(!element->IsFloating()) continue;
    89                 numvertices = element->GetNumberOfVertices();
    90                 distances=xNew<IssmDouble>(numvertices);
    91                 element->GetInputListOnVertices(&distances[0],DistanceToGroundinglineEnum);
    92                 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    93                 for(k=0; k<numvertices; k++){
    94                         if(fabs(distances[k])>dist_max){dist_max=fabs(distances[k]);}
    95                         if(fabs(distances[k])>dmax_basin[basinid]){dmax_basin[basinid]=fabs(distances[k]);}
    96                 }
    97         }
    98         //1 Define pico boxes
    99         nd=xNew<int>(num_basins);
    100         #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable.
    101         for(i=0; i<num_basins;i++){
    102                 nd[i] = round(sqrt(dmax_basin[i]/dist_max)*(maxbox-1)); //0-based
    103         }
    104         #else
    105         _error_("Don't know how to call round with AD on (present in other parts of the code as well");
    106         #endif
    107         for(int i=0;i<femmodel->elements->Size();i++){
    108                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    109                 element->PicoUpdateBoxid(nd);
    110         }
    111 
    112         //2 Get area of the boxes
    113         boxareas=xNew<IssmDouble>(num_basins*maxbox);
    114         for(i=0;i<num_basins*maxbox;i++){boxareas[i]=0.;}
    115         for(i=0;i<femmodel->elements->Size();i++){
    116                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    117                 if(!element->IsFloating()) continue;
    118                 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    119                 element->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
    120                 _error_("fix thi");
    121                 //boxareas[basinid*maxbox+boxid]+=element->GetArea();
    122         }
    123 
    124         //3 Compute variables in first box
    125         for(int p=0;p<femmodel->elements->Size();p++){
    126                 element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(p));
    127                 element->PicoUpdateFirstBox();
    128                 //Need to save the box1 average of temp, sal, and overt for input into next box.
    129         }
    130 
    131 
    132         //4 Compute variables of all other boxes
    133 
    134 //      for(int i=0;i<femmodel->elements->Size();i++){
    135 //              Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    136 //              element->PicoFloatingiceMeltingRate();
    137 //      }
    138         /*Cleanup and return */
    139         xDelete<int>(nd);
    140         xDelete<IssmDouble>(dmax_basin);
    141         xDelete<IssmDouble>(distances);
    142         xDelete<IssmDouble>(boxareas);
    143         xDelete<IssmDouble>(t_farocean);
    144    xDelete<IssmDouble>(s_farocean);
    145 }/*}}}*/

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

    r22471 r22625  
    77
    88#include "../../classes/classes.h"
     9#include "../FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.h"
    910
    1011/* local prototypes: */
    … …  
    1213void LinearFloatingiceMeltingRatex(FemModel* femmodel);
    1314void MismipFloatingiceMeltingRatex(FemModel* femmodel);
    14 void PicoFloatingiceMeltingRatex(FemModel* femmodel);
    1515
    1616#endif  /* _FloatingiceMeltingRatex_H*/

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

    r22606 r22625  
    211211                                iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.farocean_temperature");
    212212                                parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,&transparam[0],&transparam[M*(N-1)],interp,M,N));
     213                                xDelete<IssmDouble>(transparam);
    213214                                iomodel->FetchData(&transparam,&N,&M,"md.basalforcings.farocean_salinity");
    214215                                parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,&transparam[0],&transparam[M*(N-1)],interp,M,N));
     216                                xDelete<IssmDouble>(transparam);
    215217                        break;
    216218                default:

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

    r22498 r22625  
    6060#include "./Krigingx/Krigingx.h"
    6161#include "./FloatingiceMeltingRatex/FloatingiceMeltingRatex.h"
     62#include "./FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.h"
    6263#include "./Mergesolutionfromftogx/Mergesolutionfromftogx.h"
    6364#include "./MeshPartitionx/MeshPartitionx.h"

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

    r22620 r22625  
    6666        BasalforcingsMeltrateFactorEnum,
    6767        BasalforcingsNusseltEnum,
     68        BasalforcingsPicoAverageOverturningEnum,
     69        BasalforcingsPicoAverageSalinityEnum,
     70        BasalforcingsPicoAverageTemperatureEnum,
     71        BasalforcingsPicoBoxAreaEnum,
    6872        BasalforcingsPicoFarOceansalinityEnum,
     73        BasalforcingsPicoFarOceantemperatureEnum,
    6974        BasalforcingsPicoGammaTEnum,
     75        BasalforcingsPicoMaxboxcountEnum,
    7076        BasalforcingsPicoNumBasinsEnum,
    7177        BasalforcingsPicoOverturningCoeffEnum,
    … …  
    357363        BasalforcingsPicoBasinIdEnum,
    358364        BasalforcingsPicoBoxIdEnum,
    359         BasalforcingsPicoMaxboxcountEnum,
     365        BasalforcingsPicoSubShelfOceanOverturningEnum,
     366        BasalforcingsPicoSubShelfOceanSalinityEnum,
     367        BasalforcingsPicoSubShelfOceanTempEnum,
    360368        BaseEnum,
    361369        BedEnum,
    … …  
    624632        BasalCrevasseEnum,
    625633        BasalforcingsPicoEnum,
    626         BasalforcingsPicoFarOceantemperatureEnum,
    627634        BedSlopeSolutionEnum,
    628635        BoolExternalResultEnum,

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

    r22621 r22625  
    7474                case BasalforcingsMeltrateFactorEnum : return "BasalforcingsMeltrateFactor";
    7575                case BasalforcingsNusseltEnum : return "BasalforcingsNusselt";
     76                case BasalforcingsPicoAverageOverturningEnum : return "BasalforcingsPicoAverageOverturning";
     77                case BasalforcingsPicoAverageSalinityEnum : return "BasalforcingsPicoAverageSalinity";
     78                case BasalforcingsPicoAverageTemperatureEnum : return "BasalforcingsPicoAverageTemperature";
     79                case BasalforcingsPicoBoxAreaEnum : return "BasalforcingsPicoBoxArea";
    7680                case BasalforcingsPicoFarOceansalinityEnum : return "BasalforcingsPicoFarOceansalinity";
     81                case BasalforcingsPicoFarOceantemperatureEnum : return "BasalforcingsPicoFarOceantemperature";
    7782                case BasalforcingsPicoGammaTEnum : return "BasalforcingsPicoGammaT";
     83                case BasalforcingsPicoMaxboxcountEnum : return "BasalforcingsPicoMaxboxcount";
    7884                case BasalforcingsPicoNumBasinsEnum : return "BasalforcingsPicoNumBasins";
    7985                case BasalforcingsPicoOverturningCoeffEnum : return "BasalforcingsPicoOverturningCoeff";
    … …  
    363369                case BasalforcingsPicoBasinIdEnum : return "BasalforcingsPicoBasinId";
    364370                case BasalforcingsPicoBoxIdEnum : return "BasalforcingsPicoBoxId";
    365                 case BasalforcingsPicoMaxboxcountEnum : return "BasalforcingsPicoMaxboxcount";
     371                case BasalforcingsPicoSubShelfOceanOverturningEnum : return "BasalforcingsPicoSubShelfOceanOverturning";
     372                case BasalforcingsPicoSubShelfOceanSalinityEnum : return "BasalforcingsPicoSubShelfOceanSalinity";
     373                case BasalforcingsPicoSubShelfOceanTempEnum : return "BasalforcingsPicoSubShelfOceanTemp";
    366374                case BaseEnum : return "Base";
    367375                case BedEnum : return "Bed";
    … …  
    628636                case BasalCrevasseEnum : return "BasalCrevasse";
    629637                case BasalforcingsPicoEnum : return "BasalforcingsPico";
    630                 case BasalforcingsPicoFarOceantemperatureEnum : return "BasalforcingsPicoFarOceantemperature";
    631638                case BedSlopeSolutionEnum : return "BedSlopeSolution";
    632639                case BoolExternalResultEnum : return "BoolExternalResult";

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

    r22621 r22625  
    7474              else if (strcmp(name,"BasalforcingsMeltrateFactor")==0) return BasalforcingsMeltrateFactorEnum;
    7575              else if (strcmp(name,"BasalforcingsNusselt")==0) return BasalforcingsNusseltEnum;
     76              else if (strcmp(name,"BasalforcingsPicoAverageOverturning")==0) return BasalforcingsPicoAverageOverturningEnum;
     77              else if (strcmp(name,"BasalforcingsPicoAverageSalinity")==0) return BasalforcingsPicoAverageSalinityEnum;
     78              else if (strcmp(name,"BasalforcingsPicoAverageTemperature")==0) return BasalforcingsPicoAverageTemperatureEnum;
     79              else if (strcmp(name,"BasalforcingsPicoBoxArea")==0) return BasalforcingsPicoBoxAreaEnum;
    7680              else if (strcmp(name,"BasalforcingsPicoFarOceansalinity")==0) return BasalforcingsPicoFarOceansalinityEnum;
     81              else if (strcmp(name,"BasalforcingsPicoFarOceantemperature")==0) return BasalforcingsPicoFarOceantemperatureEnum;
    7782              else if (strcmp(name,"BasalforcingsPicoGammaT")==0) return BasalforcingsPicoGammaTEnum;
     83              else if (strcmp(name,"BasalforcingsPicoMaxboxcount")==0) return BasalforcingsPicoMaxboxcountEnum;
    7884              else if (strcmp(name,"BasalforcingsPicoNumBasins")==0) return BasalforcingsPicoNumBasinsEnum;
    7985              else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum;
    … …  
    131137              else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum;
    132138              else if (strcmp(name,"HydrologydcEplflipLock")==0) return HydrologydcEplflipLockEnum;
    133               else if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum;
     139         else stage=2;
     140   }
     141   if(stage==2){
     142              if (strcmp(name,"HydrologydcEplThickComp")==0) return HydrologydcEplThickCompEnum;
    134143              else if (strcmp(name,"HydrologydcIsefficientlayer")==0) return HydrologydcIsefficientlayerEnum;
    135144              else if (strcmp(name,"HydrologydcLeakageFactor")==0) return HydrologydcLeakageFactorEnum;
    … …  
    137146              else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
    138147              else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum;
    139          else stage=2;
    140    }
    141    if(stage==2){
    142               if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
     148              else if (strcmp(name,"HydrologydcRelTol")==0) return HydrologydcRelTolEnum;
    143149              else if (strcmp(name,"HydrologydcSedimentlimit")==0) return HydrologydcSedimentlimitEnum;
    144150              else if (strcmp(name,"HydrologydcSedimentlimitFlag")==0) return HydrologydcSedimentlimitFlagEnum;
    … …  
    254260              else if (strcmp(name,"SettingsSolverResidueThreshold")==0) return SettingsSolverResidueThresholdEnum;
    255261              else if (strcmp(name,"SettingsWaitonlock")==0) return SettingsWaitonlockEnum;
    256               else if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
     262         else stage=3;
     263   }
     264   if(stage==3){
     265              if (strcmp(name,"SmbAdThresh")==0) return SmbAdThreshEnum;
    257266              else if (strcmp(name,"SmbAIce")==0) return SmbAIceEnum;
    258267              else if (strcmp(name,"SmbAIdx")==0) return SmbAIdxEnum;
    … …  
    260269              else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum;
    261270              else if (strcmp(name,"SmbDelta18o")==0) return SmbDelta18oEnum;
    262          else stage=3;
    263    }
    264    if(stage==3){
    265               if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
     271              else if (strcmp(name,"SmbDelta18oSurface")==0) return SmbDelta18oSurfaceEnum;
    266272              else if (strcmp(name,"SmbDenIdx")==0) return SmbDenIdxEnum;
    267273              else if (strcmp(name,"SmbDt")==0) return SmbDtEnum;
    … …  
    369375              else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum;
    370376              else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum;
    371               else if (strcmp(name,"BasalforcingsPicoMaxboxcount")==0) return BasalforcingsPicoMaxboxcountEnum;
     377              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum;
     378              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum;
     379              else if (strcmp(name,"BasalforcingsPicoSubShelfOceanTemp")==0) return BasalforcingsPicoSubShelfOceanTempEnum;
    372380              else if (strcmp(name,"Base")==0) return BaseEnum;
    373381              else if (strcmp(name,"Bed")==0) return BedEnum;
    … …  
    375383              else if (strcmp(name,"BedSlopeY")==0) return BedSlopeYEnum;
    376384              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    377               else if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"CalvingHabFraction")==0) return CalvingHabFractionEnum;
    378389              else if (strcmp(name,"CalvinglevermannCoeff")==0) return CalvinglevermannCoeffEnum;
    379390              else if (strcmp(name,"CalvinglevermannMeltingrate")==0) return CalvinglevermannMeltingrateEnum;
    … …  
    383394              else if (strcmp(name,"CalvingrateyAverage")==0) return CalvingrateyAverageEnum;
    384395              else if (strcmp(name,"Calvingratey")==0) return CalvingrateyEnum;
    385          else stage=4;
    386    }
    387    if(stage==4){
    388               if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
     396              else if (strcmp(name,"CalvingStressThresholdFloatingice")==0) return CalvingStressThresholdFloatingiceEnum;
    389397              else if (strcmp(name,"CalvingStressThresholdGroundedice")==0) return CalvingStressThresholdGroundediceEnum;
    390398              else if (strcmp(name,"Converged")==0) return ConvergedEnum;
    … …  
    498506              else if (strcmp(name,"SmbAini")==0) return SmbAiniEnum;
    499507              else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum;
    500               else if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"SmbBMax")==0) return SmbBMaxEnum;
    501512              else if (strcmp(name,"SmbBMin")==0) return SmbBMinEnum;
    502513              else if (strcmp(name,"SmbBNeg")==0) return SmbBNegEnum;
    … …  
    506517              else if (strcmp(name,"SmbDini")==0) return SmbDiniEnum;
    507518              else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum;
    508          else stage=5;
    509    }
    510    if(stage==5){
    511               if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
     519              else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum;
    512520              else if (strcmp(name,"SmbDz")==0) return SmbDzEnum;
    513521              else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
    … …  
    621629              else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum;
    622630              else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum;
    623               else if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Adjointp")==0) return AdjointpEnum;
    624635              else if (strcmp(name,"AggressiveMigration")==0) return AggressiveMigrationEnum;
    625636              else if (strcmp(name,"AmrBamg")==0) return AmrBamgEnum;
    … …  
    629640              else if (strcmp(name,"AutodiffJacobian")==0) return AutodiffJacobianEnum;
    630641              else if (strcmp(name,"Balancethickness2Analysis")==0) return Balancethickness2AnalysisEnum;
    631          else stage=6;
    632    }
    633    if(stage==6){
    634               if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
     642              else if (strcmp(name,"Balancethickness2Solution")==0) return Balancethickness2SolutionEnum;
    635643              else if (strcmp(name,"BalancethicknessAnalysis")==0) return BalancethicknessAnalysisEnum;
    636644              else if (strcmp(name,"BalancethicknessApparentMassbalance")==0) return BalancethicknessApparentMassbalanceEnum;
    … …  
    643651              else if (strcmp(name,"BasalCrevasse")==0) return BasalCrevasseEnum;
    644652              else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum;
    645               else if (strcmp(name,"BasalforcingsPicoFarOceantemperature")==0) return BasalforcingsPicoFarOceantemperatureEnum;
    646653              else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
    647654              else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum;
    … …  
    745752              else if (strcmp(name,"GroundedAreaScaled")==0) return GroundedAreaScaledEnum;
    746753              else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum;
    747               else if (strcmp(name,"Gset")==0) return GsetEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"Gset")==0) return GsetEnum;
    748758              else if (strcmp(name,"Gsl")==0) return GslEnum;
    749759              else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
    … …  
    752762              else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum;
    753763              else if (strcmp(name,"HydrologyDCEfficientAnalysis")==0) return HydrologyDCEfficientAnalysisEnum;
    754          else stage=7;
    755    }
    756    if(stage==7){
    757               if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
     764              else if (strcmp(name,"Hydrologydc")==0) return HydrologydcEnum;
    758765              else if (strcmp(name,"HydrologydcEplColapseThickness")==0) return HydrologydcEplColapseThicknessEnum;
    759766              else if (strcmp(name,"HydrologydcEplCompressibility")==0) return HydrologydcEplCompressibilityEnum;
    … …  
    868875              else if (strcmp(name,"MinVx")==0) return MinVxEnum;
    869876              else if (strcmp(name,"MinVy")==0) return MinVyEnum;
    870               else if (strcmp(name,"MinVz")==0) return MinVzEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"MinVz")==0) return MinVzEnum;
    871881              else if (strcmp(name,"MismipFloatingMeltRate")==0) return MismipFloatingMeltRateEnum;
    872882              else if (strcmp(name,"Moulin")==0) return MoulinEnum;
    … …  
    875885              else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum;
    876886              else if (strcmp(name,"Mumps")==0) return MumpsEnum;
    877          else stage=8;
    878    }
    879    if(stage==8){
    880               if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
     887              else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    881888              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    882889              else if (strcmp(name,"Nodalvalue")==0) return NodalvalueEnum;
    … …  
    991998              else if (strcmp(name,"Outputdefinition100")==0) return Outputdefinition100Enum;
    992999              else if (strcmp(name,"P0Array")==0) return P0ArrayEnum;
    993               else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
    9941004              else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum;
    9951005              else if (strcmp(name,"P1DG")==0) return P1DGEnum;
    … …  
    9981008              else if (strcmp(name,"P1xP2")==0) return P1xP2Enum;
    9991009              else if (strcmp(name,"P1xP3")==0) return P1xP3Enum;
    1000          else stage=9;
    1001    }
    1002    if(stage==9){
    1003               if (strcmp(name,"P1xP4")==0) return P1xP4Enum;
     1010              else if (strcmp(name,"P1xP4")==0) return P1xP4Enum;
    10041011              else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
    10051012              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
    … …  
    11141121              else if (strcmp(name,"Water")==0) return WaterEnum;
    11151122              else if (strcmp(name,"WorldComm")==0) return WorldCommEnum;
    1116               else if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
     1123         else stage=10;
     1124   }
     1125   if(stage==10){
     1126              if (strcmp(name,"XTaylorHood")==0) return XTaylorHoodEnum;
    11171127              else if (strcmp(name,"XY")==0) return XYEnum;
    11181128              else if (strcmp(name,"XYZ")==0) return XYZEnum;
    … …  
    11211131              else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum;
    11221132              else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum;
    1123          else stage=10;
    1124    }
    1125    if(stage==10){
    1126               if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
     1133              else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum;
    11271134              else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum;
    11281135              else if (strcmp(name,"MeshZ")==0) return MeshZEnum;

  • issm/trunk-jpl/src/m/classes/basalforcingspico.m

    r22419 r22625  
    1111                overturning_coeff         = 0.;
    1212                gamma_T                   = 0.;
    13                 box0temperature           = NaN;
    14                 box0salinity              = NaN;
     13                farocean_temperature      = NaN;
     14                farocean_salinity         = NaN;
    1515                geothermalflux            = NaN;
     16                groundedice_melting_rate  = NaN;
    1617        end
    1718        methods
    … …  
    1920                        self.basin_id=project3d(md,'vector',self.basin_id,'type','element','layer',1);
    2021                        self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','element','layer',1); %bedrock only gets geothermal flux
     22                        self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1);
    2123                end % }}}
    2224                function self = basalforcingspico(varargin) % {{{
    … …  
    4648                                disp('      no turbulent temperature exchange velocity set, setting value to 2e-5');
    4749                        end
     50                        if isnan(self.groundedice_melting_rate),
     51                                self.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1);
     52                                disp('      no basalforcings.groundedice_melting_rate specified: values set as zero');
     53                        end
     54
    4855                end % }}}
    4956                function self = setdefaultparameters(self) % {{{
    … …  
    5764
    5865                                md = checkfield(md,'fieldname','basalforcings.num_basins','numel',1,'NaN',1,'Inf',1,'>',0);
    59                                 md = checkfield(md,'fieldname','basalforcings.basin_id','Inf',1,'>',0,'<=',md.basalforcings.num_basins,'size',[md.mesh.numberofelements 1]);
     66                                md = checkfield(md,'fieldname','basalforcings.basin_id','Inf',1,'>=',0,'<=',md.basalforcings.num_basins,'size',[md.mesh.numberofelements 1]);
    6067                                md = checkfield(md,'fieldname','basalforcings.maxboxcount','numel',1,'NaN',1,'Inf',1,'>',0);
    6168                                md = checkfield(md,'fieldname','basalforcings.overturning_coeff','numel',1,'NaN',1,'Inf',1,'>',0);
    6269                                md = checkfield(md,'fieldname','basalforcings.gamma_T','numel',1,'NaN',1,'Inf',1,'>',0);
    63                                 md = checkfield(md,'fieldname','basalforcings.box0temperature','NaN',1,'Inf',1,'>',0);
    64                                 md = checkfield(md,'fieldname','basalforcings.box0salinity','NaN',1,'Inf',1,'>',0);
     70                                md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'>',0);
     71                                md = checkfield(md,'fieldname','basalforcings.farocean_salinity','NaN',1,'Inf',1,'>',0);
    6572                                md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'>=',0,'timeseries',1);
     73                                md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1);
    6674
    6775                end % }}}
    … …  
    7381                        fielddisplay(self,'overturning_coeff','Overturning strength [m^3/s]');
    7482                        fielddisplay(self,'gamma_T','Turbulent temperature exchange velocity [m/s]');
    75                         fielddisplay(self,'box0temperature','depth averaged ocean temperature in front of the ice shelf for basin i [K]');
    76                         fielddisplay(self,'box0salinity','depth averaged ocean salinity in front of the ice shelf for basin i [psu]');
     83                        fielddisplay(self,'farocean_temperature','depth averaged ocean temperature in front of the ice shelf for basin i [K]');
     84                        fielddisplay(self,'farocean_salinity','depth averaged ocean salinity in front of the ice shelf for basin i [psu]');
    7785                        fielddisplay(self,'geothermalflux','geothermal heat flux [W/m^2]');
     86                        fielddisplay(self,'groundedice_melting_rate','basal melting rate (positive if melting) [m/yr]');
    7887
    7988                end % }}}
    … …  
    8796                        WriteData(fid,prefix,'object',self,'fieldname','overturning_coeff','format','Double');
    8897                        WriteData(fid,prefix,'object',self,'fieldname','gamma_T','format','Double');
    89                         WriteData(fid,prefix,'object',self,'fieldname','box0temperature','format','DoubleMat','name','md.basalforcings.box0temperature','timeserieslength',md.basalforcings.num_basins+1,'yts',md.constants.yts);
    90                         WriteData(fid,prefix,'object',self,'fieldname','box0salinity','format','DoubleMat','name','md.basalforcings.box0salinity','timeserieslength',md.basalforcings.num_basins+1,'yts',md.constants.yts);
    91                         WriteData(fid,prefix,'object',self,'fieldname','basin_id','format','DoubleMat','name','md.basalforcings.basin_id','mattype',2);
     98                        WriteData(fid,prefix,'object',self,'fieldname','farocean_temperature','format','DoubleMat','name','md.basalforcings.farocean_temperature','timeserieslength',md.basalforcings.num_basins+1,'yts',md.constants.yts);
     99                        WriteData(fid,prefix,'object',self,'fieldname','farocean_salinity','format','DoubleMat','name','md.basalforcings.farocean_salinity','timeserieslength',md.basalforcings.num_basins+1,'yts',md.constants.yts);
     100                        %WriteData(fid,prefix,'object',self,'fieldname','basin_id','format','DoubleMat','name','md.basalforcings.basin_id','mattype',2);
     101                        WriteData(fid,prefix,'object',self,'fieldname','basin_id','data',self.basin_id-1,'name','md.basalforcings.basin_id','format','IntMat','mattype',2);   %Change to 0-indexing
    92102                        WriteData(fid,prefix,'object',self,'fieldname','geothermalflux','format','DoubleMat','name','md.basalforcings.geothermalflux','mattype',1,'timeserieslength',md.mesh.numberofelements+1,'yts',md.constants.yts);
     103                        WriteData(fid,prefix,'object',self,'fieldname','groundedice_melting_rate','format','DoubleMat','mattype',1,'scale',1./yts,'timeserieslength',md.mesh.numberofvertices+1,'yts',md.constants.yts);
    93104
    94105                end % }}}
Note: See TracChangeset for help on using the changeset viewer.