Changeset 24010


Ignore:
Timestamp:
06/11/19 11:07:32 (6 years ago)
Author:
tpelle
Message:

NEW: ISMIP6 and PICO basal melting parameterizations working in HO, PICOP coming soon

Location:
issm/trunk-jpl/src
Files:
12 edited

Legend:

Unmodified
Added
Removed

  • TabularUnified issm/trunk-jpl/src/c/analyses/GLheightadvectionAnalysis.cpp

    r23966 r24010  
    1717
    1818        /*First fetch data: */
     19        if(iomodel->domaintype==Domain3DEnum) iomodel->FetchData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    1920        ::CreateNodes(nodes,iomodel,GLheightadvectionAnalysisEnum,P1Enum);
     21        iomodel->DeleteData(2,"md.mesh.vertexonbase","md.mesh.vertexonsurface");
    2022}/*}}}*/
    2123int  GLheightadvectionAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
    … …  
    5658ElementMatrix* GLheightadvectionAnalysis::CreateKMatrix(Element* element){/*{{{*/
    5759
     60        /* Spawn basal element */
     61        if(!element->IsOnBase()) return NULL;
     62        Element* basalelement = element->SpawnBasalElement();
     63
    5864        /* Check if ice in element */
    59         if(!element->IsIceInElement()) return NULL;
    60         //if(!element->IsFloating()) return NULL;
     65        if(!basalelement->IsIceInElement()) return NULL;
    6166
    6267        /*Intermediaries */
    … …  
    6570        IssmDouble Jdet,D_scalar,onboundary;
    6671        IssmDouble vel,vx,vy;
    67         IssmDouble* xyz_list = NULL;
     72        IssmDouble* xyz_list      = NULL;
     73        Input* vx_input           = NULL;
     74        Input* vy_input           = NULL;
     75        Input* bc_input           = NULL;
    6876
    6977        /*Get problem dimension*/
    70         element->FindParam(&domaintype,DomainTypeEnum);
     78        basalelement->FindParam(&domaintype,DomainTypeEnum);
    7179        switch(domaintype){
    7280                case Domain2DhorizontalEnum: dim = 2; break;
    … …  
    7684
    7785        /*Fetch number of nodes and dof for this finite element*/
    78         int numnodes = element->GetNumberOfNodes();
     86        int numnodes = basalelement->GetNumberOfNodes();
    7987
    8088        /*Initialize Element vector and other vectors*/
    81         ElementMatrix* Ke     = element->NewElementMatrix();
     89        ElementMatrix* Ke     = basalelement->NewElementMatrix();
    8290        IssmDouble*    basis  = xNew<IssmDouble>(numnodes);
    8391        IssmDouble*    dbasis = xNew<IssmDouble>(dim*numnodes);
    … …  
    8694
    8795        /*Retrieve all inputs and parameters*/
    88         element->GetVerticesCoordinates(&xyz_list);
    89         //Input* vx_input=element->GetInput(BaseSlopeXEnum); _assert_(vx_input);
    90         //Input* vy_input=element->GetInput(BaseSlopeYEnum); _assert_(vy_input);
    91         Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input);
    92         Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input);
    93         Input* bc_input=element->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input);
    94 
    95         IssmDouble h = element->CharacteristicLength();
     96        basalelement->GetVerticesCoordinates(&xyz_list);
     97        switch(domaintype){
     98                case Domain2DhorizontalEnum:
     99                        vx_input=basalelement->GetInput(VxEnum); _assert_(vx_input);
     100                        vy_input=basalelement->GetInput(VyEnum); _assert_(vy_input);
     101                   bc_input=basalelement->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input);
     102                break;
     103                case Domain3DEnum:
     104                        vx_input=basalelement->GetInput(VxAverageEnum); _assert_(vx_input);
     105                        vy_input=basalelement->GetInput(VyAverageEnum); _assert_(vy_input);
     106                        bc_input=basalelement->GetInput(MeshVertexonboundaryEnum); _assert_(bc_input);
     107                break;
     108                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
     109        }
     110        IssmDouble h = basalelement->CharacteristicLength();
    96111
    97112        /* Start  looping on the number of gaussian points: */
    98         Gauss* gauss=element->NewGauss(4);
     113        Gauss* gauss=basalelement->NewGauss(4);
    99114        for(int ig=gauss->begin();ig<gauss->end();ig++){
    100115                gauss->GaussPoint(ig);
    101116
    102                 element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    103                 element->NodalFunctions(basis,gauss);
    104                 element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
    105                 GetBprime(Bprime,element,dim,xyz_list,gauss);
     117                basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
     118                basalelement->NodalFunctions(basis,gauss);
     119                basalelement->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     120                GetBprime(Bprime,basalelement,dim,xyz_list,gauss);
    106121
    107122                /*Get velocity*/
    … …  
    168183        xDelete<IssmDouble>(Bprime);
    169184        delete gauss;
     185        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    170186        return Ke;
    171187
    172188}/*}}}*/
    173189ElementVector* GLheightadvectionAnalysis::CreatePVector(Element* element){/*{{{*/
    174 
    175190        return NULL;
    176191}/*}}}*/
    … …  
    203218        for(int i=0;i<femmodel->elements->Size();i++){
    204219                Element    *element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    205 
    206220                int         numnodes  = element->GetNumberOfNodes();
    207221                IssmDouble *mask      = xNew<IssmDouble>(numnodes);

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

    r23999 r24010  
    21572157}
    21582158/*}}}*/
     2159void       Element::Ismip6FloatingiceMeltingRate(){/*{{{*/
     2160
     2161        if(!this->IsIceInElement() || !this->IsFloating()) return;
     2162
     2163        int         numvertices = this->GetNumberOfVertices();
     2164        int         basinid,num_basins,M,N;
     2165        IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval;
     2166        IssmDouble  basalmeltrate[numvertices];
     2167        bool        islocal;
     2168        IssmDouble* delta_t = NULL;
     2169        IssmDouble* mean_tf = NULL;
     2170        IssmDouble* depths  = NULL;
     2171
     2172        /*Get variables*/
     2173        IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum);
     2174        IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum);
     2175        IssmDouble lf   = this->FindParam(MaterialsLatentheatEnum);
     2176        IssmDouble cp   = this->FindParam(MaterialsMixedLayerCapacityEnum);
     2177
     2178        /*Hard code sea water density to be consistent with ISMIP6 documentation*/
     2179        rhow = 1028.;
     2180
     2181        /* Get parameters and inputs */
     2182        this->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
     2183        this->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum);
     2184        this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum);
     2185        this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum);    _assert_(M==num_basins);
     2186        this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum);
     2187        if(!islocal) {
     2188                this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins);
     2189        }
     2190        Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum);            _assert_(tf_input);
     2191        delta_t_basin = delta_t[basinid];
     2192        if(!islocal) mean_tf_basin = mean_tf[basinid];
     2193
     2194        /*Compute melt rate for Local and Nonlocal parameterizations*/
     2195        Gauss* gauss=this->NewGauss();
     2196        for(int i=0;i<numvertices;i++){
     2197                gauss->GaussVertex(i);
     2198                tf_input->GetInputValue(&tf,gauss);
     2199                if(!islocal) {
     2200                        absval = mean_tf_basin+delta_t_basin;
     2201                        if (absval<0) {absval = -1.*absval;}
     2202                        basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;}
     2203                else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);}
     2204        }
     2205
     2206        /*Return basal melt rate*/
     2207        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum);
     2208
     2209        /*Cleanup and return*/
     2210        delete gauss;
     2211        xDelete<IssmDouble>(delta_t);
     2212        xDelete<IssmDouble>(mean_tf);
     2213        xDelete<IssmDouble>(depths);
     2214
     2215}/*}}}*/
    21592216bool       Element::IsWaterInElement(){/*{{{*/
    21602217        return (this->inputs->Max(MaskOceanLevelsetEnum)>0.);
    … …  
    25202577}
    25212578/*}}}*/
     2579void       Element::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/
     2580
     2581        if(!this->IsIceInElement() || !this->IsFloating()) return;
     2582
     2583        int        basin_id;
     2584        IssmDouble dist_gl,dist_cf;
     2585
     2586        inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
     2587        IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id])+1.;
     2588
     2589        Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
     2590        Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum);  _assert_(dist_cf_input);
     2591
     2592        /*Get dist_gl and dist_cf at center of element*/
     2593        Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
     2594        dist_gl_input->GetInputValue(&dist_gl,gauss);
     2595        dist_cf_input->GetInputValue(&dist_cf,gauss);
     2596        delete gauss;
     2597
     2598        /*Ensure values are positive for floating ice*/
     2599        dist_gl = fabs(dist_gl);
     2600        dist_cf = fabs(dist_cf);
     2601
     2602        /*Compute relative distance to grounding line*/
     2603        IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf);
     2604
     2605        /*Assign box numbers based on rel_dist_gl*/
     2606        int boxid = -1;
     2607        for(IssmDouble i=0.;i<boxid_max;i++){
     2608                IssmDouble lowbound  = 1. -sqrt((boxid_max-i   )/boxid_max);
     2609                IssmDouble highbound = 1. -sqrt((boxid_max-i-1.)/boxid_max);
     2610                if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){
     2611                        boxid=reCast<int>(i);
     2612                        break;
     2613                }
     2614        }
     2615        if(boxid==-1) _error_("No boxid found for element " << this->Sid() << "!");
     2616
     2617        this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));
     2618       
     2619}/*}}}*/
     2620void       Element::PicoUpdateBox(int loopboxid){/*{{{*/
     2621
     2622        if(!this->IsIceInElement() || !this->IsFloating()) return;
     2623
     2624        int boxid;
     2625        this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
     2626        if(loopboxid!=boxid) return;
     2627
     2628        int        NUMVERTICES = this->GetNumberOfVertices();
     2629        int        basinid, maxbox, num_basins, numnodes, M;
     2630        IssmDouble gamma_T, overturning_coeff, thickness;
     2631        IssmDouble pressure, T_star,p_coeff, q_coeff;
     2632        bool       isplume;
     2633
     2634        /*Get variables*/
     2635        IssmDouble rhoi       = this->FindParam(MaterialsRhoIceEnum);
     2636        IssmDouble rhow       = this->FindParam(MaterialsRhoSeawaterEnum);
     2637        IssmDouble earth_grav = this->FindParam(ConstantsGEnum);
     2638        IssmDouble rho_star   = 1033.;             // kg/m^3
     2639        IssmDouble nu         = rhoi/rhow;
     2640        IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum);
     2641        IssmDouble c_p_ocean  = this->FindParam(MaterialsMixedLayerCapacityEnum);
     2642        IssmDouble lambda     = latentheat/c_p_ocean;
     2643        IssmDouble a          = -0.0572;          // K/psu
     2644        IssmDouble b          = 0.0788 + this->FindParam(MaterialsMeltingpointEnum);  //K
     2645        IssmDouble c          = 7.77e-4;
     2646        IssmDouble alpha      = 7.5e-5;           // 1/K
     2647        IssmDouble Beta       = 7.7e-4;           // K
     2648
     2649        /* Get non-box-specific parameters and inputs */
     2650        this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);
     2651        this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
     2652        this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
     2653        this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
     2654        this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
     2655        this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum);
     2656        Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
     2657        _assert_(basinid<=num_basins);
     2658
     2659        IssmDouble* boxareas = NULL;
     2660        this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
     2661        _assert_(M==num_basins*maxbox);
     2662
     2663        IssmDouble area_boxi        = boxareas[basinid*maxbox+boxid];
     2664        IssmDouble g1               = area_boxi*gamma_T;
     2665
     2666        IssmDouble basalmeltrates_shelf[NUMVERTICES];
     2667        IssmDouble potential_pressure_melting_point[NUMVERTICES];
     2668        IssmDouble Tocs[NUMVERTICES];
     2669        IssmDouble Socs[NUMVERTICES];
     2670
     2671        /* First box calculations */
     2672        if(boxid==0){
     2673                /* Get box1 parameters and inputs */
     2674                IssmDouble time, toc_farocean, soc_farocean;
     2675                this->parameters->FindParam(&time,TimeEnum);
     2676                this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
     2677                this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
     2678                IssmDouble s1 = soc_farocean/(nu*lambda);
     2679                IssmDouble overturnings[NUMVERTICES];
     2680
     2681                /* Start looping on the number of verticies and calculate ocean vars */
     2682                Gauss* gauss=this->NewGauss();
     2683                for(int i=0;i<NUMVERTICES;i++){
     2684                        gauss->GaussVertex(i);
     2685                        thickness_input->GetInputValue(&thickness,gauss);
     2686                        pressure = (rhoi*earth_grav*1e-4)*thickness;
     2687                        T_star   = a*soc_farocean+b-c*pressure-toc_farocean;
     2688                        p_coeff  = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
     2689                        q_coeff  = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
     2690
     2691                        /* To avoid negatives under the square root */
     2692                        if((0.25*pow(p_coeff,2)-q_coeff)<0) q_coeff = 0.25*p_coeff*p_coeff;
     2693
     2694                        Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff));
     2695                        Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]);
     2696                        potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
     2697                        if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
     2698                        overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i]));
     2699                }
     2700
     2701                if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
     2702                this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
     2703                this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
     2704                this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum);
     2705
     2706                /*Cleanup and return*/
     2707                delete gauss;
     2708        }
     2709
     2710        /* Subsequent box calculations */
     2711        else {
     2712                /* Get subsequent box parameters and inputs */
     2713                IssmDouble* toc_weighted_avg         = NULL;
     2714                IssmDouble* soc_weighted_avg         = NULL;
     2715                IssmDouble* overturning_weighted_avg = NULL;
     2716                this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum);
     2717                this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);
     2718                this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);
     2719                IssmDouble mean_toc                  = toc_weighted_avg[basinid];
     2720                IssmDouble mean_soc                  = soc_weighted_avg[basinid];
     2721                IssmDouble mean_overturning          = overturning_weighted_avg[basinid];
     2722                IssmDouble g2                        = g1/(nu*lambda);
     2723
     2724                /* Start looping on the number of verticies and calculate ocean vars */
     2725                Gauss* gauss=this->NewGauss();
     2726                for(int i=0;i<NUMVERTICES;i++){
     2727                        gauss->GaussVertex(i);
     2728                        thickness_input->GetInputValue(&thickness,gauss);
     2729                        pressure = (rhoi*earth_grav*1.e-4)*thickness;
     2730                        T_star   = a*mean_soc+b-c*pressure-mean_toc;
     2731                        Tocs[i]  = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc));
     2732                        Socs[i]  = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda));
     2733                        potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
     2734                        if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
     2735                }
     2736
     2737                if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
     2738                this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
     2739                this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
     2740
     2741                /*Cleanup and return*/
     2742                xDelete<IssmDouble>(toc_weighted_avg);
     2743                xDelete<IssmDouble>(soc_weighted_avg);
     2744                xDelete<IssmDouble>(overturning_weighted_avg);
     2745                delete gauss;
     2746        }
     2747
     2748        /*Cleanup and return*/
     2749        xDelete<IssmDouble>(boxareas);
     2750       
     2751}/*}}}*/
     2752void       Element::PicoComputeBasalMelt(){/*{{{*/
     2753
     2754        if(!this->IsIceInElement() || !this->IsFloating()) return;
     2755
     2756        int        NUMVERTICES = this->GetNumberOfVertices();
     2757        IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0;
     2758        IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat;
     2759        IssmDouble alpha, zgl, Toc, Soc, z_base, yts, slopex, slopey;
     2760
     2761        /*Get variables*/
     2762        E0    = 3.6e-2;        //Entrainment coefficient
     2763        Cd    = 2.5e-3;        //Drag coefficient
     2764        CdT   = 1.1e-3;        //Turbulent heat exchange coefficient
     2765        YT    = CdT/sqrt(Cd);  //Heat exchange coefficient
     2766        lam1  = -5.73e-2;      //Freezing point-salinity coefficient (degrees C)
     2767        lam2  = 8.32e-2;       //Freezing point offset (degrees C)
     2768        lam3  = 7.61e-4;       //Freezing point-depth coefficient (K m-1)
     2769        M0    = 10.;           //Melt-rate parameter (m yr-1 C-2)
     2770        CdTS0 = 6e-4;          //Heat exchange parameter
     2771        y1    = 0.545;         //Heat exchange parameter
     2772        y2    = 3.5e-5;        //Heat exchange parameter
     2773        x0    = 0.56;          //Dimentionless scaling factor
     2774
     2775        /*Define arrays*/
     2776        IssmDouble basalmeltrates_shelf[NUMVERTICES];  //Basal melt-rate
     2777
     2778        /*Polynomial coefficients*/
     2779        IssmDouble p[12];
     2780        p[0]  =  0.1371330075095435;
     2781        p[1]  =  5.527656234709359E1;
     2782        p[2]  = -8.951812433987858E2;
     2783        p[3]  =  8.927093637594877E3;
     2784        p[4]  = -5.563863123811898E4;
     2785        p[5]  =  2.218596970948727E5;
     2786        p[6]  = -5.820015295669482E5;
     2787        p[7]  =  1.015475347943186E6;
     2788        p[8]  = -1.166290429178556E6;
     2789        p[9]  =  8.466870335320488E5;
     2790        p[10] = -3.520598035764990E5;
     2791        p[11] =  6.387953795485420E4;
     2792
     2793        /*Get inputs*/
     2794        Input* zgl_input            = this->GetInput(GroundinglineHeightEnum);                     _assert_(zgl_input);
     2795        Input* toc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum);      _assert_(toc_input);
     2796        Input* soc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum);  _assert_(soc_input);
     2797        Input* base_input           = this->GetInput(BaseEnum);                                    _assert_(base_input);
     2798        Input* baseslopex_input     = this->GetInput(BaseSlopeXEnum);                              _assert_(baseslopex_input);
     2799        Input* baseslopey_input     = this->GetInput(BaseSlopeYEnum);                              _assert_(baseslopey_input);
     2800        this->FindParam(&yts, ConstantsYtsEnum);
     2801
     2802        /*Loop over the number of vertices in this element*/
     2803        Gauss* gauss=this->NewGauss();
     2804        for(int i=0;i<NUMVERTICES;i++){
     2805                gauss->GaussVertex(i);
     2806
     2807                /*Get inputs*/
     2808                zgl_input->GetInputValue(&zgl,gauss);
     2809                toc_input->GetInputValue(&Toc,gauss); //(K)
     2810                soc_input->GetInputValue(&Soc,gauss);
     2811                base_input->GetInputValue(&z_base,gauss);
     2812                baseslopex_input->GetInputValue(&slopex,gauss);
     2813                baseslopey_input->GetInputValue(&slopey,gauss);
     2814
     2815                /*Compute ice shelf base slope (radians)*/
     2816                alpha = atan(sqrt(slopex*slopex + slopey*slopey));
     2817                if(alpha>=M_PI) alpha = M_PI - 0.001;               //ensure sin(alpha) > 0 for meltrate calculations
     2818
     2819                /*Make necessary conversions*/
     2820                Toc = Toc-273.15;
     2821                if(zgl>z_base) zgl=z_base;
     2822
     2823                /*Low bound for Toc to ensure X_hat is between 0 and 1*/
     2824                if(Toc<lam1*Soc+lam2) Toc=lam1*Soc+lam2;
     2825
     2826                /*Compute parameters needed for melt-rate calculation*/
     2827                Tf_gl = lam1*Soc+lam2+lam3*zgl;                                              //Characteristic freezing point
     2828                YTS = YT*(y1+y2*(((Toc-Tf_gl)*E0*sin(alpha))/(lam3*(CdTS0+E0*sin(alpha))))); //Effective heat exchange coefficient
     2829                CdTS = sqrt(Cd)*YTS;                                                         //Heat exchange coefficient
     2830                G1 = sqrt(sin(alpha)/(Cd+E0*sin(alpha)));                                    //Geometric factor
     2831                G2 = sqrt(CdTS/(CdTS+E0*sin(alpha)));                                        //Geometric factor
     2832                G3 = (E0*sin(alpha))/(CdTS+E0*sin(alpha));                                   //Geometric factor
     2833                g_alpha = G1*G2*G3;                                                          //Melt scaling factor
     2834                M = M0*g_alpha*pow((Toc-Tf_gl),2);                                           //Melt-rate scale
     2835                l = ((Toc-Tf_gl)*(x0*CdTS+E0*sin(alpha)))/(lam3*x0*(CdTS+E0*sin(alpha)));    //Length scale
     2836                X_hat = (z_base-zgl)/l;                                                      //Dimentionless coordinate system
     2837
     2838                /*Compute polynomial fit*/
     2839                M_hat = 0.;                                                                  //Reset summation variable for each node
     2840                for(int ii=0;ii<12;ii++) {
     2841                        M_hat += p[ii]*pow(X_hat,ii);                                               //Polynomial fit
     2842                }
     2843
     2844                /*Compute melt-rate*/
     2845                basalmeltrates_shelf[i] = (M*M_hat)/yts;                                     //Basal melt-rate (m/s)
     2846        }
     2847
     2848        /*Save computed melt-rate*/
     2849        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
     2850
     2851        /*Cleanup and return*/
     2852        delete gauss;
     2853               
     2854}/*}}}*/
    25222855void       Element::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac){/*{{{*/
    25232856

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

    r23993 r24010  
    134134                bool               IsIceInElement();
    135135                bool               IsLandInElement();
     136                void               Ismip6FloatingiceMeltingRate();
    136137                bool               IsWaterInElement();
    137138                void               LinearFloatingiceMeltingRate();
    … …  
    145146                ElementMatrix*     NewElementMatrixCoupling(int number_nodes,int approximation_enum=NoneApproximationEnum);
    146147                ElementVector*     NewElementVector(int approximation_enum=NoneApproximationEnum);
     148      void               PicoUpdateBoxid(int* pmax_boxid_basin);
     149                void               PicoUpdateBox(int loopboxid);
     150                void               PicoComputeBasalMelt();
    147151                void               PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm,bool ismungsm,bool issetpddfac);
    148152                void               PositiveDegreeDaySicopolis(bool isfirnwarming);
    … …  
    257261                virtual bool       IsFaceOnBoundary(void)=0;
    258262                virtual bool       IsIcefront(void)=0;
    259                 virtual void       Ismip6FloatingiceMeltingRate(void)=0;
    260263                virtual bool       IsNodeOnShelfFromFlags(IssmDouble* flags)=0;
    261264                virtual bool       IsOnBase()=0;
    … …  
    299302                virtual int        NumberofNodesPressure(void)=0;
    300303                virtual int        NumberofNodesVelocity(void)=0;
    301                 virtual void       PicoUpdateBoxid(int* pmax_boxid_basin)=0;
    302                 virtual void       PicoUpdateBox(int loopboxid)=0;
    303                 virtual void       PicoComputeBasalMelt(void)=0;
    304304                virtual void       PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding)=0;
    305305                virtual int        PressureInterpolation()=0;

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

    r23801 r24010  
    9595                IssmDouble     IceVolume(bool scaled);
    9696                IssmDouble     IceVolumeAboveFloatation(bool scaled);
    97                 void                            Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");};
    9897                void           InputDepthAverageAtBase(int enum_type,int average_enum_type);
    9998                void             InputExtrude(int enum_type,int start);
    … …  
    144143                int            NumberofNodesPressure(void);
    145144                int            NumberofNodesVelocity(void);
    146                 void                            PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};
    147                 void                            PicoUpdateBox(int loopboxid){_error_("not implemented yet");};
    148                 void                            PicoComputeBasalMelt(void){_error_("not implemented yet");};
    149145                void           PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    150146                int            PressureInterpolation();

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

    r23795 r24010  
    8888                bool        IsFaceOnBoundary(void){_error_("not implemented yet");};
    8989                bool               IsIcefront(void);
    90                 void               Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");};
    9190                bool        IsNodeOnShelfFromFlags(IssmDouble* flags){_error_("not implemented yet");};
    9291                bool        IsOnBase(){_error_("not implemented yet");};
    … …  
    130129                int         NumberofNodesPressure(void){_error_("not implemented yet");};
    131130                int         NumberofNodesVelocity(void){_error_("not implemented yet");};
    132                 void        PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};
    133                 void                    PicoUpdateBox(int loopboxid){_error_("not implemented yet");};
    134                 void                    PicoComputeBasalMelt(void){_error_("not implemented yet");};
    135131                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
    136132                int         PressureInterpolation(void){_error_("not implemented yet");};

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

    r23795 r24010  
    8888                IssmDouble  IceVolume(bool scaled){_error_("not implemented yet");};
    8989                IssmDouble  IceVolumeAboveFloatation(bool scaled){_error_("not implemented yet");};
    90                 void                    Ismip6FloatingiceMeltingRate(void){_error_("not implemented yet");};
    9190                bool        IsFaceOnBoundary(void){_error_("not implemented yet");};
    9291                bool               IsIcefront(void);
    … …  
    138137                int         NumberofNodesPressure(void);
    139138                int         NumberofNodesVelocity(void);
    140                 void                    PicoUpdateBoxid(int* pmax_boxid_basin){_error_("not implemented yet");};       
    141                 void                    PicoUpdateBox(int loopboxid){_error_("not implemented yet");};
    142                 void                    PicoComputeBasalMelt(void){_error_("not implemented yet");};
    143139                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding){_error_("not implemented yet");};
    144140                int         PressureInterpolation(void);

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

    r23999 r24010  
    27262726        return isicefront;
    27272727}/*}}}*/
    2728 void       Tria::Ismip6FloatingiceMeltingRate(void){/*{{{*/
    2729 
    2730         if(!this->IsIceInElement() || !this->IsFloating()) return;
    2731 
    2732         int         basinid,num_basins,M,N;
    2733         IssmDouble  tf,gamma0,base,delta_t_basin,mean_tf_basin,absval;
    2734         IssmDouble  basalmeltrate[NUMVERTICES];
    2735         bool        islocal;
    2736         IssmDouble* delta_t = NULL;
    2737         IssmDouble* mean_tf = NULL;
    2738         IssmDouble* depths  = NULL;
    2739 
    2740         /*Get variables*/
    2741         IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum);
    2742         IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum);
    2743         IssmDouble lf   = this->FindParam(MaterialsLatentheatEnum);
    2744         IssmDouble cp   = this->FindParam(MaterialsMixedLayerCapacityEnum);
    2745 
    2746         /*Hard code sea water density to be consistent with ISMIP6 documentation*/
    2747         rhow = 1028.;
    2748 
    2749         /* Get parameters and inputs */
    2750         this->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
    2751         this->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum);
    2752         this->parameters->FindParam(&gamma0,BasalforcingsIsmp6Gamma0Enum);
    2753         this->parameters->FindParam(&delta_t,&M,BasalforcingsIsmp6DeltaTEnum);    _assert_(M==num_basins);
    2754         this->parameters->FindParam(&islocal,BasalforcingsIsmp6IsLocalEnum);
    2755         if(!islocal) {
    2756                 this->parameters->FindParam(&mean_tf,&N,BasalforcingsIsmp6AverageTfEnum); _assert_(N==num_basins);
    2757         }
    2758         Input* tf_input=this->GetInput(BasalforcingsIsmp6TfShelfEnum);            _assert_(tf_input);
    2759         delta_t_basin = delta_t[basinid];
    2760         if(!islocal) mean_tf_basin = mean_tf[basinid];
    2761 
    2762         /* Compute melt rate */
    2763         Gauss* gauss=this->NewGauss();
    2764         for(int i=0;i<NUMVERTICES;i++){
    2765                 gauss->GaussVertex(i);
    2766                 tf_input->GetInputValue(&tf,gauss);
    2767                 if(!islocal) {
    2768                         absval = mean_tf_basin+delta_t_basin;
    2769                         if (absval<0) {absval = -1.*absval;}
    2770                         basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*(tf+delta_t_basin)*absval;}
    2771                 else {basalmeltrate[i] = gamma0*pow((rhow*cp)/(rhoi*lf),2)*pow(max(tf+delta_t_basin,0.),2);}
    2772         }
    2773 
    2774         /*Return basal melt rate*/
    2775         this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum);
    2776 
    2777         /*Cleanup and return*/
    2778         delete gauss;
    2779         xDelete<IssmDouble>(delta_t);
    2780         xDelete<IssmDouble>(mean_tf);
    2781         xDelete<IssmDouble>(depths);
    2782 
    2783 }/*}}}*/
    27842728bool       Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
    27852729
    … …  
    32693213}
    32703214/*}}}*/
    3271 void       Tria::PicoUpdateBoxid(int* max_boxid_basin_list){/*{{{*/
    3272 
    3273         if(!this->IsIceInElement() || !this->IsFloating()) return;
    3274 
    3275         int        basin_id;
    3276         IssmDouble dist_gl,dist_cf;
    3277 
    3278         inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
    3279         IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id])+1.;
    3280 
    3281         Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
    3282         Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum);  _assert_(dist_cf_input);
    3283 
    3284         /*Get dist_gl and dist_cf at center of element*/
    3285         Gauss* gauss=this->NewGauss(1); gauss->GaussPoint(0);
    3286         dist_gl_input->GetInputValue(&dist_gl,gauss);
    3287         dist_cf_input->GetInputValue(&dist_cf,gauss);
    3288         delete gauss;
    3289 
    3290         /*Ensure values are positive for floating ice*/
    3291         dist_gl = fabs(dist_gl);
    3292         dist_cf = fabs(dist_cf);
    3293 
    3294         /*Compute relative distance to grounding line*/
    3295         IssmDouble rel_dist_gl=dist_gl/(dist_gl+dist_cf);
    3296 
    3297         /*Assign box numbers based on rel_dist_gl*/
    3298         int boxid = -1;
    3299         for(IssmDouble i=0.;i<boxid_max;i++){
    3300                 IssmDouble lowbound  = 1. -sqrt((boxid_max-i   )/boxid_max);
    3301                 IssmDouble highbound = 1. -sqrt((boxid_max-i-1.)/boxid_max);
    3302                 if(rel_dist_gl>=lowbound && rel_dist_gl<=highbound){
    3303                         boxid=reCast<int>(i);
    3304                         break;
    3305                 }
    3306         }
    3307         if(boxid==-1) _error_("No boxid found for element " << this->Sid() << "!");
    3308 
    3309         this->inputs->AddInput(new IntInput(BasalforcingsPicoBoxIdEnum, boxid));       
    3310 }/*}}}*/
    3311 void       Tria::PicoUpdateBox(int loopboxid){/*{{{*/
    3312 
    3313         if(!this->IsIceInElement() || !this->IsFloating()) return;
    3314 
    3315         int boxid;
    3316         this->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
    3317         if(loopboxid!=boxid) return;
    3318 
    3319         int        basinid, maxbox, num_basins, numnodes, M;
    3320         IssmDouble gamma_T, overturning_coeff, thickness;
    3321         IssmDouble pressure, T_star,p_coeff, q_coeff;
    3322         bool       isplume;
    3323 
    3324         /*Get variables*/
    3325         IssmDouble rhoi       = this->FindParam(MaterialsRhoIceEnum);
    3326         IssmDouble rhow       = this->FindParam(MaterialsRhoSeawaterEnum);
    3327         IssmDouble earth_grav = this->FindParam(ConstantsGEnum);
    3328         IssmDouble rho_star   = 1033.;             // kg/m^3
    3329         IssmDouble nu         = rhoi/rhow;
    3330         IssmDouble latentheat = this->FindParam(MaterialsLatentheatEnum);
    3331         IssmDouble c_p_ocean  = this->FindParam(MaterialsMixedLayerCapacityEnum);
    3332         IssmDouble lambda     = latentheat/c_p_ocean;
    3333         IssmDouble a          = -0.0572;          // K/psu
    3334         IssmDouble b          = 0.0788 + this->FindParam(MaterialsMeltingpointEnum);  //K
    3335         IssmDouble c          = 7.77e-4;
    3336         IssmDouble alpha      = 7.5e-5;           // 1/K
    3337         IssmDouble Beta       = 7.7e-4;           // K
    3338 
    3339         /* Get non-box-specific parameters and inputs */
    3340         this->parameters->FindParam(&num_basins, BasalforcingsPicoNumBasinsEnum);
    3341         this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
    3342         this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
    3343         this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
    3344         this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    3345         this->parameters->FindParam(&isplume, BasalforcingsPicoIsplumeEnum);
    3346         Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
    3347         _assert_(basinid<=num_basins);
    3348 
    3349    IssmDouble* boxareas = NULL;
    3350         this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
    3351         _assert_(M==num_basins*maxbox);
    3352 
    3353         IssmDouble area_boxi        = boxareas[basinid*maxbox+boxid];
    3354         IssmDouble g1               = area_boxi*gamma_T;
    3355 
    3356         IssmDouble basalmeltrates_shelf[NUMVERTICES];
    3357         IssmDouble potential_pressure_melting_point[NUMVERTICES];
    3358         IssmDouble Tocs[NUMVERTICES];
    3359         IssmDouble Socs[NUMVERTICES];
    3360 
    3361         /* First box calculations */
    3362         if(boxid==0){
    3363                 /* Get box1 parameters and inputs */
    3364                 IssmDouble time, toc_farocean, soc_farocean;
    3365                 this->parameters->FindParam(&time,TimeEnum);
    3366                 this->parameters->FindParam(&toc_farocean, basinid, time, BasalforcingsPicoFarOceantemperatureEnum);
    3367                 this->parameters->FindParam(&soc_farocean, basinid, time, BasalforcingsPicoFarOceansalinityEnum);
    3368                 IssmDouble s1 = soc_farocean/(nu*lambda);
    3369                 IssmDouble overturnings[NUMVERTICES];
    3370 
    3371                 /* Start looping on the number of verticies and calculate ocean vars */
    3372                 Gauss* gauss=this->NewGauss();
    3373                 for(int i=0;i<NUMVERTICES;i++){
    3374                         gauss->GaussVertex(i);
    3375                         thickness_input->GetInputValue(&thickness,gauss);
    3376                         pressure = (rhoi*earth_grav*1e-4)*thickness;
    3377                         T_star   = a*soc_farocean+b-c*pressure-toc_farocean;
    3378                         p_coeff  = g1/(overturning_coeff*rho_star*(Beta*s1-alpha));
    3379                         q_coeff  = T_star*(g1/(overturning_coeff*rho_star*(Beta*s1-alpha)));
    3380 
    3381                         /* To avoid negatives under the square root */
    3382                         if((0.25*pow(p_coeff,2)-q_coeff)<0) q_coeff = 0.25*p_coeff*p_coeff;
    3383 
    3384                         Tocs[i] = toc_farocean-(-0.5*p_coeff+sqrt(0.25*pow(p_coeff,2)-q_coeff));
    3385                         Socs[i] = soc_farocean-(soc_farocean/(nu*lambda))*(toc_farocean-Tocs[i]);
    3386                         potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
    3387                         if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
    3388                         overturnings[i] = overturning_coeff*rho_star*(Beta*(soc_farocean-Socs[i])-alpha*(toc_farocean-Tocs[i]));
    3389                 }
    3390 
    3391                 if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
    3392                 this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
    3393                 this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
    3394                 this->AddInput(BasalforcingsPicoSubShelfOceanOverturningEnum,overturnings,P1Enum);
    3395 
    3396                 /*Cleanup and return*/
    3397                 delete gauss;
    3398         }
    3399 
    3400         /* Subsequent box calculations */
    3401         else {
    3402                 /* Get subsequent box parameters and inputs */
    3403                 IssmDouble* toc_weighted_avg         = NULL;
    3404                 IssmDouble* soc_weighted_avg         = NULL;
    3405                 IssmDouble* overturning_weighted_avg = NULL;
    3406                 this->parameters->FindParam(&toc_weighted_avg,&num_basins,BasalforcingsPicoAverageTemperatureEnum);
    3407                 this->parameters->FindParam(&soc_weighted_avg,&num_basins,BasalforcingsPicoAverageSalinityEnum);
    3408                 this->parameters->FindParam(&overturning_weighted_avg,&num_basins,BasalforcingsPicoAverageOverturningEnum);
    3409                 IssmDouble mean_toc                  = toc_weighted_avg[basinid];
    3410                 IssmDouble mean_soc                  = soc_weighted_avg[basinid];
    3411                 IssmDouble mean_overturning          = overturning_weighted_avg[basinid];
    3412                 IssmDouble g2                        = g1/(nu*lambda);
    3413 
    3414                 /* Start looping on the number of verticies and calculate ocean vars */
    3415                 Gauss* gauss=this->NewGauss();
    3416                 for(int i=0;i<NUMVERTICES;i++){
    3417                         gauss->GaussVertex(i);
    3418                         thickness_input->GetInputValue(&thickness,gauss);
    3419                         pressure = (rhoi*earth_grav*1.e-4)*thickness;
    3420                         T_star   = a*mean_soc+b-c*pressure-mean_toc;
    3421                         Tocs[i]  = mean_toc+T_star*(g1/(mean_overturning+g1-g2*a*mean_soc));
    3422                         Socs[i]  = mean_soc-mean_soc*((mean_toc-Tocs[i])/(nu*lambda));
    3423                         potential_pressure_melting_point[i] = a*Socs[i]+b-c*pressure;
    3424                         if(!isplume) basalmeltrates_shelf[i] = (-gamma_T/(nu*lambda))*(potential_pressure_melting_point[i]-Tocs[i]);
    3425                 }
    3426 
    3427                 if(!isplume) this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
    3428                 this->AddInput(BasalforcingsPicoSubShelfOceanTempEnum,Tocs,P1Enum);
    3429                 this->AddInput(BasalforcingsPicoSubShelfOceanSalinityEnum,Socs,P1Enum);
    3430 
    3431                 /*Cleanup and return*/
    3432                 xDelete<IssmDouble>(toc_weighted_avg);
    3433                 xDelete<IssmDouble>(soc_weighted_avg);
    3434                 xDelete<IssmDouble>(overturning_weighted_avg);
    3435                 delete gauss;
    3436         }
    3437 
    3438         /*Cleanup and return*/
    3439         xDelete<IssmDouble>(boxareas);
    3440 }
    3441 /*}}}*/
    3442 void       Tria::PicoComputeBasalMelt(void){/*{{{*/
    3443 
    3444         if(!this->IsIceInElement() || !this->IsFloating()) return;
    3445 
    3446    IssmDouble E0, Cd, CdT, YT, lam1, lam2, lam3, M0, CdTS0, y1, y2, x0;
    3447         IssmDouble Tf_gl, YTS, CdTS, G1, G2, G3, g_alpha, M, l, X_hat, M_hat;
    3448         IssmDouble alpha, zgl, Toc, Soc, z_base, yts, slopex, slopey;
    3449 
    3450         /*Get variables*/
    3451         E0    = 3.6e-2;        //Entrainment coefficient
    3452         Cd    = 2.5e-3;        //Drag coefficient
    3453         CdT   = 1.1e-3;        //Turbulent heat exchange coefficient
    3454         YT    = CdT/sqrt(Cd);  //Heat exchange coefficient
    3455         lam1  = -5.73e-2;      //Freezing point-salinity coefficient (degrees C)
    3456         lam2  = 8.32e-2;       //Freezing point offset (degrees C)
    3457         lam3  = 7.61e-4;       //Freezing point-depth coefficient (K m-1)
    3458         M0    = 10.;           //Melt-rate parameter (m yr-1 C-2)
    3459         CdTS0 = 6e-4;          //Heat exchange parameter
    3460         y1    = 0.545;         //Heat exchange parameter
    3461         y2    = 3.5e-5;        //Heat exchange parameter
    3462         x0    = 0.56;          //Dimentionless scaling factor
    3463 
    3464         /*Define arrays*/
    3465         IssmDouble basalmeltrates_shelf[NUMVERTICES];  //Basal melt-rate
    3466 
    3467         /*Polynomial coefficients*/
    3468         IssmDouble p[12];
    3469         p[0]  =  0.1371330075095435;
    3470         p[1]  =  5.527656234709359E1;
    3471         p[2]  = -8.951812433987858E2;
    3472         p[3]  =  8.927093637594877E3;
    3473         p[4]  = -5.563863123811898E4;
    3474         p[5]  =  2.218596970948727E5;
    3475         p[6]  = -5.820015295669482E5;
    3476         p[7]  =  1.015475347943186E6;
    3477         p[8]  = -1.166290429178556E6;
    3478         p[9]  =  8.466870335320488E5;
    3479         p[10] = -3.520598035764990E5;
    3480         p[11] =  6.387953795485420E4;
    3481 
    3482         /*Get inputs*/
    3483    Input* zgl_input            = this->GetInput(GroundinglineHeightEnum);                     _assert_(zgl_input);
    3484         Input* toc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanTempEnum);      _assert_(toc_input);
    3485         Input* soc_input            = this->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum);  _assert_(soc_input);
    3486         Input* base_input           = this->GetInput(BaseEnum);                                    _assert_(base_input);
    3487         Input* baseslopex_input     = this->GetInput(BaseSlopeXEnum);                              _assert_(baseslopex_input);
    3488         Input* baseslopey_input     = this->GetInput(BaseSlopeYEnum);                              _assert_(baseslopey_input);
    3489         this->FindParam(&yts, ConstantsYtsEnum);
    3490 
    3491         /*Loop over the number of vertices in this element*/
    3492         Gauss* gauss=this->NewGauss();
    3493         for(int i=0;i<NUMVERTICES;i++){
    3494                 gauss->GaussVertex(i);
    3495 
    3496                 /*Get inputs*/
    3497       zgl_input->GetInputValue(&zgl,gauss);
    3498                 toc_input->GetInputValue(&Toc,gauss); //(K)
    3499                 soc_input->GetInputValue(&Soc,gauss);
    3500                 base_input->GetInputValue(&z_base,gauss);
    3501                 baseslopex_input->GetInputValue(&slopex,gauss);
    3502                 baseslopey_input->GetInputValue(&slopey,gauss);
    3503 
    3504                 /*Compute ice shelf base slope (radians)*/
    3505            alpha = atan(sqrt(slopex*slopex + slopey*slopey));
    3506                 if(alpha>=M_PI) alpha = M_PI - 0.001;               //ensure sin(alpha) > 0 for meltrate calculations
    3507        
    3508                 /*Make necessary conversions*/
    3509                 Toc = Toc-273.15;
    3510            if(zgl>z_base) zgl=z_base;
    3511 
    3512                 /*Low bound for Toc to ensure X_hat is between 0 and 1*/
    3513                 if(Toc<lam1*Soc+lam2) Toc=lam1*Soc+lam2;
    3514 
    3515                 /*Compute parameters needed for melt-rate calculation*/
    3516                 Tf_gl = lam1*Soc+lam2+lam3*zgl;                                              //Characteristic freezing point
    3517                 YTS = YT*(y1+y2*(((Toc-Tf_gl)*E0*sin(alpha))/(lam3*(CdTS0+E0*sin(alpha))))); //Effective heat exchange coefficient
    3518                 CdTS = sqrt(Cd)*YTS;                                                         //Heat exchange coefficient
    3519                 G1 = sqrt(sin(alpha)/(Cd+E0*sin(alpha)));                                    //Geometric factor
    3520                 G2 = sqrt(CdTS/(CdTS+E0*sin(alpha)));                                        //Geometric factor
    3521                 G3 = (E0*sin(alpha))/(CdTS+E0*sin(alpha));                                   //Geometric factor
    3522       g_alpha = G1*G2*G3;                                                          //Melt scaling factor
    3523                 M = M0*g_alpha*pow((Toc-Tf_gl),2);                                           //Melt-rate scale
    3524                 l = ((Toc-Tf_gl)*(x0*CdTS+E0*sin(alpha)))/(lam3*x0*(CdTS+E0*sin(alpha)));    //Length scale
    3525                 X_hat = (z_base-zgl)/l;                                                      //Dimentionless coordinate system
    3526 
    3527                 /*Compute polynomial fit*/
    3528                 M_hat = 0.;                                                                  //Reset summation variable for each node
    3529                 for(int ii=0;ii<12;ii++) {
    3530                  M_hat += p[ii]*pow(X_hat,ii);                                               //Polynomial fit
    3531            }
    3532 
    3533                 /*Compute melt-rate*/
    3534                 basalmeltrates_shelf[i] = (M*M_hat)/yts;                                     //Basal melt-rate (m/s)
    3535         }
    3536          
    3537         /*Save computed melt-rate*/
    3538    this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrates_shelf,P1Enum);
    3539 
    3540    /*Cleanup and return*/
    3541         delete gauss;
    3542 }/*}}}*/
    35433215void       Tria::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
    35443216

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

    r23993 r24010  
    9898                IssmDouble  IcefrontMassFluxLevelset(bool scaled);
    9999                IssmDouble  GroundinglineMassFlux(bool scaled);
    100                 void        Ismip6FloatingiceMeltingRate(void);
    101100                void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
    102101                void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
    … …  
    116115                int         NumberofNodesPressure(void);
    117116                int         NumberofNodesVelocity(void);
    118                 void        PicoUpdateBoxid(int* pmax_boxid_basin);
    119                 void        PicoUpdateBox(int loopboxid);
    120                 void        PicoComputeBasalMelt(void);
    121117                void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
    122118                int         PressureInterpolation();

  • TabularUnified issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp

    r23255 r24010  
    1212        bool isplume;
    1313
    14    /*First, reset all melt to 0 */
    15    for(int i=0;i<femmodel->elements->Size();i++){
    16               Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    17               int numvertices = element->GetNumberOfVertices();
    18               IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
    19               element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
    20               xDelete<IssmDouble>(values);
    21            }
     14        /*First, reset all melt to 0 */
     15        for(int i=0;i<femmodel->elements->Size();i++){
     16                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     17                int numvertices = element->GetNumberOfVertices();
     18                IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
     19                element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
     20                xDelete<IssmDouble>(values);
     21        }
    2222
    2323        /*PICO melt rate parameterization (Reese et al., 2018)*/
    24    femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
    25    UpdateBoxIdsPico(femmodel);
    26    ComputeBoxAreasPico(femmodel);
    27    for(int i=0;i<maxbox;i++){
    28               UpdateBoxPico(femmodel,i);
    29               ComputeAverageOceanvarsPico(femmodel,i);
     24        femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
     25        UpdateBoxIdsPico(femmodel);
     26        ComputeBoxAreasPico(femmodel);
     27        for(int i=0;i<maxbox;i++){
     28                UpdateBoxPico(femmodel,i);
     29                ComputeAverageOceanvarsPico(femmodel,i);
    3030        }
    3131
    … …  
    103103void ComputeBoxAreasPico(FemModel* femmodel){/*{{{*/
    104104
    105         int num_basins,maxbox,basinid,boxid;
     105        int num_basins,maxbox,basinid,boxid,domaintype;
    106106        IssmDouble dist_max,area;
    107107
    … …  
    113113        for(int i=0;i<femmodel->elements->Size();i++){
    114114                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    115                 if(!element->IsIceInElement() || !element->IsFloating()) continue;
    116                 element->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
    117                 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    118                 boxareas[basinid*maxbox+boxid]+=element->GetHorizontalSurfaceArea();
     115                if(!element->IsOnBase()) continue;
     116                Element* basalelement = element->SpawnBasalElement();
     117                if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
     118                basalelement->inputs->GetInputValue(&boxid,BasalforcingsPicoBoxIdEnum);
     119                basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
     120                boxareas[basinid*maxbox+boxid]+=basalelement->GetHorizontalSurfaceArea();
     121                basalelement->FindParam(&domaintype,DomainTypeEnum);
     122                if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    119123        }
    120124
    … …  
    140144void ComputeAverageOceanvarsPico(FemModel* femmodel, int boxid){/*{{{*/
    141145
    142         int num_basins, basinid, maxbox, M;
     146        int num_basins, basinid, maxbox, M, domaintype;
    143147        IssmDouble area, toc, soc, overturning;
    144148        IssmDouble* boxareas=NULL;
    … …  
    157161        for(int i=0;i<femmodel->elements->Size();i++){
    158162                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    159                 if(!element->IsIceInElement() || !element->IsFloating()) continue;
     163                if(!element->IsOnBase()) continue;
     164                Element* basalelement = element->SpawnBasalElement();
     165                if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
    160166                int el_boxid;
    161                 element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
     167                basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
    162168                if(el_boxid!=boxid) continue;
    163169
    164                 Input* tocs_input=element->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input);
    165                 Input* socs_input=element->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
    166 
    167                 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    168                 Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
     170                Input* tocs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanTempEnum); _assert_(tocs_input);
     171                Input* socs_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanSalinityEnum); _assert_(socs_input);
     172
     173                basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
     174                Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
    169175                tocs_input->GetInputValue(&toc,gauss);
    170176                socs_input->GetInputValue(&soc,gauss);
    171177                delete gauss;
    172                 area=element->GetHorizontalSurfaceArea();
     178                area=basalelement->GetHorizontalSurfaceArea();
    173179                toc_weighted_avg[basinid]+=toc*area;
    174180                soc_weighted_avg[basinid]+=soc*area;
     181                basalelement->FindParam(&domaintype,DomainTypeEnum);
     182                if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    175183        }
    176184
    … …  
    194202                for(int i=0;i<femmodel->elements->Size();i++){
    195203                        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    196                         if(!element->IsIceInElement() || !element->IsFloating()) continue;
     204                        if(!element->IsOnBase()) continue;
     205                        Element* basalelement = element->SpawnBasalElement();
     206                        if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
    197207                        int el_boxid;
    198                         element->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
     208                        basalelement->inputs->GetInputValue(&el_boxid,BasalforcingsPicoBoxIdEnum);
    199209                        if(el_boxid!=boxid) continue;
    200210
    201                 Input* overturnings_input=element->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
    202 
    203                         element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
    204                         Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
     211                Input* overturnings_input=basalelement->GetInput(BasalforcingsPicoSubShelfOceanOverturningEnum); _assert_(overturnings_input);
     212
     213                        basalelement->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
     214                        Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
    205215                        overturnings_input->GetInputValue(&overturning,gauss);
    206216                        delete gauss;
    207                         area=element->GetHorizontalSurfaceArea();
     217                        area=basalelement->GetHorizontalSurfaceArea();
    208218                        overturning_weighted_avg[basinid]+=overturning*area;
     219                        basalelement->FindParam(&domaintype,DomainTypeEnum);
     220                        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    209221                }
    210222

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

    r23982 r24010  
    7373void FloatingiceMeltingRateIsmip6x(FemModel* femmodel){/*{{{*/
    7474
    75         int         num_basins, basinid,num_depths;
     75        int         num_basins, basinid, num_depths, domaintype;
    7676        IssmDouble  area, tf, base, time;
    7777        bool        islocal;
    … …  
    154154                for(int i=0;i<femmodel->elements->Size();i++){
    155155                        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    156                         if(!element->IsIceInElement() || !element->IsFloating()) continue;
    157                         Input* tf_input=element->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input);
    158                         element->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
    159                         Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
     156                        if(!element->IsOnBase()) continue;
     157                        /*Spawn basal element if on base to compute element area*/
     158                        Element* basalelement = element->SpawnBasalElement();
     159                        if(!basalelement->IsIceInElement() || !basalelement->IsFloating()) continue;
     160                        Input* tf_input=basalelement->GetInput(BasalforcingsIsmp6TfShelfEnum); _assert_(tf_input);
     161                        basalelement->inputs->GetInputValue(&basinid,BasalforcingsIsmp6BasinIdEnum);
     162                        Gauss* gauss=basalelement->NewGauss(1); gauss->GaussPoint(0);
    160163                        tf_input->GetInputValue(&tf,gauss);
    161164                        delete gauss;
    162                         area=element->GetHorizontalSurfaceArea();
     165                        area=basalelement->GetHorizontalSurfaceArea();
    163166                        tf_weighted_avg[basinid]+=tf*area;
    164167                        areas_summed[basinid]   +=area;
     168                        /*Delete spawned element if we are in 3D*/
     169                        basalelement->FindParam(&domaintype,DomainTypeEnum);
     170                        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    165171                }
    166172
    … …  
    174180        }
    175181
     182   /*Compute meltrates*/
     183        for(int i=0;i<femmodel->elements->Size();i++){
     184                Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     185                element->Ismip6FloatingiceMeltingRate();
     186        }
     187
    176188        /*Cleanup and return */
    177189        xDelete<IssmDouble>(tf_weighted_avg);
    … …  
    180192        xDelete<IssmDouble>(areas_summed_cpu);
    181193        xDelete<IssmDouble>(tf_depths);
    182 
    183    /*Compute meltrates*/
    184         for(int i=0;i<femmodel->elements->Size();i++){
    185                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    186                 element->Ismip6FloatingiceMeltingRate();
    187         }
    188194}
    189195/*}}}*/

  • TabularUnified issm/trunk-jpl/src/m/classes/basalforcingsismip6.m

    r23906 r24010  
    1919                function self = extrude(self,md) % {{{
    2020                        self.basin_id=project3d(md,'vector',self.basin_id,'type','element','layer',1);
    21                         self.tf=project3d(md,'vector',self.tf,'type','element','layer',1);
    22                         self.delta_t=project3d(md,'vector',self.delta_t,'type','element','layer',1);
     21                        %self.tf=project3d(md,'vector',self.tf,'type','element','layer',1);
     22                        %self.delta_t=project3d(md,'vector',self.delta_t,'type','element','layer',1);
     23                        self.tf=project3d(md,'vector',self.tf,'type','node');
    2324                        self.geothermalflux=project3d(md,'vector',self.geothermalflux,'type','element','layer',1); %bedrock only gets geothermal flux
    2425                        self.groundedice_melting_rate=project3d(md,'vector',self.groundedice_melting_rate,'type','node','layer',1);

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

    r22918 r24010  
    1313                farocean_temperature      = NaN;
    1414                farocean_salinity         = NaN;
    15                 isplume                   = NaN;
     15                isplume                   = 0;
    1616                geothermalflux            = NaN;
    1717                groundedice_melting_rate  = NaN;
    … …  
    6868                                md = checkfield(md,'fieldname','basalforcings.overturning_coeff','numel',1,'NaN',1,'Inf',1,'>',0);
    6969                                md = checkfield(md,'fieldname','basalforcings.gamma_T','numel',1,'NaN',1,'Inf',1,'>',0);
    70                                 md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]);
     70                                md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'size',[md.basalforcings.num_basins+1 NaN]);
    7171                                md = checkfield(md,'fieldname','basalforcings.farocean_salinity','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]);
    7272                                md = checkfield(md,'fieldname','basalforcings.isplume','values',[0 1]);
Note: See TracChangeset for help on using the changeset viewer.