Changeset 23797


Ignore:
Timestamp:
03/13/19 15:22:07 (6 years ago)
Author:
Mathieu Morlighem
Message:

NEW: finished Ismip6 melt implementation

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

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp

    r23795 r23797  
    191191                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    192192                                for(int kk=0;kk<K;kk++){
    193                                         element->DatasetInputAdd(BasalforcingsIsmp6TfEnum,array3d[kk],iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmp6TfEnum,7,BasalforcingsIsmp6TfEnum);
     193                                        element->DatasetInputAdd(BasalforcingsIsmp6TfEnum,array3d[kk],iomodel,Ms[kk],Ns[kk],1,BasalforcingsIsmp6TfEnum,7,kk);
    194194                                }
    195195                        }
    … …  
    536536        element->FindParam(&melt_style,GroundinglineMeltInterpolationEnum);
    537537        element->FindParam(&dt,TimesteppingTimeStepEnum);
    538         Input* gmb_input           = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
    539         Input* fmb_input           = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
    540         Input* gllevelset_input   = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
    541         Input* ms_input            = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
    542         Input* thickness_input     = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
     538        Input* gmb_input        = element->GetInput(BasalforcingsGroundediceMeltingRateEnum);  _assert_(gmb_input);
     539        Input* fmb_input        = element->GetInput(BasalforcingsFloatingiceMeltingRateEnum);  _assert_(fmb_input);
     540        Input* gllevelset_input = element->GetInput(MaskGroundediceLevelsetEnum);              _assert_(gllevelset_input);
     541        Input* ms_input         = element->GetInput(SmbMassBalanceEnum);                       _assert_(ms_input);
     542        Input* thickness_input  = element->GetInput(ThicknessEnum);                            _assert_(thickness_input);
    543543
    544544        /*Recover portion of element that is grounded*/

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

    r23644 r23797  
    19011901}
    19021902/*}}}*/
    1903 void       Element::DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code,int input_enum){/*{{{*/
     1903void       Element::DatasetInputAdd(int enum_type,IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int input_enum,int code,int input_id){/*{{{*/
    19041904    /*enum_type: the name of the DatasetInput (eg Outputdefinition1)
    19051905          * vector: information being stored (eg observations)
    19061906          * vector_type: is if by element or by vertex
    1907           * vector_enum: is the name of the vector being stored
     1907          * input_enum: is the name of the vector being stored
    19081908          * code: what type of data is in the vector (booleans, ints, doubles)
    19091909          */
    … …  
    19421942                          values[0]=vector[0];
    19431943                                switch(this->ObjectEnum()){
    1944                     case TriaEnum:  datasetinput->AddInput(new TriaInput(vector_enum,values,P0Enum),input_enum); break;
    1945                     case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P0Enum),input_enum); break;
    1946                     case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P0Enum),input_enum); break;
     1944                    case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break;
     1945                    case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break;
     1946                    case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break;
    19471947                    default: _error_("Not implemented yet");
    19481948                                }
    … …  
    19511951            for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
    19521952                                switch(this->ObjectEnum()){
    1953                     case TriaEnum:  datasetinput->AddInput(new TriaInput(vector_enum,values,P1Enum),input_enum); break;
    1954                     case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1Enum),input_enum); break;
    1955                     case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1Enum),input_enum); break;
     1953                    case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
     1954                    case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break;
     1955                    case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break;
    19561956                    default: _error_("Not implemented yet");
    19571957                                }  }
    … …  
    19601960            IssmDouble* times = xNew<IssmDouble>(N);
    19611961            for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    1962             TransientInput* transientinput=new TransientInput(vector_enum,times,N);
     1962            TransientInput* transientinput=new TransientInput(input_enum,times,N);
    19631963            for(t=0;t<N;t++){
    19641964                for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
    19651965                switch(this->ObjectEnum()){
    1966                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
    1967                     case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break;
    1968                     case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break;
     1966                    case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,values,P1Enum)); break;
     1967                    case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,values,P1Enum)); break;
     1968                    case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,values,P1Enum)); break;
    19691969                    default: _error_("Not implemented yet");
    19701970                }
    19711971            }
    1972             datasetinput->AddInput(transientinput,input_enum);
     1972            datasetinput->AddInput(transientinput,input_id);
    19731973            xDelete<IssmDouble>(times);
    19741974        }
    … …  
    19821982                          if     (N==this->GetNumberOfNodes(P1Enum)   ){
    19831983                                  switch(this->ObjectEnum()){
    1984                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(vector_enum,values,P1Enum),input_enum); break;
    1985                                           case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1Enum),input_enum); break;
    1986                                           case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1Enum),input_enum); break;
     1984                                          case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1Enum),input_id); break;
     1985                                          case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1Enum),input_id); break;
     1986                                          case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1Enum),input_id); break;
    19871987                                          default: _error_("Not implemented yet");
    19881988                                  }
    … …  
    19901990                          else if(N==this->GetNumberOfNodes(P0Enum)   ){
    19911991                                  switch(this->ObjectEnum()){
    1992                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(vector_enum,values,P0Enum),input_enum); break;
    1993                                           case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P0Enum),input_enum); break;
    1994                                           case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P0Enum),input_enum); break;
     1992                                          case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P0Enum),input_id); break;
     1993                                          case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P0Enum),input_id); break;
     1994                                          case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P0Enum),input_id); break;
    19951995                                          default: _error_("Not implemented yet");
    19961996                                  }
    … …  
    19981998                          else if(N==this->GetNumberOfNodes(P1xP2Enum)){
    19991999                                  switch(this->ObjectEnum()){
    2000                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(vector_enum,values,P1xP2Enum),input_enum); break;
    2001                                           case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1xP2Enum),input_enum); break;
    2002                                           case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1xP2Enum),input_enum); break;
     2000                                          case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1xP2Enum),input_id); break;
     2001                                          case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP2Enum),input_id); break;
     2002                                          case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP2Enum),input_id); break;
    20032003                                          default: _error_("Not implemented yet");
    20042004                                  }
    … …  
    20062006                          else if(N==this->GetNumberOfNodes(P1xP3Enum)) {
    20072007                                 switch(this->ObjectEnum()){
    2008                                           case TriaEnum:  datasetinput->AddInput(new TriaInput(vector_enum,values,P1xP3Enum),input_enum); break;
    2009                                           case PentaEnum: datasetinput->AddInput(new PentaInput(vector_enum,values,P1xP3Enum),input_enum); break;
    2010                                           case TetraEnum: datasetinput->AddInput(new TetraInput(vector_enum,values,P1xP3Enum),input_enum); break;
     2008                                          case TriaEnum:  datasetinput->AddInput(new TriaInput(input_enum,values,P1xP3Enum),input_id); break;
     2009                                          case PentaEnum: datasetinput->AddInput(new PentaInput(input_enum,values,P1xP3Enum),input_id); break;
     2010                                          case TetraEnum: datasetinput->AddInput(new TetraInput(input_enum,values,P1xP3Enum),input_id); break;
    20112011                                          default: _error_("Not implemented yet");
    20122012                                  }
    … …  
    20162016                  }
    20172017                  else{
    2018                           _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     2018                          _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
    20192019                  }
    20202020
    … …  
    20292029        if(M==iomodel->numberofelements){
    20302030            if (code==5){ //boolean
    2031                 datasetinput->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->Sid()])),input_enum);
     2031                datasetinput->AddInput(new BoolInput(input_enum,reCast<bool>(vector[this->Sid()])),input_id);
    20322032            }
    20332033            else if (code==6){ //integer
    2034                 datasetinput->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->Sid()])),input_enum);
     2034                datasetinput->AddInput(new IntInput(input_enum,reCast<int>(vector[this->Sid()])),input_id);
    20352035            }
    20362036            else if (code==7){ //IssmDouble
    2037                 datasetinput->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]),input_enum);
     2037                datasetinput->AddInput(new DoubleInput(input_enum,vector[this->Sid()]),input_id);
    20382038            }
    20392039            else _error_("could not recognize nature of vector from code " << code);
    … …  
    20432043            IssmDouble* times = xNew<IssmDouble>(N);
    20442044            for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
    2045             TransientInput* transientinput=new TransientInput(vector_enum,times,N);
     2045            TransientInput* transientinput=new TransientInput(input_enum,times,N);
    20462046            TriaInput* bof=NULL;
    20472047            for(t=0;t<N;t++){
    20482048                value=vector[N*this->Sid()+t];
    20492049                switch(this->ObjectEnum()){
    2050                     case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,&value,P0Enum)); break;
    2051                     case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,&value,P0Enum)); break;
    2052                     case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,&value,P0Enum)); break;
     2050                    case TriaEnum:  transientinput->AddTimeInput(new TriaInput( input_enum,&value,P0Enum)); break;
     2051                    case PentaEnum: transientinput->AddTimeInput(new PentaInput(input_enum,&value,P0Enum)); break;
     2052                    case TetraEnum: transientinput->AddTimeInput(new TetraInput(input_enum,&value,P0Enum)); break;
    20532053                    default: _error_("Not implemented yet");
    20542054                }
    20552055            }
    2056             datasetinput->AddInput(transientinput,input_enum);
     2056            datasetinput->AddInput(transientinput,input_id);
    20572057            xDelete<IssmDouble>(times);
    20582058        }
    2059         else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     2059        else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
    20602060    }
    20612061    else if(vector_type==3){ //element vector
    … …  
    20662066            IssmDouble* layers = xNewZeroInit<IssmDouble>(N);;
    20672067            for(t=0;t<N;t++) layers[t] = vector[N*this->Sid()+t];
    2068             DoubleArrayInput* arrayinput=new DoubleArrayInput(vector_enum,layers,N);
    2069             datasetinput->AddInput(arrayinput,input_enum);
     2068            DoubleArrayInput* arrayinput=new DoubleArrayInput(input_enum,layers,N);
     2069            datasetinput->AddInput(arrayinput,input_id);
    20702070            xDelete<IssmDouble>(layers);
    20712071        }
    2072         else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");
     2072        else _error_("element vector is either numberofelements or numberofelements+1 long. Field provided (" << EnumToStringx(input_enum) << ") is " << M << " long");
    20732073    }
    20742074    else _error_("Cannot add input for vector type " << vector_type << " (not supported)");

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

    r23795 r23797  
    24152415
    24162416        /*Get variables*/
    2417         IssmDouble rhoi     = this->FindParam(MaterialsRhoIceEnum);
    2418         IssmDouble rhow     = this->FindParam(MaterialsRhoSeawaterEnum)+5;
    2419         IssmDouble lf       = this->FindParam(MaterialsLatentheatEnum);
    2420         IssmDouble cp       = this->FindParam(MaterialsMixedLayerCapacityEnum);
     2417        IssmDouble rhoi = this->FindParam(MaterialsRhoIceEnum);
     2418        IssmDouble rhow = this->FindParam(MaterialsRhoSeawaterEnum);
     2419        IssmDouble lf   = this->FindParam(MaterialsLatentheatEnum);
     2420        IssmDouble cp   = this->FindParam(MaterialsMixedLayerCapacityEnum);
     2421
     2422        /*Hard code sea water density to be consistent with ISMIP6 documentation*/
     2423        rhow = 1028.;
    24212424
    24222425        /* Get parameters and inputs */
    … …  
    24292432        delta_t_basin = delta_t[basinid]; mean_tf_basin = mean_tf[basinid];
    24302433
    2431    /* Compute melt rate */
    2432    Gauss* gauss=this->NewGauss();
     2434        /* Compute melt rate */
     2435        Gauss* gauss=this->NewGauss();
    24332436        for(int i=0;i<NUMVERTICES;i++){
    24342437                gauss->GaussVertex(i);
    … …  
    24392442        }
    24402443
    2441         cout << "meltrate = " << basalmeltrate[0] << endl;
    2442 
    24432444        /*Return basal melt rate*/
    24442445        this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,basalmeltrate,P1Enum);
    2445         _error_("STOP HERE");
    24462446
    24472447        /*Cleanup and return*/
    … …  
    24512451        xDelete<IssmDouble>(depths);
    24522452
    2453         _error_("not implemented yet");
    2454 
    2455         }/*}}}*/
     2453}/*}}}*/
    24562454bool       Tria::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
    24572455

  • issm/trunk-jpl/src/c/classes/Inputs/DatasetInput.cpp

    r23518 r23797  
    157157/*Object functions*/
    158158void DatasetInput::Configure(Parameters* parameters){/*{{{*/
    159         /*do nothing: */
     159
     160        for(int i=0;i<this->numids;i++){
     161                Input* input=xDynamicCast<Input*>(this->inputs->GetObjectByOffset(i));
     162                input->Configure(parameters);
     163        }
     164
    160165}
    161166/*}}}*/

  • issm/trunk-jpl/src/c/classes/Inputs/TransientInput.cpp

    r23546 r23797  
    256256int  TransientInput::GetTimeInputOffset(IssmDouble time){/*{{{*/
    257257
    258         int        found;
    259         int        offset;
     258        int offset;
    260259
    261260        /*go through the timesteps, and figure out which interval we
    262261         *     *fall within. Then interpolate the values on this interval: */
    263         found=binary_search(&offset,time,this->timesteps,this->numtimesteps);
     262        int found=binary_search(&offset,time,this->timesteps,this->numtimesteps);
    264263        if(!found) _error_("Input not found (is TransientInput sorted ?)");
    265264

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

    r23795 r23797  
    7373void FloatingiceMeltingRateIsmip6x(FemModel* femmodel){/*{{{*/
    7474
     75        int         num_basins, basinid,num_depths;
     76        IssmDouble  area, tf, base, time;
     77        IssmDouble* tf_depths = NULL;
    7578
    76         int         num_basins, basinid;
    77         IssmDouble  area, tf, base, time;
    78         IssmDouble  tf_test[3];
     79        femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum);
     80        femmodel->parameters->FindParam(&tf_depths,&num_depths,BasalforcingsIsmp6TfDepthsEnum); _assert_(tf_depths);
    7981
    80    femmodel->parameters->FindParam(&num_basins,BasalforcingsIsmp6NumBasinsEnum);
     82        /*Binary search works for vectors that are sorted in increasing order only, make depths positive*/
     83        for(int i=0;i<num_depths;i++) tf_depths[i] = -tf_depths[i];
    8184
    82    IssmDouble* tf_weighted_avg     = xNewZeroInit<IssmDouble>(num_basins);
     85        IssmDouble* tf_weighted_avg     = xNewZeroInit<IssmDouble>(num_basins);
    8386        IssmDouble* tf_weighted_avg_cpu = xNewZeroInit<IssmDouble>(num_basins);
    84    IssmDouble* areas_summed        = xNewZeroInit<IssmDouble>(num_basins);
    85    IssmDouble* areas_summed_cpu    = xNewZeroInit<IssmDouble>(num_basins);
    86 
    87         /*Find and save TF at each ice shelf point*/
    88         //      element->parameters->FindParam(&time,TimeEnum);
    89         //      Input* tf_input = element->GetInput(BasalforcingsIsmp6TfEnum); _assert_(tf_input);
    90         //      tf_input->Echo();
    91         //      Input* base_input = element->GetInput(BaseEnum);               _assert_(base_input);
    92         //      Gauss* gauss=element->NewGauss(1); gauss->GaussPoint(0);
    93         //      base_input->GetInputValue(&base,gauss);
    94         //      //tf_input->GetInputValue(&tf,gauss,BasalforcingsIsmp6TfEnum);
     87        IssmDouble* areas_summed        = xNewZeroInit<IssmDouble>(num_basins);
     88        IssmDouble* areas_summed_cpu    = xNewZeroInit<IssmDouble>(num_basins);
    9589
    9690        /*TEST: Set tf=2 for all ice shelf elements*/
    9791        for(int i=0;i<femmodel->elements->Size();i++){
    98                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    99                 if(!element->IsIceInElement() || !element->IsFloating()) continue;
    100                 for(int k=0;k<3;k++) {tf_test[k] = 2.;}
     92                Element* element     = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     93                int      numvertices = element->GetNumberOfVertices();
     94
     95                /*Set melt to 0 if non floating*/
     96                if(!element->IsIceInElement() || !element->IsFloating()){
     97                        IssmDouble* values = xNewZeroInit<IssmDouble>(numvertices);
     98                        element->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);
     99                        element->AddInput(BasalforcingsIsmp6TfShelfEnum,values,P1Enum);
     100                        xDelete<IssmDouble>(values);
     101                        continue;
     102                }
     103
     104                /*Get TF on all vertices*/
     105                IssmDouble* tf_test        = xNew<IssmDouble>(numvertices);
     106                IssmDouble* depth_vertices = xNew<IssmDouble>(numvertices);
     107                Input*      tf_input = element->GetInput(BasalforcingsIsmp6TfEnum); _assert_(tf_input);
     108
     109                element->GetInputListOnVertices(&depth_vertices[0],BaseEnum);
     110
     111                Gauss* gauss=element->NewGauss();
     112                for(int iv=0;iv<numvertices;iv++){
     113                        gauss->GaussVertex(iv);
     114
     115                        /*Find out where the ice shelf base is within tf_depths*/
     116                        IssmDouble depth = -depth_vertices[iv]; /*NOTE: make sure we are dealing with depth>0*/
     117                        int offset;
     118                        int found=binary_search(&offset,depth,tf_depths,num_depths);
     119                        if(!found) _error_("depth not found");
     120
     121                        if (offset==-1){
     122                                /*get values for the first depth: */
     123                                _assert_(depth<=tf_depths[0]);
     124                                tf_input->GetInputValue(&tf_test[iv],gauss,0);
     125                        }
     126                        else if(offset==num_depths-1){
     127                                /*get values for the last time: */
     128                                _assert_(depth>=tf_depths[num_depths-1]);
     129                                tf_input->GetInputValue(&tf_test[iv],gauss,num_depths-1);
     130                        }
     131                        else {
     132                                /*get values between two times [offset:offset+1[, Interpolate linearly*/
     133                                _assert_(depth>=tf_depths[offset] && depth<tf_depths[offset+1]);
     134                                IssmDouble deltaz=tf_depths[offset+1]-tf_depths[offset];
     135                                IssmDouble alpha2=(depth-tf_depths[offset])/deltaz;
     136                                IssmDouble alpha1=(1.-alpha2);
     137                                IssmDouble tf1,tf2;
     138                                tf_input->GetInputValue(&tf1,gauss,offset);
     139                                tf_input->GetInputValue(&tf2,gauss,offset+1);
     140                                tf_test[iv] = alpha1*tf1 + alpha2*tf2;
     141                        }
     142                }
     143
    101144                element->AddInput(BasalforcingsIsmp6TfShelfEnum,tf_test,P1Enum);
     145                xDelete<IssmDouble>(tf_test);
     146                xDelete<IssmDouble>(depth_vertices);
    102147        }
    103148
    … …  
    129174        xDelete<IssmDouble>(areas_summed);
    130175        xDelete<IssmDouble>(areas_summed_cpu);
     176        xDelete<IssmDouble>(tf_depths);
    131177
    132178   /*Compute meltrates*/
Note: See TracChangeset for help on using the changeset viewer.