source: issm/oecreview/Archive/22819-23185/ISSM-23065-23066.diff

Last change on this file was 23186, checked in by Mathieu Morlighem, 7 years ago

CHG: added Archive/22819-23185
File size: 181.5 KB

  • ../trunk-jpl/src/c/modules/SetActiveNodesLSMx/SetActiveNodesLSMx.cpp

     
    9191
    9292        /* Intermediaries */
    9393        int numvertices = element->GetNumberOfVertices();
    94        
     94
    9595        if(element->IsIceInElement()){
    9696                for(int i = 0;i<numvertices;i++){
    9797                        vec_mask_ice->SetValue(element->vertices[i]->Sid(),1.,INS_VAL);

  • ../trunk-jpl/src/c/modules/OutputDefinitionsResponsex/OutputDefinitionsResponsex.cpp

     
    2121
    2222                        /*This is the object that we have been chasing for. compute the response and return: */
    2323                        IssmDouble return_value=definition->Response(femmodel);
    24                
     24
    2525                        /*cleanup: */
    2626                        xDelete<char>(name);
    2727
    … …  
    3030                }
    3131                xDelete<char>(name);
    3232        }
    33        
     33
    3434        /*If we are here, did not find the definition for this response, not good!: */
    3535        _error_("Could not find the response for output definition " << output_string << " because could not find the definition itself!");
    3636
    … …  
    4343
    4444        /*Now, go through the output definitions, and retrieve the object which corresponds to our requested response, output_enum: */
    4545        for(int i=0;i<output_definitions->Size();i++){
    46                
     46
    4747                //Definition* definition=xDynamicCast<Definition*>(output_definitions->GetObjectByOffset(i));
    4848                Definition* definition=dynamic_cast<Definition*>(output_definitions->GetObjectByOffset(i));
    4949
    … …  
    5252
    5353                        /*This is the object that we have been chasing for. compute the response and return: */
    5454                        IssmDouble return_value=definition->Response(femmodel);
    55                
     55
    5656                        /*return:*/
    5757                        return return_value;
    5858                }
    5959        }
    60        
     60
    6161        /*If we are here, did not find the definition for this response, not good!: */
    6262        _error_("Could not find the response for output definition " << EnumToStringx(output_enum)
    6363                                <<" ("<<output_enum<<")"

  • ../trunk-jpl/src/c/modules/GroundinglineMigrationx/GroundinglineMigrationx.cpp

     
    1414        IssmDouble         *phi_ungrounding                  = NULL;
    1515        Element            *element                          = NULL;
    1616
    17 
    1817        /*retrieve parameters: */
    1918        parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
    2019        parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    2120
    2221        if(migration_style==NoneEnum) return;
    23        
     22
    2423        if(VerboseModule()) _printf0_("   Migrating grounding line\n");
    2524
    2625        /*Set toolkit to default*/

  • ../trunk-jpl/src/c/modules/AllocateSystemMatricesx/AllocateSystemMatricesx.cpp

     
    7777                if(pdf) df =new Vector<IssmDouble>(flocalsize,fsize);
    7878                if(ppf) pf =new Vector<IssmDouble>(flocalsize,fsize);
    7979        }
    80        
     80
    8181        /*Free ressources: */
    8282        xDelete<char>(toolkittype);
    8383

  • ../trunk-jpl/src/c/modules/GetVectorFromControlInputsx/GetVectorFromControlInputsx.cpp

     
    2727                int size = 0;
    2828                for(int i=0;i<num_controls;i++) size+=M[i]*N[i];
    2929
    30 
    3130                /*2. Allocate vector*/
    3231                vector=new Vector<IssmDouble>(size);
    3332
    … …  
    5958                parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    6059                parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
    6160
    62 
    6361                /*2. Allocate vector*/
    6462                vector=new Vector<IssmDouble>(num_controls*vertices->NumberOfVertices());
    6563

  • ../trunk-jpl/src/c/modules/MeshProfileIntersectionx/MeshProfileIntersectionx.cpp

     
    165165        edge3=SegmentIntersect(&gamma1,&gamma2, xel,yel,xsegment,ysegment);
    166166
    167167        /*edge can be either IntersectEnum (one and only one intersection between the edge and the segment), ColinearEnum (edge and segment are collinear) and SeparateEnum (no intersection): */
    168                
     168
    169169        /*if (el==318 && contouri==9){
    170170                _printf_(edge1 << " " << edge2 << " " << edge3 << " "  << alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
    171        
     171
    172172        _printf_("Bool" << (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum));
    173173        }*/
    174174
    … …  
    198198                }
    199199                segments_dataset->AddObject(new  Segment<double>(el+1,xfinal[0],yfinal[0],xfinal[1],yfinal[1]));
    200200
    201                
    202201                /*This case is impossible: not quite! */
    203202                //_printf_(alpha1 << " " << alpha2 << " " << beta1 << " " << beta2 << " " << gamma1 << " " << gamma2 << " " << xsegment[0] << " "  << xsegment[1] << " " << ysegment[0] << " " << ysegment[1] << " " << xnodes[0] << " " << xnodes[1] << " " << xnodes[2] << " " << ynodes[0] << " " << ynodes[1] << " " << ynodes[2]);
    204203                /* _error_("error: a line cannot go through 3 different vertices!");*/
    … …  
    225224                }
    226225        }
    227226        else if(  (edge1==IntersectEnum) || (edge2==IntersectEnum) || (edge3==IntersectEnum)   ){
    228        
     227
    229228                /*if (el==318 && contouri==9){
    230229                        _printf_("hello" <<  " NodeInElement 0 " << (NodeInElement(xnodes,ynodes,xsegment[0],ysegment[0])) <<  " NodeInElement 1 " << (NodeInElement(xnodes,ynodes,xsegment[1],ysegment[1])));
    231230                }*/
    … …  
    250249                        if (IsIdenticalNode(xnodes[0],ynodes[0],xsegment[0],ysegment[0],tolerance) ||
    251250                                IsIdenticalNode(xnodes[1],ynodes[1],xsegment[0],ysegment[0],tolerance) ||
    252251                                IsIdenticalNode(xnodes[2],ynodes[2],xsegment[0],ysegment[0],tolerance)){
    253                                
     252
    254253                                /*ok, segments[0] is common to one of our vertices: */
    255254                                coord1=0;
    256255                                if(edge1==IntersectEnum){coord2=alpha1;}

  • ../trunk-jpl/src/c/modules/DistanceToMaskBoundaryx/DistanceToMaskBoundaryxt.cpp

     
    3838
    3939                IssmDouble d0=1.e+10;
    4040                IssmDouble xi,yi;
    41                
     41
    4242                //recover vertex position:
    4343                xi=x[i];  yi=y[i];
    4444

  • ../trunk-jpl/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp

     
    8383        int interpolation_type;
    8484        /*this one is special: we don't specify the type, but let the nature of the inputs dictace.
    8585         * P0 -> ElementSIdEnum, P1 ->VertexSIdEnum: */
    86        
     86
    8787        /*We go find the input of the first element, and query its interpolation type: */
    8888        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(0));
    8989        Input* input=element->GetInput(name);

  • ../trunk-jpl/src/c/modules/NodalValuex/NodalValuex.cpp

     
    2020         *element, figure out  if they hold the vertex, and the data. If so, return it: */
    2121        cpu_found=-1;
    2222        found=0;
    23        
     23
    2424        for(int i=0;i<elements->Size();i++){
    2525                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
    2626                found=element->NodalValue(&value,index,natureofdataenum);

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Control/UpdateElementsAndMaterialsControl.cpp

     
    88#include "../../MeshPartitionx/MeshPartitionx.h"
    99#include "../ModelProcessorx.h"
    1010
    11 
    1211#if !defined(_HAVE_ADOLC_)
    1312void    UpdateElementsAndMaterialsControl(Elements* elements,Parameters* parameters,Materials* materials, IoModel* iomodel){
    1413        /*Intermediary*/
    … …  
    2120        char     **controls         = NULL;
    2221        char     **cost_functions   = NULL;
    2322
    24 
    2523        /*Fetch parameters: */
    2624        iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
    2725        if(control_analysis) iomodel->FindConstant(&num_controls,"md.inversion.num_control_parameters");
    … …  
    2826
    2927        /*Now, return if no control*/
    3028        if(!control_analysis) return;
    31        
     29
    3230        /*Process controls and convert from string to enums*/
    3331        iomodel->FindConstant(&controls,&num_controls,"md.inversion.control_parameters");
    3432        if(num_controls<1) _error_("no controls found");
    … …  
    4745        }
    4846
    4947        iomodel->FetchData(3,"md.inversion.cost_functions_coefficients","md.inversion.min_parameters","md.inversion.max_parameters");
    50        
     48
    5149        /*Fetch Observations */
    5250        iomodel->FindConstant(&domaintype,"md.mesh.domain_type");
    5351        for(int i=0;i<num_cost_functions;i++){
    … …  
    153151        _assert_(M==num_independent_objects);
    154152        iomodel->FetchData(&types,NULL,NULL,"md.autodiff.independent_object_types");
    155153
    156                
    157154        M_all = xNew<int>(num_independent_objects);
    158155
    159156        /*create independent objects, and at the same time, fetch the corresponding independent variables,
    … …  
    161158        iomodel->FetchData(&independents_fullmin,&M_par,&N_par,"md.autodiff.independent_min_parameters");
    162159        iomodel->FetchData(&independents_fullmax,&M_par,&N_par,"md.autodiff.independent_max_parameters");
    163160        iomodel->FetchData(&control_sizes,NULL,NULL,"md.autodiff.independent_control_sizes");
    164        
     161
    165162        int* start_point = NULL;
    166163        start_point = xNew<int>(num_independent_objects);
    167164        int counter = 0;
    … …  
    179176                        int   input_enum;
    180177                        IssmDouble*     independents_min                        = NULL;
    181178                        IssmDouble*        independents_max                     = NULL;
    182        
     179
    183180                        FieldAndEnumFromCode(&input_enum,&iofieldname,names[i]);
    184181
    185182                        /*Fetch required data*/
    … …  
    186183                        iomodel->FetchData(&independent,&M,&N,iofieldname);
    187184                        _assert_(independent);
    188185                        _assert_(N==control_sizes[i]);
    189                
     186
    190187                        independents_min = NULL; independents_min = xNew<IssmDouble>(M*N);
    191188                        independents_max = NULL; independents_max = xNew<IssmDouble>(M*N);
    192189                        for(int m=0;m<M;m++){
    … …  
    205202                        xDelete<IssmDouble>(independent);
    206203                        xDelete<IssmDouble>(independents_min);
    207204                        xDelete<IssmDouble>(independents_max);
    208        
     205
    209206                }
    210207                else{
    211208                        _error_("Independent cannot be of size " << types[i]);

  • ../trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp

     
    6868                        /*Intermediaries*/
    6969                        int            num_independent_objects,M;
    7070                        char**         names                   = NULL;
    71                                
     71
    7272                                /*this is done somewhere else*/
    7373                                parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_independent_objects",InversionNumControlParametersEnum));
    7474                           parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.num_dependent_objects",InversionNumCostFunctionsEnum));
    75                                
     75
    7676                                /*Step1: create controls (independents)*/
    7777                                iomodel->FetchData(&num_independent_objects,"md.autodiff.num_independent_objects");
    7878                                _assert_(num_independent_objects>0);
    … …  
    9595                                }
    9696                                xDelete<char*>(cm_responses);
    9797                                parameters->AddObject(new IntVecParam(InversionCostFunctionsEnum,costfunc_enums,num_costfunc));
    98                                
     98
    9999                                iomodel->FetchData(&control_scaling_factors,NULL,NULL,"md.autodiff.independent_scaling_factors");
    100100                                parameters->AddObject(new DoubleVecParam(InversionControlScalingFactorsEnum,control_scaling_factors,num_independent_objects));
    101        
     101
    102102                                /*cleanup*/
    103103                                for(int i=0;i<num_independent_objects;i++){
    104104                                        xDelete<char>(names[i]);

  • ../trunk-jpl/src/c/modules/ModelProcessorx/CreateOutputDefinitions.cpp

     
    1010void CreateOutputDefinitions(Elements* elements, Parameters* parameters,IoModel* iomodel){
    1111
    1212        int i,j;
    13        
     13
    1414        DataSet*     output_definitions      = NULL;
    1515        int*         output_definition_enums = NULL;
    1616        int          num_output_definitions;
    … …  
    6464                        }
    6565                        else if (output_definition_enums[i]==MisfitEnum){
    6666                                /*Deal with misfits: {{{*/
    67                        
     67
    6868                                /*misfit variables: */
    6969                                int          nummisfits;
    7070                                char**       misfit_name_s                                              = NULL;   
    … …  
    114114                                        else
    115115                                         _error_("misfit weight size not supported yet");
    116116
    117 
    118117                                        /*First create a misfit object for that specific string (misfit_model_string_s[j]):*/
    119118                                        output_definitions->AddObject(new Misfit(misfit_name_s[j],StringToEnumx(misfit_definitionstring_s[j]),StringToEnumx(misfit_model_string_s[j]),StringToEnumx(misfit_observation_string_s[j]),misfit_timeinterpolation_s[j],misfit_local_s[j],StringToEnumx(misfit_weights_string_s[j])));
    120119
    … …  
    158157                        }
    159158                        else if (output_definition_enums[i]==CfsurfacesquareEnum){
    160159                                /*Deal with cfsurfacesquare: {{{*/
    161                                
     160
    162161                                /*cfsurfacesquare variables: */
    163162                                int          num_cfsurfacesquares;
    164163                                char**       cfsurfacesquare_name_s                                             = NULL;   
    … …  
    213212                                        for(int k=0;k<elements->Size();k++){
    214213
    215214                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
    216                                                
     215
    217216                                                element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_observation_s[j], iomodel,cfsurfacesquare_observation_M_s[j],cfsurfacesquare_observation_N_s[j],obs_vector_type,StringToEnumx(cfsurfacesquare_observation_string_s[j]),7,SurfaceObservationEnum);
    218217                                                element->DatasetInputAdd(StringToEnumx(cfsurfacesquare_definitionstring_s[j]),cfsurfacesquare_weights_s[j], iomodel,cfsurfacesquare_weights_M_s[j],cfsurfacesquare_weights_N_s[j],weight_vector_type,StringToEnumx(cfsurfacesquare_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
    219218
    … …  
    250249                        }
    251250                        else if (output_definition_enums[i]==CfdragcoeffabsgradEnum){
    252251                                /*Deal with cfdragcoeffabsgrad: {{{*/
    253                                
     252
    254253                                /*cfdragcoeffabsgrad variables: */
    255254                                int          num_cfdragcoeffabsgrads;
    256255                                char**       cfdragcoeffabsgrad_name_s                                          = NULL;   
    … …  
    285284                                        for(int k=0;k<elements->Size();k++){
    286285
    287286                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
    288                                                
     287
    289288                                                element->DatasetInputAdd(StringToEnumx(cfdragcoeffabsgrad_definitionstring_s[j]),cfdragcoeffabsgrad_weights_s[j], iomodel,cfdragcoeffabsgrad_weights_M_s[j],cfdragcoeffabsgrad_weights_N_s[j],weight_vector_type,StringToEnumx(cfdragcoeffabsgrad_weights_string_s[j]),7,WeightsSurfaceObservationEnum);
    290289
    291290                                        }
    … …  
    312311                        }
    313312                        else if (output_definition_enums[i]==CfsurfacelogvelEnum){
    314313                                /*Deal with cfsurfacelogvel: {{{*/
    315                                
     314
    316315                                /*cfsurfacelogvel variables: */
    317316                                int          num_cfsurfacelogvels;
    318317                                char**       cfsurfacelogvel_name                                               = NULL;   
    … …  
    368367                                        for(int k=0;k<elements->Size();k++){
    369368
    370369                                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(k));
    371                                                
     370
    372371                                                element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vxobs[j], iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vxobs_string[j]),7,VxObsEnum);
    373372                                                        element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_vyobs[j], iomodel,cfsurfacelogvel_observation_M[j],cfsurfacelogvel_observation_N[j],obs_vector_type,StringToEnumx(cfsurfacelogvel_vyobs_string[j]),7,VyObsEnum);
    374373                                                element->DatasetInputAdd(StringToEnumx(cfsurfacelogvel_definitionstring[j]),cfsurfacelogvel_weights[j], iomodel,cfsurfacelogvel_weights_M[j],cfsurfacelogvel_weights_N[j],weight_vector_type,StringToEnumx(cfsurfacelogvel_weightstring[j]),7,WeightsSurfaceObservationEnum);
    … …  
    408407                        }
    409408                        else if (output_definition_enums[i]==NodalvalueEnum){
    410409                                /*Deal with nodal values: {{{*/
    411                                
     410
    412411                                /*nodal value variables: */
    413412                                int          numnodalvalues;
    414413                                char**       nodalvalue_name_s             = NULL;   
    … …  
    427426                                        /*First create a nodalvalue object for that specific enum (nodalvalue_model_enum_s[j]):*/
    428427                                        output_definitions->AddObject(new Nodalvalue(nodalvalue_name_s[j],StringToEnumx(nodalvalue_definitionstrings[j]),StringToEnumx(nodalvalue_modelstrings[j]),nodalvalue_node_s[j]-1)); //-1 because matlab to c indexing.
    429428                                }
    430                                        
     429
    431430                                /*Free ressources:*/
    432431                                for(j=0;j<numnodalvalues;j++){
    433432                                        char* string=NULL;
    … …  
    481480                        }
    482481                        else if (output_definition_enums[i]==MassconaxpbyEnum){
    483482                                /*Deal with masscon combinations: {{{*/
    484                                
     483
    485484                                /*masscon variables: */
    486485                                char**       masscon_name_s             = NULL;   
    487486                                char**           masscon_definitionstring_s             = NULL;   
    … …  
    616615                                        }
    617616                                        output_definitions->AddObject(new Numberedcostfunction(ncf_name_s[j],StringToEnumx(ncf_definitionstring_s[j]),num_cost_functions,cost_function_enums));
    618617                                }
    619                                
     618
    620619                                /*Free data: */
    621620                                iomodel->DeleteData(2,"md.numberedcostfunction.name","md.numberedcostfunction.definitionstring");
    622621                                xDelete<int>(cost_function_enums);

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

     
    2222        char**      requestedoutputs = NULL;
    2323        char*       fieldname = NULL;
    2424        IssmDouble  time;
    25        
     25
    2626        /*parameters for mass flux:*/
    2727        int          mass_flux_num_profiles     = 0;
    2828        bool         qmu_mass_flux_present      = false;
    … …  
    6161        parameters->AddObject(iomodel->CopyConstantObject("md.inversion.type",InversionTypeEnum));
    6262        parameters->AddObject(iomodel->CopyConstantObject("md.calving.law",CalvingLawEnum));
    6363        parameters->AddObject(new IntParam(SealevelriseRunCountEnum,1)); 
    64        
    6564
    6665          {/*This is specific to ice...*/
    6766                parameters->AddObject(iomodel->CopyConstantObject("md.mesh.elementtype",MeshElementtypeEnum));

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

     
    5454                analysis->UpdateElements(elements,iomodel,i,analysis_enum);
    5555                delete analysis;
    5656
    57 
    5857                /* Update counters, because we have created more nodes, loads and
    5958                 * constraints, and ids for objects created in next call to CreateDataSets
    6059                 * will need to start at the end of the updated counters: */

  • ../trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

     
    3434
    3535        /*Create elements*/
    3636        if(control_analysis && !adolc_analysis)iomodel->FetchData(2,"md.inversion.min_parameters","md.inversion.max_parameters");
    37        
     37
    3838        switch(iomodel->meshelementtype){
    3939                case TriaEnum:
    4040                        for(i=0;i<iomodel->numberofelements;i++){
    … …  
    122122                        }
    123123                        break;
    124124                case MaterialsEnum:
    125        
     125
    126126                        //we have several types of materials. Retrieve this info first:
    127127                        iomodel->FetchData(&nature,&nnat,&dummy,"md.materials.nature");
    128128
    … …  
    234234        if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
    235235        else iomodel->FetchDataToInput(elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.);
    236236        if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long");
    237        
     237
    238238        CreateNumberNodeToElementConnectivity(iomodel,solution_type);
    239239
    240240        int lid = 0;

  • ../trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp

     
    3131                offset += M[i]*N[i];
    3232                }
    3333
    34 
    3534                xDelete<int>(control_type);
    3635        }
    3736        else{

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

     
    6666        }
    6767}
    6868/*}}}*/
    69 

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

     
    66#include "../../shared/shared.h"
    77#include "../../toolkits/toolkits.h"
    88
    9 
    109void SmbGradientsx(FemModel* femmodel){/*{{{*/
    1110
    1211        // void SurfaceMassBalancex(hd,agd,ni){
    … …  
    317316                                correct for number of seconds in a year [s/yr]*/
    318317                        smb = smb/yts + anomaly;
    319318
    320 
    321319                        /*Update array accordingly*/
    322320                        smblist[v] = smb;
    323321

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

     
    106106        /*Free ressouces:*/
    107107        xDelete(dzT);
    108108        xDelete(dzB);
    109        
     109
    110110        //---------NEED TO IMPLEMENT A PROPER GRID STRECHING ALGORITHM------------
    111111
    112112        /*assign ouput pointers: */
    … …  
    201201        re=*pre;
    202202        gdn=*pgdn;
    203203        gsp=*pgsp;
    204        
     204
    205205        /*only when aIdx = 1 or 2 do we run grainGrowth: */
    206206        if(aIdx!=1 & aIdx!=2){
    207207                /*come out as we came in: */
    … …  
    220220        //set maximum water content by mass to 9 percent (Brun, 1980)
    221221        for(int i=0;i<m;i++)if(lwc[i]>9.0-Wtol) lwc[i]=9.0;
    222222
    223 
    224223        /* Calculate temperature gradiant across grid cells
    225224         * Returns the avereage gradient across the upper and lower grid cell */
    226225
    … …  
    243242
    244243        // take absolute value of temperature gradient
    245244        for(int i=0;i<m;i++)dT[i]=fabs(dT[i]);
    246        
     245
    247246        /*Snow metamorphism. Depends on value of snow dendricity and wetness of the snowpack: */
    248247        for(int i=0;i<m;i++){
    249248                if (gdn[i]>0.0+Gdntol){
    … …  
    323322                //convert grain size back to effective grain radius:
    324323                re[i]=gsz[i]/2.0;
    325324        }
    326        
     325
    327326        /*Free ressources:*/
    328327        xDelete<IssmDouble>(gsz);
    329328        xDelete<IssmDouble>(dT);
    … …  
    347346
    348347        //// Inputs
    349348        // aIdx      = albedo method to use
    350        
     349
    351350        // Method 0
    352351        //  aValue   = direct input value for albedo, override all changes to albedo
    353352
    … …  
    396395
    397396        /*Recover pointers: */
    398397        a=*pa;
    399        
     398
    400399        //some constants:
    401400        const IssmDouble dSnow = 300;   // density of fresh snow [kg m-3]       
    402401
    … …  
    512511void thermo(IssmDouble* pEC, IssmDouble** pT, IssmDouble* dz, IssmDouble* d, IssmDouble* swf, IssmDouble dlwrf, IssmDouble Ta, IssmDouble V, IssmDouble eAir, IssmDouble pAir, IssmDouble teValue, IssmDouble Ws, IssmDouble dt0, int m, IssmDouble Vz, IssmDouble Tz, IssmDouble thermo_scaling, IssmDouble dIce, int sid) { /*{{{*/
    513512
    514513        /* ENGLACIAL THERMODYNAMICS*/
    515          
     514
    516515        /* Description:
    517516           computes new temperature profile accounting for energy absorption and
    518517           thermal diffusion.*/
    … …  
    536535        // OUTPUTS
    537536        // T: grid cell temperature [k]
    538537        // EC: evaporation/condensation [kg]
    539        
     538
    540539        /*intermediary: */
    541540        IssmDouble* K = NULL;
    542541        IssmDouble* KU = NULL;
    … …  
    579578        IssmDouble dT_ulw=0.0;
    580579        IssmDouble dlw=0.0;
    581580        IssmDouble dT_dlw=0.0;
    582        
     581
    583582        /*outputs:*/
    584583        IssmDouble EC=0.0;
    585584        IssmDouble* T=*pT;
    … …  
    596595
    597596        //initialize Evaporation - Condenstation
    598597        EC = 0.0;
    599        
     598
    600599        // check if all SW applied to surface or distributed throught subsurface
    601600        // swIdx = length(swf) > 1
    602601
    … …  
    608607
    609608        // if V = 0 goes to infinity therfore if V = 0 change
    610609        if(V<0.01-Dtol)V=0.01;
    611        
     610
    612611        // Bulk-transfer coefficient for turbulent fluxes
    613612        An =  pow(0.4,2) / pow(log(Tz/z0),2);     // Bulk-transfer coefficient
    614613        C = An * dAir * V;                        // shf & lhf common coefficient
    … …  
    626625        for(int i=0;i<m;i++) if(d[i]>=dIce-Dtol) K[i] = 9.828 * exp(-5.7E-3*T[i]);
    627626
    628627        // THERMAL DIFFUSION COEFFICIENTS
    629  
     628
    630629        // A discretization scheme which truncates the Taylor-Series expansion
    631630        // after the 3rd term is used. See Patankar 1980, Ch. 3&4
    632  
     631
    633632        // discretized heat equation:
    634  
     633
    635634        //                 Tp = (Au*Tu° + Ad*Td° + (Ap-Au-Ad)Tp° + S) / Ap
    636  
     635
    637636        // where neighbor coefficients Au, Ap, & Ad are
    638  
     637
    639638        //                   Au = [dz_u/2KP + dz_p/2KE]^-1
    640639        //                   Ad = [dz_d/2KP + dz_d/2KD]^-1
    641640        //                   Ap = d*CI*dz/Dt
    … …  
    648647        KU = xNew<IssmDouble>(m);
    649648        KD = xNew<IssmDouble>(m);
    650649        KP = xNew<IssmDouble>(m);
    651        
     650
    652651        KU[0] = UNDEF;
    653652        KD[m-1] = UNDEF;
    654653        for(int i=1;i<m;i++) KU[i]= K[i-1];
    … …  
    660659        dzD = xNew<IssmDouble>(m);
    661660        dzU[0]=UNDEF;
    662661        dzD[m-1]=UNDEF;
    663        
     662
    664663        for(int i=1;i<m;i++) dzU[i]= dz[i-1];
    665664        for(int i=0;i<m-1;i++) dzD[i] = dz[i+1];
    666665
    … …  
    692691                Ad[i] = pow((dzD[i]/2.0/KP[i] + dz[i]/2.0/KD[i]),-1.0);
    693692                Ap[i] = (d[i]*dz[i]*CI)/dt;
    694693        }
    695        
     694
    696695        // create "neighbor" coefficient matrix
    697696        Nu = xNew<IssmDouble>(m);
    698697        Nd = xNew<IssmDouble>(m);
    … …  
    702701                Nd[i] = Ad[i] / Ap[i];
    703702                Np[i]= 1.0 - Nu[i] - Nd[i];
    704703        }
    705        
     704
    706705        // specify boundary conditions: constant flux at bottom
    707706        Nu[m-1] = 0.0;
    708707        Np[m-1] = 1.0;
    709        
     708
    710709        // zero flux at surface
    711710        Np[0] = 1.0 - Nd[0];
    712        
     711
    713712        // Create neighbor arrays for diffusion calculations instead of a tridiagonal matrix
    714713        Nu[0] = 0.0;
    715714        Nd[m-1] = 0.0;
    716        
     715
    717716        /* RADIATIVE FLUXES*/
    718717
    719718        // energy supplied by shortwave radiation [J]
    720719        sw = xNew<IssmDouble>(m);
    721720        for(int i=0;i<m;i++) sw[i]= swf[i] * dt;
    722        
     721
    723722        // temperature change due to SW
    724723        dT_sw = xNew<IssmDouble>(m);
    725724        for(int i=0;i<m;i++) dT_sw[i]= sw[i] / (CI * d[i] * dz[i]);
    … …  
    745744                // PART OF ENERGY CONSERVATION CHECK
    746745                // store initial temperature
    747746                //T_init = T;
    748    
     747
    749748                // calculate temperature of snow surface (Ts)
    750749                // when incoming SW radition is allowed to penetrate the surface,
    751750                // the modeled energy balance becomes very sensitive to how Ts is
    … …  
    753752                // less when Ts is taken as the mean of the x top grid cells.
    754753                Ts = (T[0] + T[1])/2.0;
    755754                Ts = fmin(CtoK,Ts);    // don't allow Ts to exceed 273.15 K (0 degC)
    756                
     755
    757756                //TURBULENT HEAT FLUX
    758    
     757
    759758                // Monin-Obukhov Stability Correction
    760759                // Reference:
    761760                // Ohmura, A., 1982: Climate and Energy-Balance on the Arctic Tundra.
    … …  
    763762
    764763                // calculate the Bulk Richardson Number (Ri)
    765764                Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
    766                
     765
    767766                // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
    768    
     767
    769768                // do not allow Ri to exceed 0.19
    770769                Ri = fmin(Ri, 0.19);
    771770
    … …  
    777776                else {
    778777                        coefM =pow (1.0-18*Ri,-0.25);
    779778                }
    780                
     779
    781780                // calculate heat/wind 'coef_H' stability factor
    782781                if (Ri <= -0.03+Ttol) coefH = coefM/1.3;
    783782                else coefH = coefM;
    784                
     783
    785784                //// Sensible Heat
    786785                // calculate the sensible heat flux [W m-2](Patterson, 1998)
    787786                shf = C * CA * (Ta - Ts);
    … …  
    800799                }
    801800                else{
    802801                        L = LS; // latent heat of sublimation
    803                        
     802
    804803                        // for an ice surface Murphy and Koop, 2005 [Equation 7]
    805804                        eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
    806805                }
    … …  
    822821                // upward longwave contribution
    823822                ulw = - (SB * pow(Ts,4.0)* teValue) * dt ;
    824823                dT_ulw = ulw / TCs;
    825                
     824
    826825                // new grid point temperature
    827    
     826
    828827                //SW penetrates surface
    829828                for(int j=0;j<m;j++) T[j] = T[j] + dT_sw[j];
    830829                T[0] = T[0] + dT_dlw + dT_ulw + dT_turb;
    831                
     830
    832831                // temperature diffusion
    833832                for(int j=0;j<m;j++) T0[1+j]=T[j];
    834833                for(int j=0;j<m;j++) Tu[j] = T0[j];
    835834                for(int j=0;j<m;j++) Td[j] = T0[2+j];
    836835                for(int j=0;j<m;j++) T[j] = (Np[j] * T[j]) + (Nu[j] * Tu[j]) + (Nd[j] * Td[j]);
    837                
     836
    838837                // calculate cumulative evaporation (+)/condensation(-)
    839838                EC = EC + (EC_day/dts)*dt;
    840839
    … …  
    872871        xDelete<IssmDouble>(Tu);
    873872        xDelete<IssmDouble>(Td);
    874873
    875 
    876874        /*Assign output pointers:*/
    877875        *pEC=EC;
    878876        *pT=T;
    … …  
    900898        // Outputs
    901899        //   swf     = absorbed shortwave radiation [W m-2]
    902900        //
    903        
     901
    904902        /*outputs: */
    905903        IssmDouble* swf=NULL;
    906904
    … …  
    909907        /*Initialize and allocate: */
    910908        swf=xNewZeroInit<IssmDouble>(m);
    911909
    912 
    913910        // SHORTWAVE FUNCTION
    914911        if (swIdx == 0) {// all sw radation is absorbed in by the top grid cell
    915        
     912
    916913                // calculate surface shortwave radiation fluxes [W m-2]
    917914                swf[0] = (1.0 - as) * dsw;
    918915        }
    … …  
    947944                                B2_cum[i]=1.0;
    948945                        }
    949946
    950 
    951947                        // spectral albedos:
    952948                        // 0.3 - 0.8um
    953949                        IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5));
    … …  
    988984                                B2_cum[i+1]=cum2;
    989985                        }
    990986
    991 
    992987                        // flux across grid cell boundaries
    993988                        Qs1 = xNew<IssmDouble>(m+1);
    994989                        Qs2 = xNew<IssmDouble>(m+1);
    … …  
    10141009                        xDelete<IssmDouble>(exp2);
    10151010                        xDelete<IssmDouble>(Qs1);
    10161011                        xDelete<IssmDouble>(Qs2);
    1017                        
    1018                        
     1012
    10191013                }
    10201014                else{  //function of grid cell density
    10211015
    … …  
    11431137        for(int i=0;i<m;i++)massinit+=mInit[i];
    11441138
    11451139        if (P > 0.0+Dtol){
    1146                        
    11471140
    11481141                if (T_air <= CtoK+Ttol){ // if snow
    11491142
    … …  
    12091202                mass=0; for(int i=0;i<m;i++)mass+=d[i]*dz[i];
    12101203
    12111204                mass_diff = mass - massinit - P;
    1212                
     1205
    12131206                #ifndef _HAVE_ADOLC_  //avoid round operation. only check in forward mode.
    12141207                mass_diff = round(mass_diff * 100.0)/100.0;
    12151208                if (mass_diff > 0) _error_("mass not conserved in accumulation function");
    … …  
    12521245        IssmDouble* EW=NULL;
    12531246        IssmDouble* M=NULL;
    12541247        int*       D=NULL;
    1255        
     1248
    12561249        IssmDouble sumM=0.0;
    12571250        IssmDouble sumER=0.0;
    12581251        IssmDouble addE=0.0;
    … …  
    12641257        IssmDouble dm=0.0;
    12651258        IssmDouble X=0.0;
    12661259        IssmDouble Wi=0.0;
    1267    
     1260
    12681261        IssmDouble Ztot=0.0;
    12691262        IssmDouble T_bot=0.0;
    12701263        IssmDouble m_bot=0.0;
    … …  
    12781271        IssmDouble EI_bot=0.0;
    12791272        IssmDouble EW_bot=0.0;
    12801273        bool        top=false;
    1281    
     1274
    12821275        IssmDouble* Zcum=NULL;
    12831276        IssmDouble* dzMin2=NULL;
    12841277        IssmDouble zY2=1.025;
    … …  
    12851278        bool lastCellFlag = false;
    12861279        int X1=0;
    12871280        int X2=0;
    1288    
     1281
    12891282        int        D_size = 0;
    12901283
    12911284        /*outputs:*/
    … …  
    13021295        IssmDouble* gsp=*pgsp;
    13031296        int         n=*pn;
    13041297        IssmDouble* R=NULL;
    1305        
     1298
    13061299        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   melt module\n");
    13071300
    13081301        //// INITIALIZATION
    … …  
    13341327        // new grid point center temperature, T [K]
    13351328        // for(int i=0;i<n;i++) T[i]-=exsT[i];
    13361329        for(int i=0;i<n;i++) T[i]=fmin(T[i],CtoK);
    1337        
     1330
    13381331        // specify irreducible water content saturation [fraction]
    13391332        const IssmDouble Swi = 0.07;                     // assumed constant after Colbeck, 1974
    13401333
    … …  
    14921485                                flxDn[i+1] = inM - F1 - dW[i];     // meltwater out       
    14931486                                T[i] = T[i] + ((F1+F2)*(LF+(CtoK - T[i])*CI)/(m[i]*CI));// change in temperature
    14941487
    1495 
    14961488                                // check if an ice layer forms
    14971489                                if (fabs(d[i] - dIce) < Dtol){
    14981490                                        // _printf_("ICE LAYER FORMS");
    … …  
    15041496                        }
    15051497                }
    15061498
    1507 
    15081499                //// GRID CELL SPACING AND MODEL DEPTH
    15091500                for(int i=0;i<n;i++)if (W[i] < 0.0 - Wtol) _error_("negative pore water generated in melt equations");
    1510                
     1501
    15111502                // delete all cells with zero mass
    15121503                // adjust pore water
    15131504                for(int i=0;i<n;i++)W[i] += dW[i];
    … …  
    15321523                celldelete(&EW,n,D,D_size);
    15331524                n=D_size;
    15341525                xDelete<int>(D);
    1535        
     1526
    15361527                // calculate new grid lengths
    15371528                for(int i=0;i<n;i++)dz[i] = m[i] / d[i];
    15381529
    … …  
    15441535        //Merging of cells as they are burried under snow.
    15451536        Zcum=xNew<IssmDouble>(n);
    15461537        dzMin2=xNew<IssmDouble>(n);
    1547    
     1538
    15481539        Zcum[0]=dz[0]; // Compute a cumulative depth vector
    1549    
     1540
    15501541        for (int i=1;i<n;i++){
    15511542                Zcum[i]=Zcum[i-1]+dz[i];
    15521543        }
    1553    
     1544
    15541545        for (int i=0;i<n;i++){
    15551546                if (Zcum[i]<=zTop+Dtol){
    15561547                        dzMin2[i]=dzMin;
    … …  
    15681559                        break;
    15691560                }
    15701561        }
    1571    
     1562
    15721563        //Last cell has to be treated separately if has to be merged (no underlying cell so need to merge with overlying cell)
    15731564        if(X==n-1){
    15741565                lastCellFlag = true;
    … …  
    15881579                        re[i+1] = (re[i]*m[i] + re[i+1]*m[i+1]) / m_new;
    15891580                        gdn[i+1] = (gdn[i]*m[i] + gdn[i+1]*m[i+1]) / m_new;
    15901581                        gsp[i+1] = (gsp[i]*m[i] + gsp[i+1]*m[i+1]) / m_new;
    1591            
     1582
    15921583                        // merge with underlying grid cell and delete old cell
    15931584                        dz[i+1] = dz[i] + dz[i+1];                 // combine cell depths
    15941585                        d[i+1] = m_new / dz[i+1];                   // combine top densities
    15951586                        W[i+1] = W[i+1] + W[i];                     // combine liquid water
    15961587                        m[i+1] = m_new;                             // combine top masses
    1597            
     1588
    15981589                        // set cell to 99999 for deletion
    15991590                        m[i] = -99999;
    16001591                }
    … …  
    16101601                                break;
    16111602                        }
    16121603                }
    1613        
     1604
    16141605                // adjust variables as a linearly weighted function of mass
    16151606                IssmDouble m_new = m[X2] + m[X1];
    16161607                T[X1] = (T[X2]*m[X2] + T[X1]*m[X1]) / m_new;
    … …  
    16181609                re[X1] = (re[X2]*m[X2] + re[X1]*m[X1]) / m_new;
    16191610                gdn[X1] = (gdn[X2]*m[X2] + gdn[X1]*m[X1]) / m_new;
    16201611                gsp[X1] = (gsp[X2]*m[X2] + gsp[X1]*m[X1]) / m_new;
    1621        
     1612
    16221613                // merge with underlying grid cell and delete old cell
    16231614                dz[X1] = dz[X2] + dz[X1];                 // combine cell depths
    16241615                d[X1] = m_new / dz[X1];                   // combine top densities
    16251616                W[X1] = W[X1] + W[X2];                     // combine liquid water
    16261617                m[X1] = m_new;                             // combine top masses
    1627        
     1618
    16281619                // set cell to 99999 for deletion
    16291620                m[X2] = -99999;
    16301621        }
    … …  
    16471638        celldelete(&EW,n,D,D_size);
    16481639        n=D_size;
    16491640        xDelete<int>(D);
    1650    
     1641
    16511642        // check if any of the top 10 cell depths are too large
    16521643        X=0;
    16531644        for(int i=9;i>=0;i--){
    … …  
    16561647                        break;
    16571648                }
    16581649        }
    1659        
     1650
    16601651        int j=0;
    16611652        while(j<=X){
    16621653                if (dz[j] > dzMin*2.0+Dtol){
    … …  
    16851676        // INTERPRETATION
    16861677        // calculate total model depth
    16871678        Ztot = cellsum(dz,n);
    1688    
     1679
    16891680        if (Ztot < zMin-Dtol){
    16901681                // printf("Total depth < zMin %f \n", Ztot);
    16911682                // mass and energy to be added
    16921683                mAdd = m[n-1]+W[n-1];
    16931684                addE = T[n-1]*m[n-1]*CI + W[n-1]*(LF+CtoK*CI);
    1694        
     1685
    16951686                // add a grid cell of the same size and temperature to the bottom
    16961687                dz_bot=dz[n-1];
    16971688                T_bot=T[n-1];
    … …  
    17041695                gsp_bot=gsp[n-1];
    17051696                EI_bot=EI[n-1];
    17061697                EW_bot=EW[n-1];
    1707        
     1698
    17081699                dz_add=dz_bot;
    1709        
     1700
    17101701                newcell(&dz,dz_bot,top,n);
    17111702                newcell(&T,T_bot,top,n);
    17121703                newcell(&W,W_bot,top,n);
    … …  
    17261717                mAdd = -(m[n-1]+W[n-1]);
    17271718                addE = -(T[n-1]*m[n-1]*CI) - (W[n-1]*(LF+CtoK*CI));
    17281719                dz_add=-(dz[n-1]);
    1729        
     1720
    17301721                // remove a grid cell from the bottom
    17311722                D_size=n-1;
    17321723                D=xNew<int>(D_size);
    1733        
     1724
    17341725                for(int i=0;i<n-1;i++) D[i]=i;
    17351726                celldelete(&dz,n,D,D_size);
    17361727                celldelete(&T,n,D,D_size);
    … …  
    17671758        if (dm !=0  || dE !=0) _error_("mass or energy are not conserved in melt equations\n"
    17681759                        << "dm: " << dm << " dE: " << dE << "\n");
    17691760        #endif
    1770        
     1761
    17711762        /*Free ressources:*/
    17721763        if(m)xDelete<IssmDouble>(m);
    17731764        if(EI)xDelete<IssmDouble>(EI);
    … …  
    17831774        if(M)xDelete<IssmDouble>(M);
    17841775        xDelete<IssmDouble>(Zcum);
    17851776        xDelete<IssmDouble>(dzMin2);
    1786    
     1777
    17871778        /*Assign output pointers:*/
    17881779        *pM=sumM;
    17891780        *pR=Rsum;
    … …  
    18611852        /*output: */
    18621853        IssmDouble* dz=NULL;
    18631854        IssmDouble* d=NULL;
    1864        
     1855
    18651856        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   densification module\n");
    18661857
    18671858        /*Recover pointers: */
    … …  
    18701861
    18711862        // initial mass
    18721863        IssmDouble* mass_init = xNew<IssmDouble>(m);for(int i=0;i<m;i++) mass_init[i]=d[i] * dz[i];
    1873        
     1864
    18741865        /*allocations and initialization of overburden pressure and factor H: */
    18751866        IssmDouble* cumdz = xNew<IssmDouble>(m-1);
    18761867        cumdz[0]=dz[0];
    … …  
    18791870        IssmDouble* obp = xNew<IssmDouble>(m);
    18801871        obp[0]=0.0;
    18811872        for(int i=1;i<m;i++)obp[i]=cumdz[i-1]*d[i-1];
    1882        
     1873
    18831874        // calculate new snow/firn density for:
    18841875        //   snow with densities <= 550 [kg m-3]
    18851876        //   snow with densities > 550 [kg m-3]
    1886                
    1887        
     1877
    18881878        for(int i=0;i<m;i++){
    18891879                switch (denIdx){
    18901880                        case 1: // Herron and Langway (1980)
    … …  
    20031993        IssmDouble shf=0.0;
    20041994        IssmDouble lhf=0.0;
    20051995        IssmDouble EC=0.0;
    2006        
     1996
    20071997        if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_("   turbulentFlux module\n");
    20081998
    20091999        // calculated air density [kg/m3]
    … …  
    20282018        Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
    20292019
    20302020        // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
    2031        
     2021
    20322022        // do not allow Ri to exceed 0.19
    20332023        Ri = fmin(Ri, 0.19);
    20342024
    … …  
    20722062                eS = exp(9.550426 - 5723.265/Ts + 3.53068 * log(Ts) - 0.00728332 * Ts);
    20732063        }
    20742064
    2075 
    20762065        // Latent heat flux [W m-2]
    20772066        lhf = C * L * (eAir - eS) * 0.622 / pAir;
    20782067

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

     
    106106
    107107        femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
    108108        femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
    109        
     109
    110110        IssmDouble* boxareas=xNewZeroInit<IssmDouble>(num_basins*maxbox);
    111111
    112112        for(int i=0;i<femmodel->elements->Size();i++){
    … …  
    121121        IssmDouble* sumareas =xNew<IssmDouble>(num_basins*maxbox);
    122122        ISSM_MPI_Allreduce(boxareas,sumareas,num_basins*maxbox,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    123123        //if(sumareas[0]==0){_error_("No elements in box 0, basal meltrates will be 0. Consider decreasing md.basalforcings.maxboxcount or refining your mesh!");}
    124        
     124
    125125        /*Update parameters to keep track of the new areas in future calculations*/
    126126        femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoBoxAreaEnum,sumareas,maxbox*num_basins));
    127127

  • ../trunk-jpl/src/c/modules/ExpToLevelSetx/ExpToLevelSetxt.cpp

     
    5151        double x1,y1;
    5252        double x2,y2;
    5353        double mind;
    54        
     54
    5555        for(int i=i0;i<i1;i++){
    5656
    5757      /*Get current point*/
    … …  
    9292      /*Beyond the 'w' end of the segment*/
    9393      return sqrt(pow(x2-x0,2)+pow(y2-y0,2));
    9494   }
    95        
     95
    9696   /*Projection falls on segment*/
    9797        double projx= x1 + t* (x2-x1);
    9898        double projy= y1 + t* (y2-y1);

  • ../trunk-jpl/src/c/solutionsequences/solutionsequence_linear.cpp

     
    1919        Vector<IssmDouble>*  ys  = NULL;
    2020        int  configuration_type;
    2121        IssmDouble solver_residue_threshold;
    22        
     22
    2323        /*solver convergence test: */
    2424        IssmDouble nKUF;
    2525        IssmDouble nF;
    … …  
    3838        femmodel->profiler->Start(SOLVER);
    3939        Solverx(&uf, Kff, pf, NULL, df, femmodel->parameters);
    4040        femmodel->profiler->Stop(SOLVER);
    41        
     41
    4242        /*Check that the solver converged nicely: */
    43                
     43
    4444        //compute KUF = KU - F = K*U - F
    4545        KU=uf->Duplicate(); Kff->MatMult(uf,KU);
    4646        KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0);
    … …  
    5252
    5353        if(!xIsNan(solver_residue_threshold) && solver_residue>solver_residue_threshold)_error_("   solver residue too high!: norm(KU-F)/norm(F)=" << solver_residue << "\n");
    5454
    55 
    5655        //clean up
    5756        delete KU; delete KUF;
    5857        delete Kff; delete pf; delete df;
    … …  
    6261        //          ADOLC_DUMP_MACRO(uf->svector->vector[i]);
    6362        //        }
    6463        //#endif
    65        
     64
    6665        Mergesolutionfromftogx(&ug, uf,ys,femmodel->nodes,femmodel->parameters);delete uf; delete ys;
    6766        InputUpdateFromSolutionx(femmodel,ug);
    6867        delete ug; 

  • ../trunk-jpl/src/c/solutionsequences/solutionsequence_fct.cpp

     
    184184
    185185        /*Initial guess*/
    186186        VecZeroEntries(udot);
    187        
     187
    188188        /*Richardson's iterations*/
    189189        for(int i=0;i<5;i++){
    190190                MatMult(Mc,udot,temp1);
    … …  
    401401        /*Create RHS: [ML + (1 − theta) deltaT L^n] u^n */
    402402        CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
    403403        delete uf;
    404        
     404
    405405        /*Go solve lower order solution*/
    406406        femmodel->profiler->Start(SOLVER);
    407407        SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);

  • ../trunk-jpl/src/c/solutionsequences/solutionsequence_thermal_nonlinear.cpp

     
    6060        }
    6161
    6262        count=1;
    63        
     63
    6464        for(;;){
    6565                delete tf_old;tf_old=tf;
    6666

  • ../trunk-jpl/src/c/toolkits/mumps/MumpsSolve.cpp

     
    177177        rhs=xNew<IssmDouble>(pf_M);
    178178#endif
    179179
    180 
    181180        recvcounts=xNew<int>(num_procs);
    182181        displs=xNew<int>(num_procs);
    183182
    … …  
    262261        rhs=xNew<IssmDouble>(pf_M);
    263262#endif
    264263
    265 
    266264        recvcounts=xNew<int>(num_procs);
    267265        displs=xNew<int>(num_procs);
    268266

  • ../trunk-jpl/src/c/cores/controlad_core.cpp

     
    104104                case 7:  _printf0_("   <g,d> > 0  or  <y,s> <0\n"); break;
    105105                default: _printf0_("   Unknown end condition\n");
    106106        }
    107        
     107
    108108        /*Save results:*/
    109109        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,G,n,1,0,0.0));
    110110        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffXpEnum,X,intn,1,0,0.0));
    … …  
    129129        int           intn;
    130130        IssmPDouble*   X=NULL;
    131131        int            i;
    132        
     132
    133133        /*Now things get complicated. The femmodel we recovered did not initialize an AD trace, so we can't compute gradients with it. We are going to recreate
    134134         *a new femmodel, identical in all aspects to the first one, with trace on though, which will allow us to run the forward mode and get the gradient
    135135         in one run of the solution core. So first recover the filenames required for the FemModel constructor, then call a new ad tailored constructor:*/
    … …  
    142142
    143143        femmodel=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,NULL);
    144144
    145        
    146145        /*Get initial guess:*/
    147146        femmodel->parameters->FindParam(&Xd,&intn,AutodiffXpEnum);
    148147        X=xNew<IssmPDouble>(intn); for(i=0;i<intn;i++) X[i]=reCast<IssmPDouble>(Xd[i]);
    … …  
    174173        FemModel   *femmodel  =  NULL;
    175174        IssmDouble    pfd;
    176175        int            i;
    177        
     176
    178177        /*Recover Femmodel*/
    179178        femmodel  = (FemModel*)dzs;
    180179
    … …  
    194193        IssmDouble* Jlist = NULL;
    195194        femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
    196195        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    197        
     196
    198197        /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
    199198        adgradient_core(femmodel);
    200199
    … …  
    210209
    211210        /*MPI broadcast results:*/
    212211        ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    213        
     212
    214213        /*Send gradient to m1qn3 core: */
    215214        for(long i=0;i<*n;i++) G[i] = G2[i];
    216        
     215
    217216        /*Constrain X and G*/
    218217        IssmDouble  Gnorm = 0.;
    219218        for(long i=0;i<*n;i++) Gnorm += G[i]*G[i];
    … …  
    227226        /*Clean-up and return*/
    228227        xDelete<IssmDouble>(Jlist);
    229228        xDelete<IssmPDouble>(G2);
    230        
     229
    231230        xDelete<char>(rootpath);
    232231        xDelete<char>(inputfilename);
    233232        xDelete<char>(outputfilename);
    … …  
    250249        FemModel   *femmodelad  = NULL;
    251250        IssmDouble    pfd;
    252251        int            i;
    253        
     252
    254253        /*Recover Femmodel*/
    255254        femmodel  = (FemModel*)dzs;
    256255
    … …  
    270269
    271270        femmodelad=new FemModel(rootpath, inputfilename, outputfilename, toolkitsfilename, lockfilename, restartfilename,IssmComm::GetComm(), femmodel->solution_type,X);
    272271        femmodel=femmodelad; //We can do this, because femmodel is being called from outside, not by reference, so we won't erase it
    273        
     272
    274273        /*Recover some parameters*/
    275274        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    276275
    … …  
    283282        IssmDouble* Jlist = NULL;
    284283        femmodel->CostFunctionx(&pfd,&Jlist,NULL); *pf=reCast<IssmPDouble>(pfd);
    285284        _printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
    286        
     285
    287286        /*Compute gradient using AD. Gradient is in the results after the ad_core is called*/
    288287        adgradient_core(femmodel);
    289288
    … …  
    299298
    300299        /*MPI broadcast results:*/
    301300        ISSM_MPI_Bcast(G2,*n,ISSM_MPI_PDOUBLE,0,IssmComm::GetComm());
    302        
     301
    303302        /*Send gradient to m1qn3 core: */
    304303        for(long i=0;i<*n;i++) G[i] = G2[i];
    305        
     304
    306305        /*Recover Gnorm: */
    307306        IssmDouble  Gnorm = 0.;
    308307        for(int i=0;i<*n;i++) Gnorm += G[i]*G[i];
    … …  
    316315        /*Clean-up and return*/
    317316        xDelete<IssmDouble>(Jlist);
    318317        xDelete<IssmPDouble>(G2);
    319        
    320318        xDelete<char>(rootpath);
    321319        xDelete<char>(inputfilename);
    322320        xDelete<char>(outputfilename);

  • ../trunk-jpl/src/c/cores/gia_core.cpp

     
    5252                int outputs[2] = {GiaWEnum,GiadWdtEnum};
    5353                femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
    5454        }
    55        
     55
    5656        xDelete<IssmDouble>(x);
    5757        xDelete<IssmDouble>(y);
    5858}

  • ../trunk-jpl/src/c/cores/bmb_core.cpp

     
    1111
    1212void bmb_core(FemModel* femmodel){
    1313
    14 
    1514        /*First, get BMB model from parameters*/
    1615        int  basalforcing_model;
    1716        bool isplume = false;

  • ../trunk-jpl/src/c/cores/dakota_core.cpp

     
    234234/*}}}*/
    235235void dakota_core(FemModel* femmodel){  /*{{{*/
    236236
    237 
    238237        int                my_rank;
    239238        char              *dakota_input_file  = NULL;
    240239        char              *dakota_output_file = NULL;

  • ../trunk-jpl/src/c/cores/esa_core.cpp

     
    2222        int solution_type;
    2323        int        numoutputs        = 0;
    2424        char     **requested_outputs = NULL;
    25        
     25
    2626        /*additional parameters: */
    2727        int  gsize;
    2828        bool spherical=true;
    … …  
    7979                U_east = new Vector<IssmDouble>(gsize);
    8080                U_x = new Vector<IssmDouble>(gsize);
    8181                U_y = new Vector<IssmDouble>(gsize);
    82                
     82
    8383                /*call the geodetic main modlule:*/
    8484                if(domaintype==Domain3DsurfaceEnum){
    8585                        femmodel->EsaGeodetic3D(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz);
    … …  
    9494                InputUpdateFromVectorx(femmodel,U_radial,EsaUmotionEnum,VertexSIdEnum); // radial displacement
    9595                InputUpdateFromVectorx(femmodel,U_north,EsaNmotionEnum,VertexSIdEnum);  // north motion
    9696                InputUpdateFromVectorx(femmodel,U_east,EsaEmotionEnum,VertexSIdEnum);           // east motion
    97                
     97
    9898                if(save_results){
    9999                        if(VerboseSolution()) _printf0_("   saving results\n");
    100100                        femmodel->parameters->FindParam(&requested_outputs,&numoutputs,EsaRequestedOutputsEnum);
    101101                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    102102                }
    103                        
     103
    104104                if(solution_type==EsaSolutionEnum)femmodel->RequestedDependentsx();
    105105
    106106                /*Free ressources:*/   
    … …  
    114114
    115115}
    116116/*}}}*/
    117 

  • ../trunk-jpl/src/c/cores/masstransport_core.cpp

     
    3030        femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
    3131        femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
    3232        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
    33                        
     33
    3434        if(VerboseSolution()) _printf0_("   computing mass transport\n");
    3535
    3636        /*Transport mass or free surface*/

  • ../trunk-jpl/src/c/cores/sealevelrise_core.cpp

     
    1212/*cores:*/
    1313void sealevelrise_core(FemModel* femmodel){ /*{{{*/
    1414
    15 
    1615        /*Parameters, variables:*/
    1716        bool save_results;
    1817        bool isslr=0;
    1918        int solution_type;
    20                
     19
    2120        /*Retrieve parameters:*/
    2221        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
    2322        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    2423        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    25        
     24
    2625        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
    2726        if(solution_type==SealevelriseSolutionEnum)isslr=1;
    2827
    … …  
    4039
    4140        /*Run steric core for sure:*/
    4241        steric_core(femmodel);
    43        
     42
    4443        /*Save results: */
    4544        if(save_results){
    4645                int        numoutputs        = 0;
    … …  
    5857/*}}}*/
    5958void geodetic_core(FemModel* femmodel){ /*{{{*/
    6059
    61 
    6260        /*variables:*/
    6361        Vector<IssmDouble> *RSLg    = NULL;
    6462        Vector<IssmDouble> *RSLg_rate    = NULL;
    … …  
    8785
    8886        /*Should we even be here?:*/
    8987        femmodel->parameters->FindParam(&geodetic,SealevelriseGeodeticEnum); if(!geodetic)return;
    90        
     88
    9189        /*Verbose: */
    9290        if(VerboseSolution()) _printf0_("         computing geodetic sea level rise\n");
    9391
    … …  
    154152
    155153                /*recover N_esa  = U_esa + RSLg:*/
    156154                N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
    157                
     155
    158156                /*transform these values into rates (as we only run this once each frequency turn:*/
    159157                N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
    160158                U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
    … …  
    161159                N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
    162160                U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
    163161                RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
    164                
     162
    165163                /*get some results into elements:{{{*/
    166164                InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum);
    167165                InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum);
    … …  
    181179                } /*}}}*/
    182180        }
    183181
    184 
    185182        if(iscoupler){
    186183                /*transfer sea level back to ice caps:*/
    187184                TransferSealevel(femmodel,SealevelNEsaRateEnum);
    … …  
    227224        Vector<IssmDouble> *N_esa_rate= NULL;
    228225        Vector<IssmDouble> *U_gia_rate= NULL;
    229226        Vector<IssmDouble> *N_gia_rate= NULL;
    230                
     227
    231228        /*parameters: */
    232229        bool isslr=0;
    233230        int  solution_type;
    … …  
    242239
    243240        /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
    244241        if(solution_type==SealevelriseSolutionEnum)isslr=1;
    245        
     242
    246243        /*Should we be here?:*/
    247244        if(!isslr)return;
    248245
    … …  
    290287}
    291288/*}}}*/
    292289Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel){ /*{{{*/
    293  
     290
    294291        /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
    295292
    296293        Vector<IssmDouble> *RSLgi    = NULL;
    … …  
    317314
    318315        /*Figure out size of g-set deflection vector and allocate solution vector: */
    319316        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
    320        
     317
    321318        /*Initialize:*/
    322319        RSLgi = new Vector<IssmDouble>(gsize);
    323320
    … …  
    330327        /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the
    331328         * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
    332329        RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
    333        
     330
    334331        /*save eustatic value for results: */
    335332        femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic));
    336333
    337 
    338334        /*clean up and return:*/
    339335        xDelete<IssmDouble>(latitude);
    340336        xDelete<IssmDouble>(longitude);
    … …  
    346342        /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
    347343          non eustatic core of the SLR solution */
    348344
    349 
    350345        Vector<IssmDouble> *RSLg    = NULL;
    351346        Vector<IssmDouble> *RSLg_old    = NULL;
    352347
    … …  
    379374        femmodel->parameters->FindParam(&eps_abs,SealevelriseAbstolEnum);
    380375        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
    381376        femmodel->parameters->FindParam(&loop,SealevelriseLoopIncrementEnum);
    382        
     377
    383378        /*computational flag: */
    384379        femmodel->parameters->FindParam(&rotation,SealevelriseRotationEnum);
    385380
    … …  
    388383
    389384        /*Figure out size of g-set deflection vector and allocate solution vector: */
    390385        gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum);
    391        
     386
    392387        /*Initialize:*/
    393388        RSLg = new Vector<IssmDouble>(gsize);
    394389        RSLg->Assemble();
    … …  
    399394
    400395        count=1;
    401396        converged=false;
    402        
     397
    403398        /*Start loop: */
    404399        for(;;){
    405400
    … …  
    412407
    413408                /*call the non eustatic module: */
    414409                femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, latitude, longitude, radius,verboseconvolution,loop);
    415        
     410
    416411                /*assemble solution vector: */
    417412                RSLgo->Assemble();
    418413
    … …  
    422417                        RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
    423418                        femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz,latitude,longitude,radius);
    424419                        RSLgo_rot->Assemble();
    425                        
     420
    426421                        /*save changes in inertia tensor as results: */
    427422                        femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
    428423                        femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz));
    … …  
    430425
    431426                        RSLgo->AXPY(RSLgo_rot,1);
    432427                }
    433                
     428
    434429                /*we need to average RSLgo over the ocean: RHS term  5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
    435430                RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo);
    436        
     431
    437432                /*RSLg is the sum of the eustatic term, and the ocean terms: */
    438433                RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1);
    439434                RSLg->Shift(-RSLgo_oceanaverage);
    … …  
    454449                        converged=true;
    455450                        break;
    456451                }       
    457                
     452
    458453                /*some minor verbosing adjustment:*/
    459454                if(count>1)verboseconvolution=false;
    460                
     455
    461456        }
    462457        if(VerboseConvergence()) _printf0_("\n              total number of iterations: " << count-1 << "\n");
    463458
    … …  
    473468        Vector<IssmDouble> *U_esa  = NULL;
    474469        Vector<IssmDouble> *U_north_esa   = NULL;
    475470        Vector<IssmDouble> *U_east_esa    = NULL;
    476                
     471
    477472        /*parameters: */
    478473        int configuration_type;
    479474        int  gsize;
    … …  
    506501        /*retrieve geometric information: */
    507502        VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
    508503        VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
    509        
     504
    510505        /*call the elastic main modlule:*/
    511506        femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg,latitude,longitude,radius,xx,yy,zz,loop,horiz);
    512507
    … …  
    531526        /*variables:*/
    532527        Vector<IssmDouble> *U_gia  = NULL;
    533528        Vector<IssmDouble> *N_gia  = NULL;
    534        
     529
    535530        /*parameters:*/
    536531        int                                     frequency;
    537532        IssmDouble          dt;
    … …  
    569564        IssmDouble*  forcing=NULL;
    570565        Vector<IssmDouble>* forcingglobal=NULL;
    571566        int*         nvs=NULL;
    572        
     567
    573568        /*transition vectors:*/
    574569        IssmDouble** transitions=NULL;
    575570        int          ntransitions;
    … …  
    576571        int*         transitions_m=NULL;
    577572        int*         transitions_n=NULL;
    578573        int          nv;
    579        
     574
    580575        /*communicators:*/
    581576        ISSM_MPI_Comm tocomm;
    582577        ISSM_MPI_Comm* fromcomms=NULL;
    … …  
    590585        femmodel->parameters->FindParam(&earthid,EarthIdEnum);
    591586        femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
    592587        my_rank=IssmComm::GetRank();
    593        
     588
    594589        /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
    595590        if(modelid==earthid){
    596591                GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
    … …  
    619614                                forcings[i]=xNew<IssmDouble>(nvs[i]);
    620615                                ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
    621616                        }
    622                        
     617
    623618                }
    624619                else{
    625620                        ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
    … …  
    630625
    631626        /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
    632627        if(modelid==earthid){
    633                
     628
    634629                /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions.
    635630                 *First, build the global delta thickness vector in the earth model: */
    636631                nv=femmodel->vertices->NumberOfVertices();
    … …  
    652647                                /*build index to plug values: */
    653648                                int*        index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
    654649
    655 
    656650                                /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular
    657651                                 * ice cap: */
    658652                                forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
    … …  
    662656
    663657                /*Assemble vector:*/
    664658                forcingglobal->Assemble();
    665                
     659
    666660                /*Plug into elements:*/
    667661                InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
    668662        }
    … …  
    695689        /*forcing being transferred from earth to ice caps: */
    696690        IssmDouble*  forcing=NULL;
    697691        IssmDouble*  forcingglobal=NULL;
    698        
     692
    699693        /*transition vectors:*/
    700694        IssmDouble** transitions=NULL;
    701695        int          ntransitions;
    … …  
    702696        int*         transitions_m=NULL;
    703697        int*         transitions_n=NULL;
    704698        int          nv;
    705        
     699
    706700        /*communicators:*/
    707701        ISSM_MPI_Comm fromcomm;
    708702        ISSM_MPI_Comm* tocomms=NULL;
    … …  
    717711        femmodel->parameters->FindParam(&earthid,EarthIdEnum);
    718712        femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
    719713        my_rank=IssmComm::GetRank();
    720        
     714
    721715        /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
    722716        if(modelid==earthid){
    723717                GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
    … …  
    732726                //femmodel->parameters->FindParam((int*)(&fromcomm), IcecapToEarthCommEnum);
    733727        }
    734728
    735 
    736729        /*Retrieve sea-level on earth model: */
    737730        if(modelid==earthid){
    738731                nv=femmodel->vertices->NumberOfVertices();
    … …  
    741734
    742735        /*Send the forcing to the ice caps:{{{*/
    743736        if(my_rank==0){
    744                
     737
    745738                if(modelid==earthid){
    746                        
     739
    747740                        /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
    748741                        femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
    749                        
     742
    750743                        if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
    751744
    752745                        for(int i=0;i<earthid;i++){
    … …  
    803796        Vector<IssmDouble> *deltathickness    = NULL;
    804797        Vector<IssmDouble> *cumdeltathickness    = NULL;
    805798        int nv;
    806        
     799
    807800        if(VerboseSolution()) _printf0_("              computing earth mass transport\n");
    808801
    809802        /*This mass transport module for the Earth is because we might have thickness variations as spcs
    … …  
    839832
    840833} /*}}}*/
    841834void slrconvergence(bool* pconverged, Vector<IssmDouble>* RSLg,Vector<IssmDouble>* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs){ /*{{{*/
    842        
     835
    843836        bool converged=true;
    844837        IssmDouble ndS,nS;
    845838        Vector<IssmDouble> *dRSLg    = NULL;
    … …  
    847840        //compute norm(du) and norm(u) if requested
    848841        dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
    849842        ndS=dRSLg->Norm(NORM_TWO);
    850        
     843
    851844        if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
    852        
     845
    853846        if(!xIsNan<IssmDouble>(eps_rel)){
    854847                nS=RSLg_old->Norm(NORM_TWO);
    855848                if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
    856849        }
    857850
    858 
    859851        //clean up
    860852        delete dRSLg;
    861853

  • ../trunk-jpl/src/c/cores/stressbalance_core.cpp

     
    2121        int        numoutputs        = 0;
    2222        char     **requested_outputs = NULL;
    2323        Analysis  *analysis          = NULL;
    24                        
    2524
    2625        /* recover parameters:*/
    2726        femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
    … …  
    3534        femmodel->parameters->FindParam(&numoutputs,StressbalanceNumRequestedOutputsEnum);
    3635        femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
    3736        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum);
    38        
     37
    3938        if(VerboseSolution()) _printf0_("   computing new velocity\n");
    4039
    4140        /*Compute slopes if necessary */
    … …  
    7877                delete analysis;
    7978        }
    8079
    81 
    8280        if(save_results){
    8381                if(VerboseSolution()) _printf0_("   saving results\n");
    8482                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    8583        }
    86        
     84
    8785        if(solution_type==StressbalanceSolutionEnum && !control_analysis)femmodel->RequestedDependentsx();
    8886
    8987        /*Free ressources:*/   

  • ../trunk-jpl/src/c/cores/adgradient_core.cpp

     
    3838        /*intermediary: */
    3939        IssmPDouble *aWeightVector=NULL;
    4040        IssmPDouble *weightVectorTimesJac=NULL;
    41        
     41
    4242        /*output: */
    4343        IssmPDouble *totalgradient=NULL;
    4444
    … …  
    5252
    5353                        /*First, stop tracing: */
    5454                        trace_off();
    55                        
     55
    5656                        if(VerboseAutodiff()){ /*{{{*/
    5757                                size_t  tape_stats[15];
    5858                                tapestats(my_rank,tape_stats); //reading of tape statistics
    … …  
    9595                        }
    9696                        delete [] sstats;
    9797                } /*}}}*/
    98                
     98
    9999                        /*retrieve parameters: */
    100100                        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    101101                        femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
    102                        
     102
    103103                        /*if no dependents, no point in running a driver: */
    104104                        if(!(num_dependents*num_independents)) return;
    105105
    … …  
    108108                        if (my_rank!=0){
    109109                                num_dependents=0; num_independents=0;
    110110                        }
    111                        
     111
    112112                        /*retrieve state variable: */
    113113                        femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
    114                        
     114
    115115                        /* driver argument */
    116116                        xp=xNew<double>(num_independents);
    117117                        for(i=0;i<num_independents;i++){
    … …  
    130130
    131131                        /*Initialize outputs: */
    132132                        totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
    133                        
     133
    134134                        for(aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
    135135
    136136                                /*initialize direction index in the weights vector: */
    137137                                aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
    138138                                if (my_rank==0) aWeightVector[aDepIndex]=1.0;
    139                                
     139
    140140                                /*initialize output gradient: */
    141141                                weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
    142142
    … …  
    161161                                xDelete(weightVectorTimesJac);
    162162                                xDelete(aWeightVector);
    163163                        }
    164                
     164
    165165                        /*add totalgradient to results*/
    166166                        femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0));
    167167
    168168                        if(VerboseAutodiff())_printf0_("   end ad core\n");
    169                        
     169
    170170                        /* delete the allocated space for the parameters and free ressources:{{{*/
    171171                        xDelete(anEDF_for_solverx_p->dp_x);
    172172                        xDelete(anEDF_for_solverx_p->dp_X);

  • ../trunk-jpl/src/c/cores/controladm1qn3_core.cpp

     
    125125        ISSM_MPI_Bcast(aX,intn,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    126126        SetControlInputsFromVectorx(femmodel,aX);
    127127        xDelete<IssmDouble>(aX);
    128        
     128
    129129        /*Compute solution (forward)*/
    130130        void (*solutioncore)(FemModel*)=NULL;
    131131        CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
    … …  
    391391        /*Optimization criterions*/
    392392        long niter = long(maxsteps); /*Maximum number of iterations*/
    393393        long nsim  = long(maxiter);/*Maximum number of function calls*/
    394        
     394
    395395        /*Get initial guess*/
    396396        Vector<double> *Xpetsc = NULL;
    397397
    … …  
    400400        Xpetsc->GetSize(&intn);
    401401        delete Xpetsc;
    402402        //_assert_(intn==numberofvertices*num_controls);
    403  
     403
    404404        /*Get problem dimension and initialize gradient and initial guess*/
    405405        long n = long(intn);
    406406        G = xNew<double>(n);
    … …  
    414414                }
    415415                N_add+=N[c];
    416416        }
    417        
     417
    418418        /*Allocate m1qn3 working arrays (see documentation)*/
    419419        long      m   = 100;
    420420        long      ndz = 4*n+m*(2*n+1);
    … …  
    470470                }
    471471                N_add +=N[c];
    472472        }
    473        
     473
    474474        /*Set X as our new control*/
    475475        IssmDouble* aX=xNew<IssmDouble>(intn);
    476476        IssmDouble* aG=xNew<IssmDouble>(intn);
    … …  
    482482
    483483        ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
    484484        SetControlInputsFromVectorx(femmodel,aX);
    485        
     485
    486486        xDelete(aX);
    487487
    488 
    489488        if (solution_type == TransientSolutionEnum){
    490489                int step = 1;
    491490                femmodel->parameters->SetParam(step,StepEnum);
    … …  
    505504                        /*Add to results*/
    506505                        femmodel->results->AddObject(G_output);
    507506                        femmodel->results->AddObject(X_output);
    508                        
     507
    509508                        offset += N[i]*numberofvertices;
    510509                }
    511510        }

  • ../trunk-jpl/src/c/cores/love_core.cpp

     
    2727
    2828        /*parameters: */
    2929        bool save_results;
    30        
     30
    3131        if(VerboseSolution()) _printf0_("   computing LOVE numbers\n");
    3232
    3333        /*Recover some parameters: */
    3434        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    35        
     35
    3636        /*recover love number parameters: */
    3737        femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
    3838        femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy);
    … …  
    4444        femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
    4545        femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
    4646        femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
    47        
     47
    4848        /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
    4949        Matlitho* matlitho=NULL;
    5050        for (int i=femmodel->materials->Size()-1;i>=0;i--){

  • ../trunk-jpl/src/c/cores/damage_core.cpp

     
    1818        int    numoutputs          = 0;
    1919        char   **requested_outputs = NULL;
    2020
    21                        
    2221        //first recover parameters common to all solutions
    2322        femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
    2423        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);

  • ../trunk-jpl/src/c/cores/smb_core.cpp

     
    2828        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
    2929        femmodel->parameters->FindParam(&numoutputs,SmbNumRequestedOutputsEnum);
    3030        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SmbRequestedOutputsEnum);
    31                        
     31
    3232        if(VerboseSolution()) _printf0_("   computing smb \n");
    33  
     33
    3434        analysis = new SmbAnalysis();
    3535        analysis->Core(femmodel);
    3636        delete analysis;

  • ../trunk-jpl/src/c/cores/ad_core.cpp

     
    4545
    4646                        /*First, stop tracing: */
    4747                        trace_off();
    48                        
     48
    4949                        /*Print tape statistics so that user can kill this run if something is off already:*/
    5050                        if(VerboseAutodiff()){ /*{{{*/
    5151                                tapestats(my_rank,tape_stats); //reading of tape statistics
    … …  
    9292                        /*retrieve parameters: */
    9393                        femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
    9494                        femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
    95        
     95
    9696                        /*if no dependents, no point in running a driver: */
    9797                        if(!(num_dependents*num_independents)) return;
    9898
    … …  
    100100                        if (my_rank!=0){
    101101                                num_dependents=0; num_independents=0;
    102102                        }
    103                        
     103
    104104                        /*retrieve state variable: */
    105105                        femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
    106106
    … …  
    116116                        /*Branch according to AD driver: */
    117117                        femmodel->parameters->FindParam(&driver,AutodiffDriverEnum);
    118118
    119 
    120119                        if (strcmp(driver,"fos_forward")==0){ /*{{{*/
    121120
    122121                                int     anIndepIndex;
    … …  
    140139                                anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
    141140#endif
    142141
    143 
    144142                                /*call driver: */
    145143                                fos_forward(my_rank,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir );
    146144
    … …  
    184182#endif
    185183                                // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
    186184
    187 
    188185                                /*seed matrix: */
    189186                                seed=xNewZeroInit<double>(num_independents,tangentDirNum);
    190187
    … …  
    243240                                anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
    244241#endif
    245242
    246 
    247243                                /*call driver: */
    248244                                fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
    249245
    … …  
    284280                                anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
    285281                                #endif
    286282
    287 
    288283                                /*seed matrix: */
    289284                                weights=xNewZeroInit<double>(weightNum,num_dependents);
    290285
    … …  
    316311                        } /*}}}*/
    317312                        else _error_("driver: " << driver << " not yet supported!");
    318313
    319 
    320314                        if(VerboseAutodiff())_printf0_("   end AD core\n");
    321315
    322316                        /*Free resources: */

  • ../trunk-jpl/src/c/classes/Materials/Matestar.cpp

     
    108108void      Matestar::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    109109
    110110        if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook();
    111        
     111
    112112        MARSHALLING_ENUM(MatestarEnum);
    113113        MARSHALLING(mid);
    114114        this->helement->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);

  • ../trunk-jpl/src/c/classes/Materials/Matlitho.cpp

     
    3636
    3737        this->radius=xNew<IssmDouble>(this->numlayers+1);
    3838        xMemCpy<IssmDouble>(this->radius, iomodel->Data("md.materials.radius"),this->numlayers+1);
    39        
     39
    4040        this->viscosity=xNew<IssmDouble>(this->numlayers);
    4141        xMemCpy<IssmDouble>(this->viscosity, iomodel->Data("md.materials.viscosity"),this->numlayers);
    42        
     42
    4343        this->lame_lambda=xNew<IssmDouble>(this->numlayers);
    4444        xMemCpy<IssmDouble>(this->lame_lambda, iomodel->Data("md.materials.lame_lambda"),this->numlayers);
    45        
     45
    4646        this->lame_mu=xNew<IssmDouble>(this->numlayers);
    4747        xMemCpy<IssmDouble>(this->lame_mu, iomodel->Data("md.materials.lame_mu"),this->numlayers);
    4848
    … …  
    6060
    6161        this->issolid=xNew<IssmDouble>(this->numlayers);
    6262        xMemCpy<IssmDouble>(this->issolid, iomodel->Data("md.materials.issolid"),this->numlayers);
    63        
     63
    6464        /*isburgersd= xNew<IssmDouble>(this->numlayers);
    6565        this->isburgers=xNew<bool>(this->numlayers);
    6666        xMemCpy<IssmDouble>(isburgersd, iomodel->Data("md.materials.isburgers"),this->numlayers);
    6767        for (int i=0;i<this->numlayers;i++)this->isburgers[i]=reCast<bool,IssmDouble>(isburgersd[i]);
    68        
     68
    6969        issolidd= xNew<IssmDouble>(this->numlayers);
    7070        this->issolid=xNew<bool>(this->numlayers);
    7171        xMemCpy<IssmDouble>(issolidd, iomodel->Data("md.materials.issolid"),this->numlayers);
    7272        for (int i=0;i<this->numlayers;i++)this->issolid[i]=reCast<bool,IssmDouble>(issolidd[i]);*/
    73        
     73
    7474        /*free ressources: */
    7575        xDelete<IssmDouble>(isburgersd);
    7676        xDelete<IssmDouble>(issolidd);
    … …  
    7777}
    7878/*}}}*/
    7979Matlitho::~Matlitho(){/*{{{*/
    80        
     80
    8181        xDelete<IssmDouble>(radius);
    8282        xDelete<IssmDouble>(viscosity);
    8383        xDelete<IssmDouble>(lame_lambda);

  • ../trunk-jpl/src/c/classes/Materials/Matice.cpp

     
    134134        _printf_("              note: helement not printed to avoid recursion.\n");
    135135        //if(helement) helement->DeepEcho();
    136136        //else _printf_("   helement = NULL\n");
    137        
     137
    138138        _printf_("   element:\n");
    139139        _printf_("     note: element not printed to avoid recursion.\n");
    140140        //if(element) element->DeepEcho();
    … …  
    142142}               
    143143/*}}}*/
    144144void      Matice::Echo(void){/*{{{*/
    145        
     145
    146146        _printf_("Matice:\n");
    147147        _printf_("   mid: " << mid << "\n");
    148148        _printf_("   isdamaged: " << isdamaged << "\n");
    149149        _printf_("   isenhanced: " << isenhanced << "\n");
    150        
     150
    151151        /*helement and element Echo were commented to avoid recursion.*/
    152152        /*Example: element->Echo calls matice->Echo which calls element->Echo etc*/
    153153        _printf_("   helement:\n");
    … …  
    154154        _printf_("     note: helement not printed to avoid recursion.\n");
    155155        //if(helement) helement->Echo();
    156156        //else _printf_("   helement = NULL\n");
    157        
     157
    158158        _printf_("   element:\n");
    159159        _printf_("     note: element not printed to avoid recursion.\n");
    160160        //if(element) element->Echo();
    … …  
    166166void      Matice::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    167167
    168168        if(marshall_direction==MARSHALLING_BACKWARD)helement=new Hook();
    169        
     169
    170170        MARSHALLING_ENUM(MaticeEnum);
    171171        MARSHALLING(mid);
    172172        MARSHALLING(isdamaged);
    … …  
    540540        /*input strain rate: */
    541541        IssmDouble exx,eyy,exy,exz,eyz;
    542542
    543 
    544543        if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) &&
    545544                                (epsilon[3]==0) && (epsilon[4]==0)){
    546545                mu_prime=0.5*pow(10.,14);

  • ../trunk-jpl/src/c/classes/Numberedcostfunction.cpp

     
    88#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
    99#endif
    1010
    11 
    1211/*Headers:*/
    1312//#include "./Definition.h"
    1413//#include "../datastructures/datastructures.h"
    … …  
    2928#include "../modules/DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h"
    3029
    3130/*}}}*/
    32                
     31
    3332/*Numberedcostfunction constructors, destructors :*/
    3433Numberedcostfunction::Numberedcostfunction(){/*{{{*/
    3534
    … …  
    109108}
    110109/*}}}*/
    111110IssmDouble Numberedcostfunction::Response(FemModel* femmodel){/*{{{*/
    112        
     111
    113112         _assert_(number_cost_functions>0 && number_cost_functions<1e3);
    114113
    115114         /*output:*/
    … …  
    156155        return value_sum;
    157156 }
    158157 /*}}}*/
    159 

  • ../trunk-jpl/src/c/classes/Params/BoolParam.cpp

     
    6262
    6363}
    6464/*}}}*/
    65                

  • ../trunk-jpl/src/c/classes/Params/StringParam.cpp

     
    5454        int size = 0;
    5555
    5656        if(marshall_direction==MARSHALLING_FORWARD || marshall_direction == MARSHALLING_SIZE)size=strlen(value)+1;
    57        
     57
    5858        MARSHALLING_ENUM(StringParamEnum);
    5959        MARSHALLING(enum_type);
    6060        MARSHALLING(size);

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

     
    115115
    116116/*DoubleVecParam specific routines:*/
    117117void  DoubleVecParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM){/*{{{*/
    118        
     118
    119119        /*Assign output pointers:*/
    120120        if(pM) *pM=M;
    121121        *pIssmDoublearray=values;

  • ../trunk-jpl/src/c/classes/Params/DoubleMatParam.cpp

     
    115115
    116116/*DoubleMatParam specific routines:*/
    117117void  DoubleMatParam::GetParameterValueByPointer(IssmDouble** pIssmDoublearray,int* pM,int* pN){/*{{{*/
    118        
     118
    119119        /*Assign output pointers:*/
    120120        if(pM) *pM=M;
    121121        if(pN) *pN=N;

  • ../trunk-jpl/src/c/classes/Params/DataSetParam.cpp

     
    5454void DataSetParam::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    5555
    5656        if(marshall_direction==MARSHALLING_BACKWARD)value=new DataSet();
    57        
     57
    5858        MARSHALLING_ENUM(DataSetParamEnum);
    5959        MARSHALLING(enum_type);
    6060        value->Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);

  • ../trunk-jpl/src/c/classes/Params/IntParam.cpp

     
    6363
    6464}
    6565/*}}}*/
    66 

  • ../trunk-jpl/src/c/classes/Params/StringArrayParam.cpp

     
    8282        }       
    8383
    8484        MARSHALLING_ENUM(StringArrayParamEnum);
    85        
     85
    8686        MARSHALLING(enum_type);
    8787        MARSHALLING(numstrings);
    8888

  • ../trunk-jpl/src/c/classes/Params/DoubleMatArrayParam.cpp

     
    232232
    233233}
    234234/*}}}*/
    235                

  • ../trunk-jpl/src/c/classes/AmrBamg.cpp

     
    3333        this->deviatoricerror_threshold                 = -1;
    3434        this->deviatoricerror_groupthreshold    = -1;
    3535        this->deviatoricerror_maximum                           = -1;
    36        
     36
    3737        /*Geometry and mesh as NULL*/
    3838        this->geometry                                                                  = NULL;
    3939        this->fathermesh                                                                = NULL;
    … …  
    136136        this->options->hmaxVertices = NULL;
    137137        this->options->hmaxVerticesSize[0] = 0;
    138138        this->options->hmaxVerticesSize[1] = 0;
    139        
     139
    140140        /*verify if the metric will be reseted or not*/
    141141        if(this->keepmetric==0){
    142142                if(this->options->metric) xDelete<IssmDouble>(this->options->metric);
    … …  
    154154        int *datalist           = xNew<int>(2);
    155155        IssmDouble *xylist= xNew<IssmDouble>(nbv*2);
    156156        int* elementslist = xNew<int>(nbt*3);
    157        
     157
    158158        datalist[0] = nbv;
    159159        datalist[1] = nbt;
    160        
     160
    161161        for(int i=0;i<nbv;i++){
    162162                xylist[2*i]             = meshout->Vertices[i*3+0];
    163163                xylist[2*i+1]   = meshout->Vertices[i*3+1];
    164164        }
    165        
     165
    166166        for(int i=0;i<nbt;i++){
    167167                elementslist[3*i+0] = reCast<int>(meshout->Triangles[4*i+0]);
    168168                elementslist[3*i+1] = reCast<int>(meshout->Triangles[4*i+1]);
    … …  
    173173        *pdatalist              = datalist;
    174174        *pxylist                        = xylist;
    175175        *pelementslist = elementslist;
    176        
     176
    177177        /*Cleanup and return*/
    178178        delete geomout;
    179179}/*}}}*/
    … …  
    180180void AmrBamg::SetBamgOpts(IssmDouble hmin_in,IssmDouble hmax_in,IssmDouble err_in,IssmDouble gradation_in){/*{{{*/
    181181
    182182        if(!this->options) _error_("AmrBamg->options is NULL!");
    183        
     183
    184184        this->options->hmin     = hmin_in;
    185185        this->options->hmax     = hmax_in;
    186186        this->options->err[0]   = err_in;

  • ../trunk-jpl/src/c/classes/kriging/Observations.cpp

     
    740740        /*Assign output pointer*/
    741741        xDelete<IssmPDouble>(counter);
    742742}/*}}}*/
    743 

  • ../trunk-jpl/src/c/classes/Cfsurfacelogvel.cpp

     
    2222#include "../classes/Inputs/Input.h"
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Cfsurfacelogvel constructors, destructors :*/
    2727Cfsurfacelogvel::Cfsurfacelogvel(){/*{{{*/
    2828
    … …  
    3939Cfsurfacelogvel::Cfsurfacelogvel(char* in_name, int in_definitionenum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/
    4040
    4141        this->definitionenum=in_definitionenum;
    42        
     42
    4343        this->name              = xNew<char>(strlen(in_name)+1);
    4444        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    4545
    4646        this->datatime=in_datatime;
    4747        this->timepassedflag=in_timepassedflag;
    48        
     48
    4949        this->misfit=0;
    5050        this->lock=0;
    5151}
    … …  
    9999}
    100100/*}}}*/
    101101IssmDouble Cfsurfacelogvel::Response(FemModel* femmodel){/*{{{*/
    102                  
     102
    103103         /*diverse: */
    104104         IssmDouble time;
    105          
     105
    106106         /*recover time parameters: */
    107107         femmodel->parameters->FindParam(&time,TimeEnum);
    108108         if(time < last_time) timepassedflag = false;
    … …  
    111111                 int i;
    112112                 IssmDouble J=0.;
    113113                 IssmDouble J_sum=0.;
    114        
     114
    115115         if(datatime<=time && !timepassedflag){
    116116                 for(i=0;i<femmodel->elements->Size();i++){
    117117                         Element* element=(Element*)femmodel->elements->GetObjectByOffset(i);
    … …  
    121121                 ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    122122                 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    123123                 J=J_sum;
    124                
     124
    125125                 timepassedflag = true;
    126126                 return J;
    127127                }
    … …  
    138138        IssmDouble misfit,Jdet;
    139139        IssmDouble vx,vy,vxobs,vyobs,weight;
    140140        IssmDouble* xyz_list = NULL;
    141        
     141
    142142        /*Get basal element*/
    143143        if(!element->IsOnSurface()) return 0.;
    144144
    … …  
    159159
    160160        /* Get node coordinates*/
    161161        topelement->GetVerticesCoordinates(&xyz_list);
    162        
     162
    163163        /*Get model values*/
    164164        Input* vx_input     =topelement->GetInput(VxEnum);                                 _assert_(vx_input);
    165165        Input* vy_input   =NULL;
    … …  
    174174        if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
    175175        datasetinput = (DatasetInput*)tempinput;
    176176
    177 
    178177        /* Start  looping on the number of gaussian points: */
    179178        Gauss* gauss=topelement->NewGauss(2);
    180179        for(int ig=gauss->begin();ig<gauss->end();ig++){
    … …  
    221220        delete gauss;
    222221        return Jelem;
    223222}/*}}}*/
    224 

  • ../trunk-jpl/src/c/classes/Cfdragcoeffabsgrad.cpp

     
    2222#include "../classes/Inputs/Input.h"
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Cfdragcoeffabsgrad constructors, destructors :*/
    2727Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(){/*{{{*/
    2828
    … …  
    3939Cfdragcoeffabsgrad::Cfdragcoeffabsgrad(char* in_name, int in_definitionenum, int in_weights_enum, bool in_timepassedflag){/*{{{*/
    4040
    4141        this->definitionenum=in_definitionenum;
    42        
     42
    4343        this->name              = xNew<char>(strlen(in_name)+1);
    4444        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    4545
    4646        this->weights_enum=in_weights_enum;
    4747        this->timepassedflag=in_timepassedflag;
    48        
     48
    4949        this->misfit=0;
    5050        this->lock=0;
    5151}
    … …  
    101101IssmDouble Cfdragcoeffabsgrad::Response(FemModel* femmodel){/*{{{*/
    102102         /*diverse: */
    103103         IssmDouble time;
    104          
     104
    105105         /*recover time parameters: */
    106106                 int i;
    107107                 IssmDouble J=0.;
    … …  
    115115                 ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    116116                 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    117117                 J=J_sum;
    118                
     118
    119119                 timepassedflag = true;
    120120                 return J;
    121121 }
    … …  
    182182        delete gauss;
    183183        return Jelem;
    184184}/*}}}*/
    185 

  • ../trunk-jpl/src/c/classes/Loads/Neumannflux.cpp

     
    3030/*}}}*/
    3131Neumannflux::Neumannflux(int neumannflux_id,int i,IoModel* iomodel,int* segments,int in_analysis_type){/*{{{*/
    3232
    33 
    3433        /*Some sanity checks*/
    3534        _assert_(segments);
    3635

  • ../trunk-jpl/src/c/classes/Loads/Moulin.cpp

     
    177177        this->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
    178178
    179179        switch(analysis_type){
    180                
     180
    181181        case HydrologyShaktiAnalysisEnum:
    182182                pe = this->CreatePVectorHydrologyShakti();
    183183                break;
    … …  
    394394        if(!this->node->IsActive()) return NULL;
    395395        IssmDouble moulin_load,dt;
    396396        ElementVector* pe=new ElementVector(&node,1,this->parameters);
    397        
     397
    398398        this->element->GetInputValue(&moulin_load,node,HydrologydcBasalMoulinInputEnum);
    399399        parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    400        
     400
    401401        pe->values[0]=moulin_load*dt;
    402402        /*Clean up and return*/
    403403        return pe;

  • ../trunk-jpl/src/c/classes/Cfsurfacesquare.cpp

     
    2222#include "../classes/Inputs/Input.h"
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Cfsurfacesquare constructors, destructors :*/
    2727Cfsurfacesquare::Cfsurfacesquare(){/*{{{*/
    2828
    … …  
    4242Cfsurfacesquare::Cfsurfacesquare(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, int in_weights_enum, IssmDouble in_datatime, bool in_timepassedflag){/*{{{*/
    4343
    4444        this->definitionenum=in_definitionenum;
    45        
     45
    4646        this->name              = xNew<char>(strlen(in_name)+1);
    4747        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    4848
    … …  
    5151        this->weights_enum=in_weights_enum;
    5252        this->datatime=in_datatime;
    5353        this->timepassedflag=in_timepassedflag;
    54        
     54
    5555        this->misfit=0;
    5656        this->lock=0;
    5757}
    … …  
    110110IssmDouble Cfsurfacesquare::Response(FemModel* femmodel){/*{{{*/
    111111         /*diverse: */
    112112         IssmDouble time;
    113          
     113
    114114         /*recover time parameters: */
    115115         femmodel->parameters->FindParam(&time,TimeEnum);
    116116         if(time < last_time) timepassedflag = false;
    … …  
    129129                 ISSM_MPI_Allreduce ( (void*)&J,(void*)&J_sum,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    130130                 ISSM_MPI_Bcast(&J_sum,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    131131                 J=J_sum;
    132                
     132
    133133                 timepassedflag = true;
    134134                 return J;
    135135        }
    … …  
    171171
    172172        /*Get input if it already exists*/
    173173        Input*  tempinput = topelement->GetInput(definitionenum);
    174        
     174
    175175        /*Cast it to a Datasetinput*/
    176176        if(tempinput->ObjectEnum()!=DatasetInputEnum) _error_("don't know what to do");
    177177        datasetinput = (DatasetInput*)tempinput;
    … …  
    213213        delete gauss;
    214214        return Jelem;
    215215}/*}}}*/
    216 

  • ../trunk-jpl/src/c/classes/Profiler.cpp

     
    9999        _assert_(tag<MAXIMUMSIZE);
    100100        if(this->running[tag]) _error_("Tag "<<tag<<" has not been stopped");
    101101
    102 
    103102        #ifdef _HAVE_MPI_
    104103        return this->time[tag];
    105104        #else
    … …  
    143142        _assert_(tag<MAXIMUMSIZE);
    144143        if(this->running[tag]) _error_("Tag "<<tag<<" is already running");
    145144
    146 
    147145        /*If mpisync requested, make sure all the cpus are at the same point in the execution: */
    148146        if(!dontmpisync){
    149147                ISSM_MPI_Barrier(IssmComm::GetComm());
    … …  
    183181        _assert_(tag<MAXIMUMSIZE);
    184182        if(!this->running[tag]) _error_("Tag "<<tag<<" is not running");
    185183
    186 
    187184        /*If mpisync requested, make sure all the cpus are at the same point in the execution: */
    188185        if(!dontmpisync){
    189186                ISSM_MPI_Barrier(IssmComm::GetComm());

  • ../trunk-jpl/src/c/classes/AdaptiveMeshRefinement.cpp

     
    100100
    101101/*Mesh refinement methods*/
    102102void AdaptiveMeshRefinement::ExecuteRefinement(double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror,int** pdatalist,double** pxylist,int** pelementslist){/*{{{*/
    103        
     103
    104104        /*IMPORTANT! pelementslist are in Matlab indexing*/
    105105        /*NEOPZ works only in C indexing*/
    106106        if(!this->fathermesh || !this->previousmesh) _error_("Impossible to execute refinement: fathermesh or previousmesh is NULL!\n");
    … …  
    117117
    118118        /*Intermediaries*/
    119119        bool verbose=VerboseSolution();
    120        
     120
    121121        /*Execute refinement*/
    122122        this->RefineMeshOneLevel(verbose,gl_distance,if_distance,deviatoricerror,thicknesserror);
    123    
     123
    124124        /*Get new geometric mesh in ISSM data structure*/
    125125        this->GetMesh(this->previousmesh,pdatalist,pxylist,pelementslist);
    126126
    … …  
    129129}
    130130/*}}}*/
    131131void AdaptiveMeshRefinement::RefineMeshOneLevel(bool &verbose,double* gl_distance,double* if_distance,double* deviatoricerror,double* thicknesserror){/*{{{*/
    132        
     132
    133133        /*Intermediaries*/
    134134        int nelem                                                       =-1;
    135135        int side2D                                                      = 6;
    … …  
    203203                if(criteria) index.push_back(gmesh->Element(this->specialelementsindex[i])->FatherIndex());
    204204        }
    205205        /*}}}*/
    206        
     206
    207207        /*Now, detele the special elements{{{*/
    208208        if(this->refinement_type==1) this->DeleteSpecialElements(verbose,gmesh);
    209209        else this->specialelementsindex.clear();
    … …  
    215215                gmesh=this->fathermesh;
    216216        }
    217217        /*}}}*/
    218        
     218
    219219        /*Unrefinement process: loop over level of refinements{{{*/
    220220        if(verbose) _printf_("\tunrefinement process...\n");
    221221        if(verbose) _printf_("\ttotal: ");
    222222        count=0;
    223          
     223
    224224        nelem=gmesh->NElements();//must keep here
    225225        for(int i=0;i<nelem;i++){
    226226                if(!gmesh->Element(i)) continue;
    … …  
    265265        /*Adjust the connectivities before continue*/
    266266        gmesh->BuildConnectivity();
    267267        /*}}}*/
    268        
     268
    269269        /*Refinement process: loop over level of refinements{{{*/
    270270        if(verbose) _printf_("\trefinement process (level max = "<<this->level_max<<")\n");
    271271        if(verbose) _printf_("\ttotal: ");
    … …  
    319319         * 2 : uniform refinment
    320320         * */
    321321        if(!geoel) _error_("geoel is NULL!\n");
    322        
     322
    323323        /*Output*/
    324324        int type=0;
    325325        int count=0;
    326        
     326
    327327        /*Intermediaries*/
    328328        TPZVec<TPZGeoEl*> sons;
    329        
     329
    330330        /*Loop over neighboors (sides 3, 4 and 5)*/
    331331        for(int j=3;j<6;j++){
    332332                sons.clear();
    … …  
    334334                if(sons.size()) count++; //if neighbour was refined
    335335                if(sons.size()>4) count++; //if neighbour's level is > element level+1
    336336        }
    337        
     337
    338338        /*Verify and return*/
    339339        if(count>1) type=2;
    340340        else type=count;
    341        
     341
    342342        return type;
    343343}
    344344/*}}}*/
    … …  
    390390}
    391391/*}}}*/
    392392void AdaptiveMeshRefinement::RefineMeshToAvoidHangingNodes(bool &verbose,TPZGeoMesh* gmesh){/*{{{*/
    393    
     393
    394394        /*Now, insert special elements to avoid hanging nodes*/
    395395        if(verbose) _printf_("\trefining to avoid hanging nodes (total: ");
    396        
     396
    397397        /*Intermediaries*/
    398398        int nelem=-1;
    399399        int count= 1;
    400        
     400
    401401        while(count>0){
    402402                nelem=gmesh->NElements();//must keep here
    403403                count=0;
    … …  
    467467        long* vertex_index2sid  = xNew<long>(ntotalvertices);
    468468        this->index2sid.clear(); this->index2sid.resize(gmesh->NElements());
    469469        this->sid2index.clear();
    470        
     470
    471471        /*Fill in the vertex_index2sid vector with non usual index value*/
    472472        for(int i=0;i<gmesh->NNodes();i++) vertex_index2sid[i]=-1;
    473        
     473
    474474        /*Fill in the this->index2sid vector with non usual index value*/
    475475        for(int i=0;i<gmesh->NElements();i++) this->index2sid[i]=-1;
    476        
     476
    477477        /*Get elements without sons and fill in the vertex_index2sid with used vertices (indexes) */
    478478        sid=0;
    479479        for(int i=0;i<gmesh->NElements();i++){//over gmesh elements index
    … …  
    510510                newmeshXY[2*sid]                = coords[0]; // X
    511511                newmeshXY[2*sid+1]      = coords[1]; // Y
    512512        }
    513                
     513
    514514        for(int i=0;i<this->sid2index.size();i++){//over the sid (fill the ISSM elements)
    515515                for(int j=0;j<this->GetNumberOfNodes();j++) {
    516516                        geoel   = gmesh->ElementVec()[this->sid2index[i]];
    … …  
    532532                xc=newmeshXY[2*c]; yc=newmeshXY[2*c+1];
    533533
    534534                detJ=(xb-xa)*(yc-ya)-(xc-xa)*(yb-ya);
    535        
     535
    536536                /*verify and swap, if necessary*/
    537537                if(detJ<0) {
    538538                        newelements[i*this->GetNumberOfNodes()+0]=b+1;//a->b
    … …  
    544544        *pdata          = newdata;
    545545        *pxy               = newmeshXY;
    546546        *pelements      = newelements;
    547    
     547
    548548        /*Cleanup*/
    549549        xDelete<long>(vertex_index2sid);
    550550}
    … …  
    553553
    554554        /* IMPORTANT! elements come in Matlab indexing
    555555                NEOPZ works only in C indexing*/
    556        
     556
    557557        if(nvertices<=0) _error_("Impossible to create initial mesh: nvertices is <= 0!\n");
    558558   if(nelements<=0) _error_("Impossible to create initial mesh: nelements is <= 0!\n");
    559559        if(this->refinement_type!=0 && this->refinement_type!=1) _error_("Impossible to create initial mesh: refinement type is not defined!\n");
    … …  
    560560
    561561    /*Verify and creating initial mesh*/
    562562   if(this->fathermesh || this->previousmesh) _error_("Initial mesh already exists!");
    563    
     563
    564564   this->fathermesh = new TPZGeoMesh();
    565565        this->fathermesh->NodeVec().Resize(nvertices);
    566566
    … …  
    575575      this->fathermesh->NodeVec()[i].SetCoord(coord);
    576576                this->fathermesh->NodeVec()[i].SetNodeId(i);
    577577        }
    578        
     578
    579579        /*Generate the elements*/
    580580        long index;
    581581   const int mat = this->GetElemMaterialID();
    … …  
    603603}
    604604/*}}}*/
    605605TPZGeoMesh* AdaptiveMeshRefinement::CreateRefPatternMesh(TPZGeoMesh* gmesh){/*{{{*/
    606        
     606
    607607        TPZGeoMesh *newgmesh = new TPZGeoMesh();
    608608   newgmesh->CleanUp();
    609    
     609
    610610   int nnodes  = gmesh->NNodes();
    611611        int nelem   = gmesh->NElements();
    612612   int mat     = this->GetElemMaterialID();;
    … …  
    616616        //nodes
    617617        newgmesh->NodeVec().Resize(nnodes);
    618618   for(int i=0;i<nnodes;i++) newgmesh->NodeVec()[i] = gmesh->NodeVec()[i];
    619    
     619
    620620   //elements
    621621   for(int i=0;i<nelem;i++){
    622622        TPZGeoEl * geoel = gmesh->Element(i);
    623                
     623
    624624                if(!geoel){
    625625                        index=newgmesh->ElementVec().AllocateNewElement();
    626626                        newgmesh->ElementVec()[index] = NULL;
    627627                        continue;
    628628                }
    629      
     629
    630630                TPZManVector<long> elem(3,0);
    631631      for(int j=0;j<3;j++) elem[j] = geoel->NodeIndex(j);
    632      
     632
    633633      newgmesh->CreateGeoElement(ETriangle,elem,mat,index,reftype);
    634634                newgmesh->ElementVec()[index]->SetId(geoel->Id());
    635        
     635
    636636      TPZGeoElRefPattern<TPZGeoTriangle>* newgeoel = dynamic_cast<TPZGeoElRefPattern<TPZGeoTriangle>*>(newgmesh->ElementVec()[index]);
    637        
     637
    638638      //old neighbourhood
    639639      const int nsides = TPZGeoTriangle::NSides;
    640640      TPZVec< std::vector<TPZGeoElSide> > neighbourhood(nsides);
    … …  
    653653            neighS = neighS.Neighbour();
    654654         }
    655655      }
    656        
     656
    657657      //inserting in new element
    658658      for(int s = 0; s < nsides; s++){
    659659        TPZGeoEl * tempEl = newgeoel;
    … …  
    667667         }
    668668         tempEl->SetNeighbour(byside, tempSide);
    669669      }
    670        
     670
    671671      long fatherindex = geoel->FatherIndex();
    672672      if(fatherindex>-1) newgeoel->SetFather(fatherindex);
    673        
     673
    674674      if(!geoel->HasSubElement()) continue;
    675        
     675
    676676      int nsons = geoel->NSubElements();
    677677
    678678      TPZAutoPointer<TPZRefPattern> ref = gRefDBase.GetUniformRefPattern(ETriangle);
    679679      newgeoel->SetRefPattern(ref);
    680        
     680
    681681      for(int j=0;j<nsons;j++){
    682682        TPZGeoEl* son = geoel->SubElement(j);
    683683         if(!son){
    … …  
    689689
    690690        /*Now, build connectivities*/
    691691        newgmesh->BuildConnectivity();
    692    
     692
    693693        return newgmesh;
    694694}
    695695/*}}}*/
    … …  
    704704}
    705705/*}}}*/
    706706void AdaptiveMeshRefinement::PrintGMeshVTK(TPZGeoMesh* gmesh,std::ofstream &file,bool matColor){/*{{{*/
    707    
     707
    708708        file.clear();
    709709        long nelements = gmesh->NElements();
    710710        TPZGeoEl *gel;
    … …  
    730730          }
    731731          MElementType elt = gel->Type();
    732732          int elNnodes = MElementType_NNodes(elt);
    733          
     733
    734734          size += (1+elNnodes);
    735735          connectivity << elNnodes;
    736          
     736
    737737          for(int t = 0; t < elNnodes; t++)
    738738          {
    739739                        for(int c = 0; c < 3; c++)
    … …  
    742742                                 node << coord << " ";
    743743                        }
    744744                        node << std::endl;
    745                        
     745
    746746                        actualNode++;
    747747                        connectivity << " " << actualNode;
    748748          }
    749749          connectivity << std::endl;
    750          
     750
    751751          int elType = this->GetVTK_ElType(gel);
    752752          type << elType << std::endl;
    753          
     753
    754754          if(matColor == true)
    755755          {
    756756                        material << gel->MaterialId() << std::endl;
    … …  
    759759          {
    760760                        material << gel->Index() << std::endl;
    761761          }
    762          
     762
    763763          nVALIDelements++;
    764764        }
    765765        node << std::endl;
    … …  
    789789}
    790790/*}}}*/
    791791int AdaptiveMeshRefinement::GetVTK_ElType(TPZGeoEl * gel){/*{{{*/
    792    
     792
    793793        MElementType pzElType = gel->Type();
    794    
     794
    795795    int elType = -1;
    796796    switch (pzElType)
    797797    {
    … …  
    848848        std::cout << "MIGHT BE CURVED ELEMENT (quadratic or quarter point)" << std::endl;
    849849        DebugStop();
    850850    }
    851    
     851
    852852    return elType;
    853853}
    854854/*}}}*/
    855 

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

     
    105105                for(i=0;i<num_nodes;i++) if(this->nodes[i]) tria->nodes[i]=this->nodes[i]; else tria->nodes[i] = NULL;
    106106        }
    107107        else tria->nodes = NULL;
    108        
     108
    109109        tria->vertices = (Vertex**)this->hvertices->deliverp();
    110110        tria->material = (Material*)this->hmaterial->delivers();
    111111        tria->matpar   = (Matpar*)this->hmatpar->delivers();
    … …  
    114114}
    115115/*}}}*/
    116116void Tria::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    117        
     117
    118118        MARSHALLING_ENUM(TriaEnum);
    119        
     119
    120120        /*Call parent classes: */
    121121        ElementHook::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
    122122        Element::MarshallElement(pmarshalled_data,pmarshalled_data_size,marshall_direction,this->numanalyses);
    … …  
    134134
    135135        /*Call inputs method*/
    136136        _assert_(this->inputs);
    137        
     137
    138138        int domaintype;
    139139        parameters->FindParam(&domaintype,DomainTypeEnum);
    140140        switch(domaintype){
    … …  
    184184        if(N!=num_inputs) _error_("sizes are not consistent");
    185185
    186186        int        tria_vertex_ids[3];
    187        
     187
    188188        for(int k=0;k<3;k++){
    189189                tria_vertex_ids[k]=reCast<int>(iomodel->elements[3*this->Sid()+k]); //ids for vertices are in the elements array from Matlab
    190190        }
    … …  
    327327}
    328328/*}}}*/
    329329void       Tria::CalvingCrevasseDepth(){/*{{{*/
    330        
     330
    331331        IssmDouble  xyz_list[NUMVERTICES][3];
    332332        IssmDouble  calvingrate[NUMVERTICES];
    333333        IssmDouble  vx,vy,vel;
    … …  
    339339
    340340        /* Get node coordinates and dof list: */
    341341        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    342                
     342
    343343        /*Get the critical fraction of thickness surface and basal crevasses penetrate for calving onset*/
    344344        this->parameters->FindParam(&critical_fraction,CalvingCrevasseDepthEnum);
    345                
     345
    346346        IssmDouble rho_ice        = this->GetMaterialParameter(MaterialsRhoIceEnum);
    347347        IssmDouble rho_seawater   = this->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    348348        IssmDouble rho_freshwater = this->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    … …  
    360360        Input*   s_xx_input              = inputs->GetInput(DeviatoricStressxxEnum);     _assert_(s_xx_input);
    361361        Input*   s_xy_input              = inputs->GetInput(DeviatoricStressxyEnum);     _assert_(s_xy_input);
    362362        Input*   s_yy_input              = inputs->GetInput(DeviatoricStressyyEnum);     _assert_(s_yy_input);
    363        
     363
    364364        /*Loop over all elements of this partition*/
    365365        GaussTria* gauss=new GaussTria();
    366366        for (int iv=0;iv<NUMVERTICES;iv++){
    367367                gauss->GaussVertex(iv);
    368        
     368
    369369                H_input->GetInputValue(&thickness,gauss);
    370370                bed_input->GetInputValue(&bed,gauss);
    371371                surface_input->GetInputValue(&float_depth,gauss);
    … …  
    377377                s_xx_input->GetInputValue(&s_xx,gauss);
    378378                s_xy_input->GetInputValue(&s_xy,gauss);
    379379                s_yy_input->GetInputValue(&s_yy,gauss);
    380                
     380
    381381                vel=sqrt(vx*vx+vy*vy)+1.e-14;
    382382
    383383                s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
    384384                s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
    385385                if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}
    386                
     386
    387387                Ho = thickness - (rho_seawater/rho_ice) * (-bed);
    388388                if(Ho<0.)  Ho=0.;
    389389
    … …  
    398398                //if (surface_crevasse[iv]<water_height){
    399399                //      water_height = surface_crevasse[iv];
    400400                //}
    401                
     401
    402402                /*basal crevasse*/
    403403                //basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (rheology_B * strainparallel * pow(straineffective,((1/rheology_n)-1)) / (rho_ice*constant_g) - Ho);
    404404                basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
    405405                if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
    406406                if (bed>0.) basal_crevasse[iv] = 0.;
    407        
     407
    408408                H_surf = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height - critical_fraction*float_depth;
    409409                H_surfbasal = (surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv])-(critical_fraction*thickness);
    410                
     410
    411411                crevasse_depth[iv] = max(H_surf,H_surfbasal);
    412412        }
    413        
     413
    414414        this->inputs->AddInput(new TriaInput(SurfaceCrevasseEnum,&surface_crevasse[0],P1Enum));
    415415        this->inputs->AddInput(new TriaInput(BasalCrevasseEnum,&basal_crevasse[0],P1Enum));
    416416        this->inputs->AddInput(new TriaInput(CrevasseDepthEnum,&crevasse_depth[0],P1Enum));
    … …  
    460460                }
    461461                else
    462462                        calvingrate[iv]=0.;
    463                
     463
    464464                calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
    465465                calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
    466466        }
    … …  
    561561        IssmDouble  strain_xy[NUMVERTICES];
    562562        IssmDouble  vorticity_xy[NUMVERTICES];
    563563        GaussTria*  gauss=NULL;
    564        
     564
    565565        /* Get node coordinates and dof list: */
    566566        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    567567
    … …  
    568568        /*Retrieve all inputs we will be needing: */
    569569        Input* vx_input=this->GetInput(EsaXmotionEnum); _assert_(vx_input);
    570570        Input* vy_input=this->GetInput(EsaYmotionEnum); _assert_(vy_input);
    571        
     571
    572572        /* Start looping on the number of vertices: */
    573573        gauss=new GaussTria();
    574574        for (int iv=0;iv<NUMVERTICES;iv++){
    … …  
    772772                }
    773773        }
    774774
    775                
    776775}/*}}}*/
    777776void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
    778777
    … …  
    13991398}
    14001399/*}}}*/
    14011400void       Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    1402        
     1401
    14031402        /* Intermediaries */
    14041403        int i, dir,nrfrontnodes;
    14051404        IssmDouble  levelset[NUMVERTICES];
    … …  
    14881487
    14891488}/*}}}*/
    14901489void       Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/
    1491        
     1490
    14921491        /* GetLevelsetIntersection computes:
    14931492         * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
    14941493         * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/
    … …  
    15641563}
    15651564/*}}}*/
    15661565void       Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
    1567        
     1566
    15681567        /*Computeportion of the element that has a positive levelset*/
    15691568
    15701569        bool               negative=true;
    … …  
    16781677        Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    16791678
    16801679        parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
    1681        
     1680
    16821681        /*Cast to Controlinput*/
    16831682        if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    16841683        ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
    1685        
     1684
    16861685        if(strcmp(data,"value")==0){
    16871686                input  = controlinput->values;
    16881687        }
    … …  
    16991698                _error_("Data " << data << " not supported yet");
    17001699        }
    17011700        /*Check what input we are dealing with*/
    1702        
     1701
    17031702        switch(input->ObjectEnum()){
    17041703                case TriaInputEnum:
    17051704                          {
    … …  
    17161715                                vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
    17171716                                break;
    17181717                          }
    1719                        
     1718
    17201719                case TransientInputEnum:
    17211720                                {
    17221721                                        TransientInput* transientinput = xDynamicCast<TransientInput*>(input);
    … …  
    19811980        surface_input->GetInputAverage(&surface);
    19821981        base_input->GetInputAverage(&bed);
    19831982        bed_input->GetInputAverage(&bathymetry);
    1984        
     1983
    19851984        /*Return: */
    19861985        return base*(surface-bed+min(rho_water/rho_ice*bathymetry,0.));
    19871986}
    … …  
    20262025        if(control_analysis && ad_analysis) iomodel->FindConstant(&num_control_type,"md.autodiff.num_independent_objects");
    20272026        if(control_analysis && ad_analysis) iomodel->FindConstant(&num_responses,"md.autodiff.num_dependent_objects");
    20282027
    2029 
    2030 
    20312028        /*Recover vertices ids needed to initialize inputs*/
    20322029        for(i=0;i<3;i++){
    20332030                tria_vertex_ids[i]=reCast<int>(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab
    … …  
    23932390/*}}}*/
    23942391IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
    23952392
    2396 
    23972393        /*intermediary: */
    23982394        IssmDouble* values=NULL;
    23992395        Input*      thickness_input=NULL;
    … …  
    24062402        int         point1;
    24072403        IssmDouble  fraction1,fraction2;
    24082404        bool        mainlynegative=true;
    2409        
     2405
    24102406        /*Output:*/
    24112407        volume=0;
    24122408
    … …  
    24152411
    24162412        /*Retrieve inputs required:*/
    24172413        thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
    2418        
     2414
    24192415        /*Retrieve material parameters: */
    24202416        rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    24212417
    … …  
    24242420        for(int i=0;i<NUMVERTICES;i++){
    24252421                values[i]=levelset[this->vertices[i]->Sid()];
    24262422        }
    2427                
     2423
    24282424        /*Ok, use the level set values to figure out where we put our gaussian points:*/
    24292425        this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
    24302426        Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
    … …  
    28762872        dist_gl_input->GetInputValue(&dist_gl,gauss);
    28772873        dist_cf_input->GetInputValue(&dist_cf,gauss);
    28782874        delete gauss;
    2879        
     2875
    28802876        /*Ensure values are positive for floating ice*/
    28812877        dist_gl = fabs(dist_gl);
    28822878        dist_cf = fabs(dist_cf);
    … …  
    31963192                }
    31973193        }
    31983194
    3199 
    32003195}
    32013196/*}}}*/
    32023197void       Tria::ResetFSBasalBoundaryCondition(void){/*{{{*/
    … …  
    32773272        IssmDouble  values[NUMVERTICES];
    32783273        int         idlist[NUMVERTICES],control_init;
    32793274
    3280 
    32813275        /*Get Domain type*/
    32823276        int domaintype;
    32833277        parameters->FindParam(&domaintype,DomainTypeEnum);
    … …  
    32973291
    32983292        /*Get out if this is not an element input*/
    32993293        if(!IsInputEnum(control_enum)) return;
    3300        
     3294
    33013295        Input* input     = (Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    33023296        if(input->ObjectEnum()!=ControlInputEnum){
    33033297                _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    … …  
    33303324        IssmDouble  values[NUMVERTICES];
    33313325        int         vertexpidlist[NUMVERTICES],control_init;
    33323326
    3333 
    33343327        /*Get Domain type*/
    33353328        int domaintype;
    33363329        parameters->FindParam(&domaintype,DomainTypeEnum);
    … …  
    41664159        IssmDouble* U_elastic_precomputed = NULL;
    41674160        IssmDouble* H_elastic_precomputed = NULL;
    41684161        int         M, hemi;
    4169        
     4162
    41704163        /*computation of Green functions:*/
    41714164        IssmDouble* U_elastic= NULL;
    41724165        IssmDouble* N_elastic= NULL;
    … …  
    41734166        IssmDouble* E_elastic= NULL;
    41744167        IssmDouble* X_elastic= NULL;
    41754168        IssmDouble* Y_elastic= NULL;
    4176        
     4169
    41774170        /*optimization:*/
    41784171        bool store_green_functions=false;
    41794172
    … …  
    41814174        Input*  deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum);
    41824175        if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
    41834176        deltathickness_input->GetInputAverage(&I);
    4184                
     4177
    41854178        /*early return if we are not on the (ice) loading point: */
    41864179        if(I==0) return;
    41874180
    … …  
    41944187
    41954188        /*which hemisphere? for north-south, east-west components*/
    41964189        this->parameters->FindParam(&hemi,EsaHemisphereEnum);
    4197        
     4190
    41984191        /*compute area of element:*/
    41994192        area=GetArea();
    42004193
    … …  
    42534246                U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i];
    42544247                Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i];
    42554248                X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i];
    4256                
     4249
    42574250                /*North-south, East-west components */
    42584251                if (hemi == -1) {
    42594252                        ang2 = PI/2 - atan2(yy[i],xx[i]);
    … …  
    42764269        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
    42774270        pX->SetValues(gsize,indices,X_values,ADD_VAL);
    42784271        pY->SetValues(gsize,indices,Y_values,ADD_VAL);
    4279        
     4272
    42804273        /*free ressources:*/
    42814274        xDelete<int>(indices);
    42824275        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
    … …  
    43064299        IssmDouble* U_elastic_precomputed = NULL;
    43074300        IssmDouble* H_elastic_precomputed = NULL;
    43084301        int         M;
    4309        
     4302
    43104303        /*computation of Green functions:*/
    43114304        IssmDouble* U_elastic= NULL;
    43124305        IssmDouble* N_elastic= NULL;
    43134306        IssmDouble* E_elastic= NULL;
    4314        
     4307
    43154308        /*optimization:*/
    43164309        bool store_green_functions=false;
    43174310
    … …  
    43194312        Input*  deltathickness_input=inputs->GetInput(EsaDeltathicknessEnum);
    43204313        if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
    43214314        deltathickness_input->GetInputAverage(&I);
    4322                
     4315
    43234316        /*early return if we are not on the (ice) loading point: */
    43244317        if(I==0) return;
    43254318
    … …  
    44184411                dx = x_element-x; dy = y_element-y; dz = z_element-z;
    44194412                N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    44204413                E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    4421                
     4414
    44224415                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    44234416                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    44244417                U_elastic[i] += U_elastic_precomputed[index];
    … …  
    44334426        pUp->SetValues(gsize,indices,U_values,ADD_VAL);
    44344427        pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
    44354428        pEast->SetValues(gsize,indices,E_values,ADD_VAL);
    4436        
     4429
    44374430        /*free ressources:*/
    44384431        xDelete<int>(indices);
    44394432        xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
    … …  
    44544447IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/
    44554448
    44564449        if(IsWaterInElement()){
    4457                
     4450
    44584451                IssmDouble area;
    44594452
    44604453                /*Compute area of element:*/
    … …  
    45274520
    45284521        if(IsWaterInElement()){
    45294522                IssmDouble rho_water, S;
    4530                
     4523
    45314524                /*recover material parameters: */
    45324525                rho_water=matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    4533        
     4526
    45344527                /*From Sg_old, recover water sea level rise:*/
    45354528                S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;
    4536                
     4529
    45374530                /* Perturbation terms for moment of inertia (moi_list):
    45384531                 * computed analytically (see Wu & Peltier, eqs 10 & 32)
    45394532                 * also consistent with my GMD formulation!
    … …  
    45454538        }
    45464539        else if(this->inputs->Max(MaskIceLevelsetEnum)<0){
    45474540                IssmDouble rho_ice, I;
    4548                
     4541
    45494542                /*recover material parameters: */
    45504543                rho_ice=matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    4551        
     4544
    45524545                /*Compute ice thickness change: */
    45534546                Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    45544547                if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
    45554548                deltathickness_input->GetInputAverage(&I);
    4556                
     4549
    45574550                dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
    45584551                dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
    45594552                dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
    45604553        }
    4561        
     4554
    45624555        return;
    45634556}/*}}}*/
    45644557void    Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/
    … …  
    45914584        bool computerigid = true;
    45924585        bool computeelastic= true;
    45934586        bool scaleoceanarea= false;
    4594        
     4587
    45954588        /*early return if we are not on an ice cap:*/
    45964589        if(!(this->inputs->Max(MaskIceLevelsetEnum)<=0)){
    45974590                constant=0; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
    … …  
    46054598                *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
    46064599                return;
    46074600        }
    4608        
     4601
    46094602        /*If we are here, we are on ice that is fully grounded or half-way to floating: */
    46104603        if ((this->inputs->Min(MaskGroundediceLevelsetEnum))<0){
    46114604                notfullygrounded=true; //used later on.
    46124605        }
    4613                
     4606
    46144607        /*Inform mask: */
    46154608        constant=1; this->inputs->AddInput(new TriaInput(SealevelEustaticMaskEnum,&constant,P0Enum));
    46164609
    … …  
    46354628        this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
    46364629
    46374630        /* Where is the centroid of this element?:{{{*/
    4638        
     4631
    46394632        /*retrieve coordinates: */
    46404633        ::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);
    4641        
     4634
    46424635        IssmDouble minlong=400;
    46434636        IssmDouble maxlong=-20;
    46444637        for (int i=0;i<NUMVERTICES;i++){
    … …  
    46524645                if (llr_list[1][1]==0)llr_list[1][1]=360;
    46534646                if (llr_list[2][1]==0)llr_list[2][1]=360;
    46544647        }
    4655        
     4648
    46564649        // correction at the north pole
    46574650        if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    46584651        if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    46594652        if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;
    4660                        
     4653
    46614654        //correction at the south pole
    46624655        if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
    46634656        if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
    … …  
    46834676                phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
    46844677                area*=phi;
    46854678        }
    4686        
     4679
    46874680        /*Compute ice thickness change: */
    46884681        Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    46894682        if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
    … …  
    46954688                bool mainlyfloating = true;
    46964689                int         point1;
    46974690                IssmDouble  fraction1,fraction2;
    4698                
     4691
    46994692                /*Recover portion of element that is grounded*/
    47004693                this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
    47014694                Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
    … …  
    47494742                        values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
    47504743                }
    47514744                pSgi->SetValues(gsize,indices,values,ADD_VAL);
    4752                
     4745
    47534746                /*free ressources:*/
    47544747                xDelete<IssmDouble>(values);
    47554748                xDelete<int>(indices);
    47564749        }
    4757        
     4750
    47584751        /*Assign output pointer:*/
    47594752        _assert_(!xIsNan<IssmDouble>(eustatic));
    47604753        _assert_(!xIsInf<IssmDouble>(eustatic));
    … …  
    47794772        /*precomputed elastic green functions:*/
    47804773        IssmDouble* G_elastic_precomputed = NULL;
    47814774        int         M;
    4782        
     4775
    47834776        /*computation of Green functions:*/
    47844777        IssmDouble* G_elastic= NULL;
    47854778        IssmDouble* G_rigid= NULL;
    … …  
    48564849        late=late/180*PI;
    48574850        longe=longe/180*PI;
    48584851        /*}}}*/
    4859        
     4852
    48604853        if(computeelastic){
    4861        
     4854
    48624855                /*recover elastic green function:*/
    48634856                DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
    48644857                parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
    … …  
    48964889                        values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
    48974890                }
    48984891        }
    4899        
     4892
    49004893        pSgo->SetValues(gsize,indices,values,ADD_VAL);
    49014894
    49024895        /*free ressources:*/
    … …  
    49294922        IssmDouble* U_elastic_precomputed = NULL;
    49304923        IssmDouble* H_elastic_precomputed = NULL;
    49314924        int         M;
    4932        
     4925
    49334926        /*computation of Green functions:*/
    49344927        IssmDouble* U_elastic= NULL;
    49354928        IssmDouble* N_elastic= NULL;
    … …  
    49564949        /*recover computational flags: */
    49574950        this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
    49584951        this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
    4959        
     4952
    49604953        /*early return if elastic not requested:*/
    49614954        if(!computeelastic) return;
    49624955
    … …  
    50255018
    50265019        /*From Sg, recover water sea level rise:*/
    50275020        S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;
    5028        
     5021
    50295022        /*Compute ice thickness change: */
    50305023        Input*  deltathickness_input=inputs->GetInput(SealevelriseDeltathicknessEnum);
    50315024        if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
    50325025        deltathickness_input->GetInputAverage(&I);
    5033                
     5026
    50345027        /*initialize: */
    50355028        U_elastic=xNewZeroInit<IssmDouble>(gsize);
    50365029        if(horiz){
    … …  
    50725065                        N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    50735066                        E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
    50745067                }
    5075                
     5068
    50765069                /*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
    50775070                int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
    50785071                U_elastic[i] += U_elastic_precomputed[index];

  • ../trunk-jpl/src/c/classes/Elements/Penta.cpp

     
    286286        IssmDouble  calvingratey[NUMVERTICES];
    287287        IssmDouble  calvingrate[NUMVERTICES];
    288288
    289 
    290289        /* Get node coordinates and dof list: */
    291290        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    292291
    … …  
    10471046}
    10481047/*}}}*/
    10491048void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    1050        
     1049
    10511050        /* Intermediaries */
    10521051        const int dim=3;
    10531052        int i, dir,nrfrontnodes;
    … …  
    11811180        if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");
    11821181        Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
    11831182
    1184 
    11851183        /*Cast to Controlinput*/
    11861184        if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
    11871185        ControlInput* controlinput = xDynamicCast<ControlInput*>(input);
    … …  
    25942592
    25952593        Penta* penta=this;
    25962594        for(;;){
    2597        
     2595
    25982596                IssmDouble  xyz_list[NUMVERTICES][3];
    25992597                /* Get node coordinates and dof list: */
    26002598                ::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
    … …  
    26032601                Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
    26042602                Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
    26052603                Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
    2606        
     2604
    26072605                /*Retrieve all inputs we will need*/
    26082606                Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
    26092607                Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
    … …  
    26142612                Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum);             _assert_(deviayy_input);
    26152613                Input* surface_input=inputs->GetInput(SurfaceEnum);                                                             _assert_(surface_input);
    26162614                Input* thickness_input=inputs->GetInput(ThicknessEnum);                                                 _assert_(thickness_input);
    2617                
     2615
    26182616                /* Start looping on the number of 2D vertices: */
    26192617                for(int ig=0;ig<3;ig++){
    26202618                        GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
    … …  
    26432641                        }
    26442642                        delete gauss;
    26452643                }
    2646                        
     2644
    26472645                /*Stop if we have reached the surface/base*/
    26482646                if(penta->IsOnSurface()) break;
    2649                
     2647
    26502648                /*get upper Penta*/
    26512649                penta=penta->GetUpperPenta();
    26522650                _assert_(penta->Id()!=this->id);

  • ../trunk-jpl/src/c/classes/Elements/ElementHook.cpp

     
    103103        MARSHALLING_DYNAMIC(hnodesi_null,bool,numanalyses);
    104104
    105105        if(marshall_direction==MARSHALLING_BACKWARD){
    106                
     106
    107107                if (!hnodes_null)this->hnodes = new Hook*[numanalyses];
    108108                else this->hnodes=NULL;
    109109                this->hvertices   = new Hook();
    … …  
    139139
    140140        _printf_(" ElementHook DeepEcho:\n");
    141141        _printf_("  numanalyses : "<< this->numanalyses <<"\n");
    142        
     142
    143143        _printf_("  hnodes:\n");
    144144        if(hnodes){
    145145                for(int i=0;i<this->numanalyses;i++) {
    … …  
    148148                }
    149149        }
    150150        else _printf_("  hnodes = NULL\n");
    151        
     151
    152152        _printf_("  hvertices:\n");
    153153        if(hvertices) hvertices->DeepEcho();
    154154        else _printf_("  hvertices = NULL\n");
    155        
     155
    156156        _printf_("  hmaterial:\n");
    157157        if(hmaterial) hmaterial->DeepEcho();
    158158   else _printf_("  hmaterial = NULL\n");
    … …  
    169169}
    170170/*}}}*/
    171171void ElementHook::Echo(){/*{{{*/
    172        
     172
    173173        _printf_(" ElementHook Echo:\n");
    174174        _printf_("  numanalyses : "<< this->numanalyses <<"\n");
    175        
     175
    176176        _printf_("  hnodes:\n");
    177177        if(hnodes){
    178178                for(int i=0;i<this->numanalyses;i++) {
    … …  
    180180                }
    181181        }
    182182        else _printf_("  hnodes = NULL\n");
    183        
     183
    184184        _printf_("  hvertices:\n");
    185185        if(hvertices) hvertices->Echo();
    186186        else _printf_("  hvertices = NULL\n");
    187        
     187
    188188        _printf_("  hmaterial:\n");
    189189        if(hmaterial) hmaterial->Echo();
    190190   else _printf_("  hmaterial = NULL\n");

  • ../trunk-jpl/src/c/classes/Elements/Seg.cpp

     
    103103
    104104        return seg;
    105105
    106 
    107106}
    108107/*}}}*/
    109108void Seg::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    … …  
    140139}
    141140/*}}}*/
    142141void       Seg::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
    143        
     142
    144143        /* Intermediaries */
    145144        int nrfrontnodes,index;
    146145        IssmDouble  levelset[NUMVERTICES];

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

     
    11041104/*}}}*/
    11051105void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/
    11061106
    1107 
    11081107        /*Get number of nodes for this element*/
    11091108        int numnodes = this->GetNumberOfNodes();
    11101109
    … …  
    11211120                if(ug[nodes[i]->GetDof(0,GsetEnum)] > input_max) input_max = ug[nodes[i]->GetDof(0,GsetEnum)];
    11221121        }
    11231122
    1124 
    11251123        /*Second loop to reassign min and max with local extrema*/
    11261124        for(int i=0;i<numnodes;i++){
    11271125                if(min[nodes[i]->GetDof(0,GsetEnum)]>input_min) min[nodes[i]->GetDof(0,GsetEnum)] = input_min;
    … …  
    13881386        default:
    13891387                _error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
    13901388        }
    1391        
     1389
    13921390        /*Clean up*/
    13931391        xDelete<int>(doflist);
    13941392        xDelete<IssmDouble>(values);
    … …  
    14761474        parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
    14771475        parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
    14781476
    1479 
    14801477        /*Get number of vertices*/
    14811478        int numvertices = this->GetNumberOfVertices();
    14821479        if(isautodiff){
    … …  
    18071804           this->inputs->AddInput(datasetinput);
    18081805        }
    18091806
    1810 
    18111807    /*Branch on type of vector: nodal or elementary: */
    18121808    if(vector_type==1){ //nodal vector
    18131809
    … …  
    20662062        this->GetInputListOnVertices(deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);
    20672063        this->GetInputListOnVertices(deepwaterel,BasalforcingsDeepwaterElevationEnum);
    20682064        this->GetInputListOnVertices(upperwaterel,BasalforcingsUpperwaterElevationEnum);
    2069        
     2065
    20702066        for(int i=0;i<numvertices;i++){
    20712067                if(base[i]>upperwaterel[i])      values[i]=0;
    20722068                else if (base[i]<deepwaterel[i]) values[i]=deepwatermelt[i];
    … …  
    28692865        /*only compute SMB at the surface: */
    28702866        if (!IsOnSurface()) return;
    28712867
    2872 
    28732868        /*Retrieve material properties and parameters:{{{ */
    28742869        rho_ice = matpar->GetMaterialParameter(MaterialsRhoIceEnum);
    28752870        rho_water = matpar->GetMaterialParameter(MaterialsRhoFreshwaterEnum);
    … …  
    30713066                /*Snow, firn and ice albedo:*/
    30723067                if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
    30733068
    3074 
    30753069                /*Distribution of absorbed short wave radation with depth:*/
    30763070                if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
    30773071

  • ../trunk-jpl/src/c/classes/gauss/GaussSeg.cpp

     
    3030        /*Get gauss points*/
    3131        this->numgauss = order;
    3232        GaussLegendreLinear(&pcoords1,&pweights,order);
    33        
     33
    3434        this->coords1=xNew<IssmDouble>(numgauss);
    3535        this->weights=xNew<IssmDouble>(numgauss);
    3636

  • ../trunk-jpl/src/c/classes/Regionaloutput.cpp

     
    154154        return val_t;
    155155}
    156156/*}}}*/
    157 

  • ../trunk-jpl/src/c/classes/Dakota/IssmParallelDirectApplicInterface.cpp

     
    1919
    2020                int world_rank;
    2121                ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank);
    22                
     22
    2323                /*Build an femmodel if you are a slave, using the corresponding communicator:*/
    2424                if(world_rank!=0){
    2525                        femmodel_init= new FemModel(argc,argv,evaluation_comm);
    … …  
    3232
    3333                int world_rank;
    3434                ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&world_rank);
    35                
     35
    3636                if(world_rank!=0){
    3737
    3838                        /*Wrap up: */

  • ../trunk-jpl/src/c/classes/Node.cpp

     
    497497}
    498498/*}}}*/
    499499void  Node::SetApproximation(int in_approximation){/*{{{*/
    500        
     500
    501501        this->approximation = in_approximation;
    502502}
    503503/*}}}*/

  • ../trunk-jpl/src/c/classes/Constraints/SpcStatic.cpp

     
    3636
    3737/*Object virtual functions definitions:*/
    3838Object* SpcStatic::copy() {/*{{{*/
    39        
     39
    4040        SpcStatic* spcstat = new SpcStatic(*this);
    4141
    4242        spcstat->id=this->id;

  • ../trunk-jpl/src/c/classes/IoModel.cpp

     
    156156
    157157        bool autodiff=false;
    158158        bool iscontrol=false;
    159        
     159
    160160        /*First, keep track of the file handle: */
    161161        this->fid=iomodel_handle;
    162162
    … …  
    421421        /*Initialize array detecting whether data[i] is an independent AD mode variable: */
    422422        this->FetchData(&autodiff,"md.autodiff.isautodiff");
    423423        this->FetchData(&iscontrol,"md.inversion.iscontrol");
    424        
     424
    425425        if(trace || (autodiff && !iscontrol)){
    426426
    427427                #ifdef _HAVE_ADOLC_
    … …  
    505505void  IoModel::DeleteData(char*** pstringarray, int numstrings,const char* data_name){/*{{{*/
    506506
    507507        char** stringarray=*pstringarray;
    508        
     508
    509509        if(numstrings){
    510510                for(int i=0;i<numstrings;i++){
    511511                        char* string=stringarray[i];
    … …  
    596596                                case 2:
    597597                                        /*Read the integer and broadcast it to other cpus:*/
    598598                                        if(fread(&integer,sizeof(int),1,this->fid)!=1) _error_("could not read integer ");
    599                                        
     599
    600600                                        /*Convert codes to Enums if needed*/
    601601                                        if(strcmp(record_name,"md.smb.model")==0) integer = IoCodeToEnumSMB(integer);
    602602                                        if(strcmp(record_name,"md.basalforcings.model")==0) integer = IoCodeToEnumBasal(integer);
    … …  
    962962
    963963        /*recover my_rank:*/
    964964        int my_rank=IssmComm::GetRank();
    965        
     965
    966966        /*Set file pointer to beginning of the data: */
    967967        fid=this->SetFilePointerToData(&code,NULL,data_name);
    968968
    … …  
    11331133                        return;
    11341134                }
    11351135        }
    1136          
     1136
    11371137        /*output: */
    11381138        int          M,N;
    11391139        IssmPDouble *matrix = NULL;
    … …  
    18241824                        if(my_rank==0){
    18251825                                /*check we are indeed finding a string, not something else: */
    18261826                                if(codes[i]!=4)_error_("expecting a string for \""<<data_name<<"\"");
    1827                
     1827
    18281828                                /*We have to read a string from disk. First read the dimensions of the string, then the string: */
    18291829                                fsetpos(fid,file_positions+i);
    18301830                                if(fread(&string_size,sizeof(int),1,fid)!=1) _error_("could not read length of string ");
    … …  
    18531853        /*Free ressources:*/
    18541854        xDelete<int>(codes);
    18551855        xDelete<fpos_t>(file_positions);
    1856        
     1856
    18571857        /*Assign output pointers: */
    18581858        *pstrings=strings;
    18591859        *pnumstrings=num_instances;
    … …  
    18741874
    18751875        /*recover my_rank:*/
    18761876        int my_rank=IssmComm::GetRank();
    1877        
     1877
    18781878        /*Get file pointers to beginning of the data (multiple instances of it): */
    18791879        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    18801880
    … …  
    18891889                                code=codes[i];
    18901890
    18911891                                if(code!=2)_error_("expecting an integer for \""<<data_name<<"\"");
    1892                                
     1892
    18931893                                /*We have to read a integer from disk. First read the dimensions of the integer, then the integer: */
    18941894                                fsetpos(fid,file_positions+i);
    18951895                                if(my_rank==0){ 
    … …  
    19021902                        vector[i]=integer;
    19031903                }
    19041904        }
    1905                        
     1905
    19061906        /*Free ressources:*/
    19071907        xDelete<fpos_t>(file_positions);
    19081908        xDelete<int>(codes);
    … …  
    19271927
    19281928        /*recover my_rank:*/
    19291929        int my_rank=IssmComm::GetRank();
    1930        
     1930
    19311931        /*Get file pointers to beginning of the data (multiple instances of it): */
    19321932        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    19331933
    … …  
    19421942                                code=codes[i];
    19431943
    19441944                                if(code!=3)_error_("expecting a double for \""<<data_name<<"\"");
    1945                                
     1945
    19461946                                /*We have to read a double from disk: */
    19471947                                fsetpos(fid,file_positions+i);
    19481948                                if(my_rank==0){ 
    … …  
    19551955                        vector[i]=scalar;
    19561956                }
    19571957        }
    1958                        
     1958
    19591959        /*Free ressources:*/
    19601960        xDelete<fpos_t>(file_positions);
    19611961        xDelete<int>(codes);
    … …  
    19841984
    19851985        /*recover my_rank:*/
    19861986        int my_rank=IssmComm::GetRank();
    1987        
     1987
    19881988        /*Get file pointers to beginning of the data (multiple instances of it): */
    19891989        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    19901990
    … …  
    20142014                        }
    20152015                        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    20162016
    2017 
    20182017                        /*Now allocate matrix: */
    20192018                        if(M*N){
    20202019                                pmatrix=xNew<IssmPDouble>(M*N);
    … …  
    20382037                        }
    20392038                        else
    20402039                                matrix=NULL;
    2041                        
    2042                        
     2040
    20432041                        /*Assign: */
    20442042                        mdims[i]=M;
    20452043                        matrices[i]=matrix;
    … …  
    20462044                        ndims[i]=N;
    20472045                }
    20482046        }
    2049                        
     2047
    20502048        /*Free ressources:*/
    20512049        xDelete<fpos_t>(file_positions);
    20522050        xDelete<int>(codes);
    … …  
    20882086
    20892087        /*recover my_rank:*/
    20902088        int my_rank=IssmComm::GetRank();
    2091        
     2089
    20922090        /*Get file pointers to beginning of the data (multiple instances of it): */
    20932091        file_positions=this->SetFilePointersToData(&codes,NULL,&num_instances,data_name);
    20942092
    … …  
    21182116                        }
    21192117                        ISSM_MPI_Bcast(&N,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    21202118
    2121 
    21222119                        /*Now allocate matrix: */
    21232120                        if(M*N){
    21242121                                pmatrix=xNew<IssmPDouble>(M*N);
    … …  
    21432140                        }
    21442141                        else
    21452142                                integer_matrix=NULL;
    2146                        
    2147                        
     2143
    21482144                        /*Assign: */
    21492145                        mdims[i]=M;
    21502146                        matrices[i]=integer_matrix;
    … …  
    21512147                        ndims[i]=N;
    21522148                }
    21532149        }
    2154                        
     2150
    21552151        /*Free ressources:*/
    21562152        xDelete<fpos_t>(file_positions);
    21572153        xDelete<int>(codes);
    … …  
    23752371                        codes          = xNew<int>(num_instances);
    23762372                        vector_types   = xNew<int>(num_instances);
    23772373                }
    2378        
     2374
    23792375                /*Reset FILE* position to the beginning of the file, and start again, this time saving the data information
    23802376                 * as we find it: */
    23812377                counter=0;
    … …  
    24182414                                codes[counter]        = record_code;
    24192415                                vector_types[counter] = vector_type;
    24202416                                fgetpos(fid,file_positions+counter);
    2421                                
     2417
    24222418                                /*backup and skip over the record, as we have more work to do: */
    24232419                                if(5<=record_code && record_code<=7) fseek(fid,-sizeof(int),SEEK_CUR);
    24242420                                fseek(fid,-sizeof(int),SEEK_CUR);
    24252421                                fseek(fid,-sizeof(int),SEEK_CUR);
    2426                                
     2422
    24272423                                /*increment counter: */
    24282424                                counter++;
    24292425                        }

  • ../trunk-jpl/src/c/classes/Nodalvalue.cpp

     
    8686}
    8787/*}}}*/
    8888IssmDouble Nodalvalue::Response(FemModel* femmodel){/*{{{*/
    89        
     89
    9090         /*output:*/
    9191         IssmDouble value;
    9292

  • ../trunk-jpl/src/c/classes/Misfit.cpp

     
    2222#include "../classes/Inputs/Input.h"
    2323#include "../classes/gauss/Gauss.h"
    2424/*}}}*/
    25                
     25
    2626/*Misfit constructors, destructors :*/
    2727Misfit::Misfit(){/*{{{*/
    2828
    … …  
    4141Misfit::Misfit(char* in_name, int in_definitionenum, int in_model_enum, int in_observation_enum, char* in_timeinterpolation, int in_local, int in_weights_enum){/*{{{*/
    4242
    4343        this->definitionenum=in_definitionenum;
    44        
     44
    4545        this->name              = xNew<char>(strlen(in_name)+1);
    4646        xMemCpy<char>(this->name,in_name,strlen(in_name)+1);
    4747
    4848        this->timeinterpolation = xNew<char>(strlen(in_timeinterpolation)+1);
    4949        xMemCpy<char>(this->timeinterpolation,in_timeinterpolation,strlen(in_timeinterpolation)+1);
    50                                
     50
    5151        this->model_enum=in_model_enum;
    5252        this->observation_enum=in_observation_enum;
    5353        this->weights_enum=in_weights_enum;
    5454        this->local=in_local;
    55        
     55
    5656        this->misfit=0;
    5757        this->lock=0;
    5858}
    … …  
    110110}
    111111/*}}}*/
    112112IssmDouble Misfit::Response(FemModel* femmodel){/*{{{*/
    113                  
     113
    114114         /*diverse: */
    115115         IssmDouble time,starttime,finaltime;
    116116         IssmDouble dt;
    117          
     117
    118118         /*recover time parameters: */
    119119         femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    120120         femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
    … …  
    129129                 IssmDouble area_t=0.;
    130130                 IssmDouble all_area_t;
    131131
    132        
    133132                 /*If we are locked, return time average: */
    134133                 if(this->lock)return misfit/(time-starttime);
    135134
    … …  
    143142                 ISSM_MPI_Allreduce ( (void*)&area_t,(void*)&all_area_t,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
    144143                 area_t=all_area_t;
    145144                 misfit_t=all_misfit_t;
    146                  
     145
    147146                 /*Divide by surface area if not nill!: */
    148147                 if (area_t!=0) misfit_t=misfit_t/area_t;
    149148
    … …  
    163162                 IssmDouble* observation= NULL;
    164163                 IssmDouble* weights= NULL;
    165164                 int msize,osize,wsize;
    166                  
     165
    167166                 /*Are we transient?:*/
    168167                 if (time==0){
    169168                         IssmDouble misfit_t=0.;
    170                          
     169
    171170                         /*get global vectors: */
    172171                         GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
    173172                         GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
    … …  
    189188                         return misfit;
    190189                 }
    191190                 else{
    192                          
     191
    193192                         IssmDouble misfit_t=0.;
    194193                         IssmDouble all_misfit_t=0.;
    195194
    … …  
    200199                         GetVectorFromInputsx(&model,&msize,femmodel,model_enum);
    201200                         GetVectorFromInputsx(&observation,&osize,femmodel,observation_enum);_assert_(msize==osize);
    202201                         GetVectorFromInputsx(&weights,&wsize,femmodel,weights_enum); _assert_(wsize==msize);
    203                          
     202
    204203                         int count=0;
    205204                         for (int i=0;i<msize;i++){
    206205                                 misfit_t += pow(model[i]-observation[i],2)*weights[i];
    … …  
    214213
    215214                         /*Do we lock? i.e. are we at final_time? :*/
    216215                         if(time==finaltime)this->lock=1;
    217                          
     216
    218217                         /*Free ressources:*/
    219218                         xDelete<IssmDouble>(model);
    220219                         xDelete<IssmDouble>(observation);
    … …  
    226225
    227226         } /*}}}*/
    228227         else{ /*global computation: {{{ */
    229                  
     228
    230229                 IssmDouble model, observation;
    231                  
     230
    232231                 /*If we are locked, return time average: */
    233232                 if(this->lock)return misfit/(time-starttime);
    234233
    … …  
    238237                 Element* element=(Element*)femmodel->elements->GetObjectByOffset(0); _assert_(element);
    239238                 Input* input = element->GetInput(observation_enum); _assert_(input);
    240239                 input->GetInputAverage(&observation);
    241                  
     240
    242241                 /*Add this time's contribution to curent misfit: */
    243242                 misfit+=dt*(model-observation);
    244                  
     243
    245244                 /*Do we lock? i.e. are we at final_time? :*/
    246245                 if(time==finaltime)this->lock=1;
    247                  
     246
    248247                 /*What we return is the value of misfit / time: */
    249248                 return misfit/(time-starttime);
    250249         } /*}}}*/

  • ../trunk-jpl/src/c/datastructures/DataSet.cpp

     
    9191
    9292/*Specific methods*/
    9393void  DataSet::Marshall(char** pmarshalled_data, int* pmarshalled_data_size, int marshall_direction){ /*{{{*/
    94        
     94
    9595        vector<Object*>::iterator obj;
    9696        int obj_size=0;
    9797        int obj_enum=0;

  • ../trunk-jpl/src/c/analyses/BalancethicknessAnalysis.cpp

     
    593593        xDelete<int>(responses);
    594594        delete gauss;
    595595
    596 
    597596}/*}}}*/
    598597void           BalancethicknessAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    599598

  • ../trunk-jpl/src/c/analyses/HydrologyDCEfficientAnalysis.cpp

     
    260260                                                basis,1,numnodes,0,
    261261                                                &Ke->values[0],1);
    262262
    263 
    264263                        /*Transfer EPL part*/
    265264                        transfer=GetHydrologyKMatrixTransfer(basalelement);
    266265                        D_scalar=dt*transfer*gauss->weight*Jdet;

  • ../trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

     
    11521152
    11531153        /*Clean up and return*/
    11541154        xDelete<int>(responses);
    1155                          
     1155
    11561156}/*}}}*/
    11571157void           AdjointHorizAnalysis::GradientJBbarFS(Element* element,Vector<IssmDouble>* gradient,int control_index){/*{{{*/
    11581158        /*WARNING: We use SSA as an estimate for now*/
    … …  
    21112111        IssmDouble vx,vy,lambda,mu;
    21122112        IssmDouble *xyz_list= NULL;
    21132113
    2114 
    21152114        /*Fetch number of vertices for this finite element*/
    21162115        int numvertices = basalelement->GetNumberOfVertices();
    21172116
    … …  
    21312130        Input* adjointx_input  = basalelement->GetInput(AdjointxEnum);    _assert_(adjointx_input);
    21322131        Input* adjointy_input  = basalelement->GetInput(AdjointyEnum);    _assert_(adjointy_input);
    21332132
    2134 
    2135 
    21362133        IssmDouble  q_exp;
    21372134        IssmDouble  C_param;
    21382135        IssmDouble  As;
    … …  
    21812178                vmag  = sqrt(vx*vx + vy*vy);
    21822179                Chi   = vmag/(pow(C_param,n)*pow(Neff,n)*As);
    21832180                Gamma = (Chi/(1.+alpha*pow(Chi,q_exp)));
    2184                
     2181
    21852182                Uder =Neff*C_param/(vmag*vmag*n) *
    21862183                        (Gamma-alpha*q_exp*pow(Chi,q_exp-1.)*Gamma*Gamma* pow(Gamma,(1.-n)/n) -
    21872184                         n* pow(Gamma,1./n));
    2188                
     2185
    21892186                /*Build gradient vector (actually -dJ/dD): */
    21902187                for(int i=0;i<numvertices;i++){
    21912188                        ge[i]+=-dalpha2dk*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
    … …  
    23182315                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    23192316        }
    23202317
    2321 
    23222318        /*Fetch number of nodes and dof for this finite element*/
    23232319        int vnumnodes = element->NumberofNodesVelocity();
    23242320        int pnumnodes = element->NumberofNodesPressure();

  • ../trunk-jpl/src/c/analyses/LevelsetAnalysis.cpp

     
    4242                        counter++;
    4343                }
    4444        }
    45        
     45
    4646        iomodel->FetchDataToInput(elements,"md.mask.ice_levelset",MaskIceLevelsetEnum);
    4747        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
    4848        iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
    … …  
    176176        }
    177177
    178178        /*Calving threshold*/
    179        
     179
    180180        /*Fetch number of nodes and dof for this finite element*/
    181181        int numnodes    = basalelement->GetNumberOfNodes();
    182182
    … …  
    325325                                         c[i]=0.; m[i]=0.;
    326326                                 }
    327327                                break;
    328                                
     328
    329329                        case CalvingLevermannEnum:
    330330                                calvingratex_input->GetInputValue(&c[0],gauss);
    331331                                if(dim==2) calvingratey_input->GetInputValue(&c[1],gauss);
    … …  
    356356                                         m[i]=0.;
    357357                                 }
    358358                                break;
    359                        
     359
    360360                        case CalvingHabEnum:
    361361                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    362362                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    … …  
    377377                                         m[i]=0.;
    378378                                 }
    379379                                break;
    380                        
     380
    381381                        case CalvingCrevasseDepthEnum:
    382382                                lsf_slopex_input->GetInputValue(&dlsf[0],gauss);
    383383                                if(dim==2) lsf_slopey_input->GetInputValue(&dlsf[1],gauss);
    … …  
    384384                                meltingrate_input->GetInputValue(&meltingrate,gauss);
    385385
    386386                                if(groundedice<0) meltingrate = 0.;
    387                                
     387
    388388                                norm_dlsf=0.;
    389389                                for(i=0;i<dim;i++) norm_dlsf+=pow(dlsf[i],2);
    390390                                norm_dlsf=sqrt(norm_dlsf);
    … …  
    515515        return Ke;
    516516}/*}}}*/
    517517ElementVector* LevelsetAnalysis::CreatePVector(Element* element){/*{{{*/
    518        
     518
    519519        if(!element->IsOnBase()) return NULL;
    520520        Element* basalelement = element->SpawnBasalElement();
    521521
    … …  
    524524        IssmDouble  Jdet,dt;
    525525        IssmDouble  lsf;
    526526        IssmDouble* xyz_list = NULL;
    527        
     527
    528528        /*Fetch number of nodes and dof for this finite element*/
    529529        int numnodes = basalelement->GetNumberOfNodes();
    530530
    … …  
    531531        /*Initialize Element vector*/
    532532        ElementVector* pe = basalelement->NewElementVector();
    533533        basalelement->FindParam(&dt,TimesteppingTimeStepEnum);
    534        
     534
    535535        if(dt!=0.){
    536536                /*Initialize basis vector*/
    537537                IssmDouble*    basis = xNew<IssmDouble>(numnodes);
    … …  
    622622        // returns distance d of point q to straight going through points s0, s1
    623623        // d=|a x b|/|b|
    624624        // with a=q-s0, b=s1-s0
    625        
     625
    626626        /* Intermediaries */
    627627        const int dim=2;
    628628        int i;
    … …  
    633633                a[i]=q[i]-s0[i];
    634634                b[i]=s1[i]-s0[i];
    635635        }
    636        
     636
    637637        norm_b=0.;
    638638        for(i=0;i<dim;i++)
    639639                norm_b+=b[i]*b[i];
    … …  
    670670        IssmDouble  rho_ice,rho_water;
    671671        IssmDouble  bed,water_depth;
    672672        IssmDouble  levelset;
    673        
     673
    674674        femmodel->parameters->FindParam(&calvinglaw,CalvingLawEnum);
    675675
    676676        if(calvinglaw==CalvingMinthicknessEnum){
    … …  
    702702                        delete gauss;
    703703                }
    704704        }
    705        
     705
    706706        if(calvinglaw==CalvingHabEnum){
    707707
    708708                /*Get the fraction of the flotation thickness at the terminus*/
    709709                femmodel->elements->InputDuplicate(MaskIceLevelsetEnum, DistanceToCalvingfrontEnum);
    710710                femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
    711                
     711
    712712                /*Loop over all elements of this partition*/
    713713                for(int i=0;i<femmodel->elements->Size();i++){
    714714                        Element* element  = xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    715                        
     715
    716716                        rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
    717717                        rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
    718718
    … …  
    743743                        delete gauss;
    744744                }
    745745        }
    746        
     746
    747747        if(calvinglaw==CalvingCrevasseDepthEnum){
    748        
     748
    749749                int                 nflipped,local_nflipped;
    750750                Vector<IssmDouble>* vec_constraint_nodes = NULL;
    751751                IssmDouble* constraint_nodes = NULL;
    752                
     752
    753753                /*Get the DistanceToCalvingfront*/
    754754                femmodel->elements->InputDuplicate(MaskIceLevelsetEnum,DistanceToCalvingfrontEnum);
    755755                femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0,DistanceToCalvingfrontEnum);
    … …  
    841841                        for(int in=0;in<numnodes;in++){
    842842                                gauss->GaussNode(element->GetElementType(),in);
    843843                                Node* node=element->GetNode(in);
    844                                
     844
    845845                                if(constraint_nodes[node->Sid()]>0.){
    846846                                        node->ApplyConstraint(0,+1.);
    847847                                }

  • ../trunk-jpl/src/c/analyses/StressbalanceVerticalAnalysis.cpp

     
    2929        else
    3030         iscoupling = false;
    3131
    32 
    3332        /*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
    3433        if(!iscoupling){
    3534                IoModelToConstraintsx(constraints,iomodel,"md.stressbalance.spcvz",StressbalanceVerticalAnalysisEnum,P1Enum,0);
    … …  
    158157}/*}}}*/
    159158ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixBase(Element* element){/*{{{*/
    160159
    161 
    162160        if(!element->IsOnBase()) return NULL;
    163161
    164162        /*Intermediaries*/
    … …  
    199197}/*}}}*/
    200198ElementMatrix* StressbalanceVerticalAnalysis::CreateKMatrixSurface(Element* element){/*{{{*/
    201199
    202 
    203200        if(!element->IsOnSurface()) return NULL;
    204201
    205202        /*Intermediaries*/

  • ../trunk-jpl/src/c/analyses/HydrologyShaktiAnalysis.cpp

     
    294294                B_input->GetInputValue(&B,gauss);
    295295                n_input->GetInputValue(&n,gauss);
    296296                A=pow(B,-n);
    297                
     297
    298298                /*Compute beta term*/
    299299                if(gap<br)
    300300                 beta = (br-gap)/lr;
    … …  
    323323
    324324        meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
    325325                _assert_(meltrate>0.);
    326        
     326
    327327                for(int i=0;i<numnodes;i++) pe->values[i]+=Jdet*gauss->weight*
    328328                 (
    329329                  meltrate*(1/rho_water-1/rho_ice)
    … …  
    397397
    398398                /*Calculate effective pressure*/
    399399                eff_pressure[i] = rho_ice*g*thickness[i] - rho_water*g*(values[i]-bed[i]);
    400        
     400
    401401                if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
    402402                if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector");
    403403        }
    … …  
    443443        Input* gap_input      = element->GetInput(HydrologyGapHeightEnum); _assert_(gap_input);
    444444        reynolds_input->GetInputAverage(&reynolds);
    445445        gap_input->GetInputAverage(&gap);
    446        
     446
    447447        /*Compute conductivity*/
    448448        IssmDouble conductivity = pow(gap,3)*g/(12.*NU*(1+OMEGA*reynolds));
    449449        _assert_(conductivity>0);
    … …  
    453453}/*}}}*/
    454454void HydrologyShaktiAnalysis::UpdateGapHeight(FemModel* femmodel){/*{{{*/
    455455
    456 
    457456        for(int j=0;j<femmodel->elements->Size();j++){
    458457                Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
    459458                UpdateGapHeight(element);
    … …  
    546545                IssmDouble pressure_ice   = rho_ice*g*thickness;    _assert_(pressure_ice>0.);
    547546                IssmDouble pressure_water = rho_water*g*(head-bed);
    548547                if(pressure_water>pressure_ice) pressure_water = pressure_ice;
    549      
     548
    550549      /* Compute change in sensible heat due to changes in pressure melting point*/
    551550           dpressure_water[0] = rho_water*g*(dh[0] - dbed[0]);
    552551                dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
    … …  
    580579
    581580        /*Add new gap as an input*/
    582581        element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
    583  
     582
    584583        /*Divide by connectivity, add basal flux as an input*/
    585584        q = q/totalweights;
    586585        element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum);

  • ../trunk-jpl/src/c/analyses/EsaAnalysis.cpp

     
    3939        int         nl;
    4040        IssmDouble* love_h=NULL;
    4141        IssmDouble* love_l=NULL;
    42        
     42
    4343        IssmDouble* U_elastic = NULL;
    4444        IssmDouble* U_elastic_local = NULL;
    4545        IssmDouble* H_elastic = NULL;
    … …  
    6868        M=reCast<int,IssmDouble>(180./degacc+1.);
    6969        U_elastic=xNew<IssmDouble>(M);
    7070        H_elastic=xNew<IssmDouble>(M);
    71        
     71
    7272        /*compute combined legendre + love number (elastic green function:*/
    7373        m=DetermineLocalSize(M,IssmComm::GetComm());
    7474        GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    … …  
    9595                        IssmDouble deltalove_U;
    9696
    9797                        deltalove_U = (love_h[n]-love_h[nl-1]);
    98        
     98
    9999                        /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
    100100                        if(n==0){
    101101                                Pn=1;
    … …  
    198198void           EsaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
    199199        /*Default, do nothing*/
    200200        return;
    201        
     201
    202202}/*}}}*/
    203203void           EsaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
    204204        /*Default, do nothing*/

  • ../trunk-jpl/src/c/analyses/SmoothAnalysis.cpp

     
    152152                default: input = element->GetInput(input_enum);
    153153        }
    154154
    155 
    156155        /* Start  looping on the number of gaussian points: */
    157156        Gauss* gauss=element->NewGauss(2);
    158157        for(int ig=gauss->begin();ig<gauss->end();ig++){
    … …  
    161160                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    162161                element->NodalFunctions(basis,gauss);
    163162
    164 
    165163                switch(input_enum){
    166164                        case DrivingStressXEnum:
    167165                        case DrivingStressYEnum:{

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

     
    164164        }
    165165        iomodel->FetchDataToInput(elements,"md.initialization.vx",VxEnum);
    166166        iomodel->FetchDataToInput(elements,"md.initialization.vy",VyEnum);
    167        
     167
    168168        /*Initialize cumdeltalthickness input*/
    169169        InputUpdateFromConstantx(elements,0.,SealevelriseCumDeltathicknessEnum);
    170170
    … …  
    202202                iomodel->FetchDataToInput(elements,"md.mesh.vertexonbase",MeshVertexonbaseEnum);
    203203                iomodel->FetchDataToInput(elements,"md.mesh.vertexonsurface",MeshVertexonsurfaceEnum);
    204204        }
    205        
     205
    206206}/*}}}*/
    207207void MasstransportAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    208208
    … …  
    220220        parameters->AddObject(new IntParam(MasstransportNumRequestedOutputsEnum,numoutputs));
    221221        if(numoutputs)parameters->AddObject(new StringArrayParam(MasstransportRequestedOutputsEnum,requestedoutputs,numoutputs));
    222222        iomodel->DeleteData(&requestedoutputs,numoutputs,"md.masstransport.requested_outputs");
    223        
    224        
    225223
    226224}/*}}}*/
    227225
    … …  
    804802        xDelete<IssmDouble>(phi);
    805803        xDelete<IssmDouble>(sealevel);
    806804        xDelete<IssmDouble>(bed);
    807        
     805
    808806        xDelete<int>(doflist);
    809807        if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
    810808}/*}}}*/

  • ../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

     
    6868        IssmDouble* love_h=NULL;
    6969        IssmDouble* love_k=NULL;
    7070        IssmDouble* love_l=NULL;
    71        
     71
    7272        bool elastic=false;
    7373        IssmDouble* G_elastic = NULL;
    7474        IssmDouble* G_elastic_local = NULL;
    … …  
    131131                H_elastic=xNew<IssmDouble>(M);
    132132                #endif
    133133
    134                
    135134                /*compute combined legendre + love number (elastic green function:*/
    136135                m=DetermineLocalSize(M,IssmComm::GetComm());
    137136                GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm());
    … …  
    167166
    168167                                deltalove_G = (love_k[n]-love_k[nl-1]-love_h[n]+love_h[nl-1]);
    169168                                deltalove_U = (love_h[n]-love_h[nl-1]);
    170                
     169
    171170                                /*compute legendre polynomials: P_n(cos\theta) & d P_n(cos\theta)/ d\theta: */
    172171                                if(n==0){
    173172                                        Pn=1;
    … …  
    229228                xDelete<IssmDouble>(H_elastic);
    230229                xDelete<IssmDouble>(H_elastic_local);
    231230        }
    232        
     231
    233232        /*Transitions: */
    234233        iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,"md.slr.transitions");
    235234        if(transitions){
    … …  
    249248        if(numoutputs)parameters->AddObject(new StringArrayParam(SealevelriseRequestedOutputsEnum,requestedoutputs,numoutputs));
    250249        iomodel->DeleteData(&requestedoutputs,numoutputs,"md.slr.requested_outputs");
    251250
    252 
    253251}/*}}}*/
    254252
    255253/*Finite Element Analysis*/

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

     
    198198                xDelete<IssmDouble>(lset);
    199199        }
    200200
    201 
    202 
    203201        return;
    204202}/*}}}*/
    205203

  • ../trunk-jpl/src/c/analyses/HydrologyShreveAnalysis.cpp

     
    378378        return;
    379379}/*}}}*/
    380380
    381 
    382381/*Needed changes to switch to the Johnson formulation*//*{{{*/
    383382/*All the changes are to be done in the velocity computation.
    384383        The new velocity needs some new parameter that should be introduce in the hydrologyshreve class:

  • ../trunk-jpl/src/c/analyses/BalancevelocityAnalysis.cpp

     
    199199                Ny[i] = -H[i]*Ny[i]/norm;
    200200        }
    201201
    202 
    203202        /* Start  looping on the number of gaussian points: */
    204203        Gauss* gauss=basalelement->NewGauss(2);
    205204        for(int ig=gauss->begin();ig<gauss->end();ig++){

  • ../trunk-jpl/src/c/analyses/HydrologyDCInefficientAnalysis.cpp

     
    213213                return NULL;
    214214        }
    215215
    216 
    217216        /*Intermediaries */
    218217        bool        active_element,isefficientlayer;
    219218        IssmDouble  D_scalar,Jdet,dt;
    … …  
    503502                default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
    504503        }
    505504
    506 
    507505        /*Fetch number of nodes for this finite element*/
    508506        int numnodes = basalelement->GetNumberOfNodes();
    509507

  • ../trunk-jpl/src/c/analyses/SmbAnalysis.cpp

     
    1818        return 1;
    1919}/*}}}*/
    2020void SmbAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
    21        
     21
    2222        int    smb_model;
    2323        bool   isdelta18o,ismungsm,isd18opd,issetpddfac,isprecipscaled,istemperaturescaled;
    24        
     24
    2525        /*Update elements: */
    2626        int counter=0;
    2727        for(int i=0;i<iomodel->numberofelements;i++){
    … …  
    3131                        counter++;
    3232                }
    3333        }
    34        
     34
    3535        /*Figure out smb model: */
    3636        iomodel->FindConstant(&smb_model,"md.smb.model");
    37                        
     37
    3838        switch(smb_model){
    3939                case SMBforcingEnum:
    4040                        iomodel->FetchDataToInput(elements,"md.smb.mass_balance",SmbMassBalanceEnum,0.);
    … …  
    141141                default:
    142142                        _error_("Surface mass balance model "<<EnumToStringx(smb_model)<<" not supported yet");
    143143        }
    144        
    145        
    146144
    147145}/*}}}*/
    148146void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
    … …  
    153151        int     smb_model;
    154152        IssmDouble *temp = NULL;
    155153        int         N,M;
    156        
     154
    157155        parameters->AddObject(iomodel->CopyConstantObject("md.smb.model",SmbEnum));
    158        
     156
    159157        iomodel->FindConstant(&smb_model,"md.smb.model");
    160158        iomodel->FindConstant(&interp,"md.timestepping.interp_forcings");
    161        
     159
    162160        switch(smb_model){
    163161                case SMBforcingEnum:
    164162                        /*Nothing to add to parameters*/
    … …  
    197195                          iomodel->FetchData(&temp,&N,&M,"md.smb.Pfac"); _assert_(N==2);
    198196                          parameters->AddObject(new TransientParam(SmbPfacEnum,&temp[0],&temp[M],interp,M));
    199197                          iomodel->DeleteData(temp,"md.smb.Pfac");
    200                        
     198
    201199                          iomodel->FetchData(&temp,&N,&M,"md.smb.Tdiff"); _assert_(N==2);
    202200                          parameters->AddObject(new TransientParam(SmbTdiffEnum,&temp[0],&temp[M],interp,M));
    203201                          iomodel->DeleteData(temp,"md.smb.Tdiff");
    … …  
    265263
    266264        /*Figure out smb model: */
    267265        femmodel->parameters->FindParam(&smb_model,SmbEnum);
    268        
     266
    269267        /*branch to correct module*/
    270268        switch(smb_model){
    271269                case SMBforcingEnum:

  • ../trunk-jpl/src/c/analyses/DamageEvolutionAnalysis.cpp

     
    125125
    126126        /*Add input*/
    127127        element->AddInput(DamageFEnum,f,element->GetElementType());
    128        
     128
    129129        /*Clean up and return*/
    130130        xDelete<IssmDouble>(f);
    131131}/*}}}*/
    … …  
    169169        Gauss* gauss=element->NewGauss();
    170170        for (int i=0;i<numnodes;i++){
    171171                gauss->GaussNode(element->GetElementType(),i);
    172                
     172
    173173                eps_xx_input->GetInputValue(&eps_xx,gauss);
    174174                eps_xy_input->GetInputValue(&eps_xy,gauss);
    175175                eps_yy_input->GetInputValue(&eps_yy,gauss);
    … …  
    176176                B_input->GetInputValue(&B,gauss);
    177177                n_input->GetInputValue(&n,gauss);
    178178                damage_input->GetInputValue(&damage,gauss);
    179        
     179
    180180                /*Calculate principal effective strain rates*/
    181181                eps1=(eps_xx+eps_yy)/2.+sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
    182182                eps2=(eps_xx+eps_yy)/2.-sqrt(pow((eps_xx-eps_yy)/2.,2)+pow(eps_xy,2));
    … …  
    194194
    195195        /*Add input*/
    196196        element->AddInput(DamageFEnum,f,element->GetElementType());
    197        
     197
    198198        /*Clean up and return*/
    199199        xDelete<IssmDouble>(f);
    200200        delete gauss;
    … …  
    261261        Gauss* gauss=element->NewGauss();
    262262        for (int i=0;i<numnodes;i++){
    263263                gauss->GaussNode(element->GetElementType(),i);
    264                
     264
    265265                damage_input->GetInputValue(&damage,gauss);
    266266                tau_xx_input->GetInputValue(&tau_xx,gauss);
    267267                tau_xy_input->GetInputValue(&tau_xy,gauss);
    … …  
    313313        }
    314314        /*Add input*/
    315315        element->AddInput(DamageFEnum,f,element->GetElementType());
    316        
     316
    317317        /*Clean up and return*/
    318318        xDelete<IssmDouble>(f);
    319319        delete gauss;
    … …  
    374374
    375375                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    376376                element->NodalFunctions(basis,gauss);
    377                
     377
    378378                vx_input->GetInputValue(&vx,gauss);
    379379                vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
    380380                vy_input->GetInputValue(&vy,gauss);
    … …  
    537537                damaged_input = element->GetInput(DamageDEnum); _assert_(damaged_input);
    538538        }
    539539
    540 
    541540        /* Start  looping on the number of gaussian points: */
    542541        Gauss* gauss=element->NewGauss(2);
    543542        for(int ig=gauss->begin();ig<gauss->end();ig++){

  • ../trunk-jpl/src/c/analyses/Balancethickness2Analysis.cpp

     
    77/*Model processing*/
    88void Balancethickness2Analysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
    99
    10 
    1110        int finiteelement = P1Enum;
    1211        IoModelToConstraintsx(constraints,iomodel,"md.balancethickness.spcthickness",Balancethickness2AnalysisEnum,finiteelement);
    1312

  • ../trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp

     
    102102        bool dakota_analysis,ismovingfront,isenthalpy;
    103103        int frictionlaw,basalforcing_model,materialstype;
    104104        int FrictionCoupling;
    105        
     105
    106106        /*Now, is the model 3d? otherwise, do nothing: */
    107107        if(iomodel->domaintype==Domain2DhorizontalEnum)return;
    108108
    … …  
    188188                default:
    189189                        _error_("not supported");
    190190        }
    191        
     191
    192192        /*Friction law variables*/
    193193        switch(frictionlaw){
    194194                case 1:
    … …  
    804804
    805805                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
    806806                element->NodalFunctions(basis,gauss);
    807                
     807
    808808                /*viscous dissipation*/
    809809                element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input);
    810810
    … …  
    822822                        enthalpypicard_input->GetInputDerivativeValue(&d1enthalpypicard[0],xyz_list,gauss);
    823823                        pressure_input->GetInputDerivativeValue(&d1pressure[0],xyz_list,gauss);
    824824                        d2pressure=0.; // for linear elements, 2nd derivative is zero
    825                        
     825
    826826                        d1H_d1P=0.;
    827827                        for(i=0;i<3;i++) d1H_d1P+=d1enthalpypicard[i]*d1pressure[i];
    828828                        d1P2=0.;
    … …  
    10551055        IssmDouble dt;
    10561056        int* basalnodeindices=NULL;
    10571057        Element* element= NULL;
    1058        
     1058
    10591059        femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
    10601060
    10611061        /*depth-integrate the drained water fraction */
    … …  
    13871387        for(int i=0;i<numindices;i++){
    13881388                gauss->GaussNode(element->GetElementType(),indices[i]);
    13891389                gaussup->GaussNode(element->GetElementType(),indicesup[i]);
    1390                
     1390
    13911391                enthalpy_input->GetInputValue(&enthalpy,gauss);
    13921392                enthalpy_input->GetInputValue(&enthalpyup,gaussup);
    13931393                pressure_input->GetInputValue(&pressure,gauss);

  • ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp

     
    144144                workelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
    145145                GetB(B,workelement,xyz_list,gauss,dim);
    146146                GetBprime(Bprime,workelement,xyz_list,gauss,dim);
    147                
     147
    148148                D_scalar=gauss->weight*Jdet;
    149149
    150150                if(extrapolatebydiffusion){
    … …  
    173173                                for(i=0;i<dim;i++)      normal[i]=dlsf[i]/norm_dlsf;
    174174                        else
    175175                                for(i=0;i<dim;i++)      normal[i]=0.;
    176                        
     176
    177177                        for(row=0;row<dim;row++)
    178178                                for(col=0;col<dim;col++)
    179179                                        if(row==col)
    … …  
    336336
    337337        /* Get parameters */
    338338        element->FindParam(&extvar_enum, ExtrapolationVariableEnum);
    339        
     339
    340340        Input* active_input=element->GetInput(IceMaskNodeActivationEnum); _assert_(active_input);
    341341        Input* extvar_input=element->GetInput(extvar_enum); _assert_(extvar_input);
    342342

  • ../trunk-jpl/src/c/main/issm_slr.cpp

     
    2727
    2828        /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
    2929        worldcomm=EnvironmentInit(argc,argv);
    30        
     30
    3131        /*What is my rank?:*/
    3232        ISSM_MPI_Comm_rank(worldcomm,&my_rank);
    3333
    … …  
    3939        rankzeros=xNew<int>(nummodels);
    4040        for(int i=0;i<nummodels;i++){
    4141                char* string=NULL;
    42                
     42
    4343                string=xNew<char>(strlen(argv[5+3*i])+1);
    4444                xMemCpy<char>(string,argv[5+3*i],strlen(argv[5+3*i])+1);
    4545                dirnames[i]=string;
    46                
     46
    4747                string=xNew<char>(strlen(argv[5+3*i+1])+1);
    4848                xMemCpy<char>(string,argv[5+3*i+1],strlen(argv[5+3*i+1])+1);
    4949                modelnames[i]=string;
    … …  
    9393
    9494        /*Initialize femmodel from arguments provided command line: */
    9595        FemModel *femmodel = new FemModel(4,arguments,modelcomm);
    96        
     96
    9797        /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */
    9898        femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(worldcomm,WorldCommEnum));
    9999        femmodel->parameters->AddObject(new IntParam(NumModelsEnum,nummodels));

  • ../trunk-jpl/src/c/main/issm_dakota.cpp

     
    1515
    1616int main(int argc,char **argv){
    1717
    18 
    1918        #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6
    2019
    2120        bool parallel=true;
    … …  
    2827
    2928        /*Initialize MPI: */
    3029        ISSM_MPI_Init(&argc, &argv); // initialize MPI
    31        
     30
    3231        /*Recover file name for dakota input file:*/
    3332        dakota_input_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.in")+2));
    3433        sprintf(dakota_input_file,"%s/%s%s",argv[2],argv[3],".qmu.in");
    35        
     34
    3635        dakota_output_file=xNew<char>((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.out")+2));
    3736        sprintf(dakota_output_file,"%s/%s%s",argv[2],argv[3],".qmu.out");
    3837
    … …  
    5655                        << std::endl;
    5756                Dakota::abort_handler(-1);
    5857        }
    59        
     58
    6059        Dakota::ProblemDescDB& problem_db = env.problem_description_db();
    6160        Dakota::ModelLIter ml_iter;
    6261        size_t model_index = problem_db.get_db_model_node(); // for restoration

  • ../trunk-jpl/src/c/main/issm_ocean.cpp

     
    2121
    2222        /*Initialize environment (MPI, PETSC, MUMPS, etc ...)*/
    2323        worldcomm=EnvironmentInit(argc,argv);
    24        
     24
    2525        /*What is my rank?:*/
    2626        ISSM_MPI_Comm_rank(worldcomm,&my_rank);
    2727        ISSM_MPI_Comm_size(worldcomm,&my_size);
    … …  
    3838        ISSM_MPI_Intercomm_create( modelcomm, 0, worldcomm, my_local_size, 0, &tomitgcmcomm);
    3939
    4040        FemModel *femmodel = new FemModel(argc,argv,modelcomm);
    41        
     41
    4242        /*Now that the models are initialized, keep communicator information in the parameters datasets of each model: */
    4343        femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(worldcomm,WorldCommEnum));
    4444        femmodel->parameters->AddObject(new GenericParam<ISSM_MPI_Comm>(tomitgcmcomm,ToMITgcmCommEnum));

  • ../trunk-jpl/src/c/main/esmfbinders.cpp

     
    22 * \brief: ESMF binders for ISSM. Binders developed initially for the GEOS-5 framework.
    33 */
    44
    5 
    65#include "./issm.h"
    76
    87/*GCM specific declarations:*/
    … …  
    3534
    3635                /*Some specific code here for the binding: */
    3736                femmodel->parameters->SetParam(SMBgcmEnum,SmbEnum); //bypass SMB model, will be provided by GCM!
    38        
     37
    3938                /*Restart file: */
    4039                femmodel->Restart();
    4140
    … …  
    112111                                                {
    113112
    114113                                                IssmDouble surface;
    115                                                
     114
    116115                                                /*Recover surface from the ISSM element: */
    117116                                                Input* surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
    118117                                                surface_input->GetInputAverage(&surface);
    119                        
     118
    120119                                                *(issmoutputs+f*numberofelements+i) = surface;
    121120
    122121                                                }
    … …  
    146145
    147146                /*Output results: */
    148147                OutputResultsx(femmodel);
    149                        
     148
    150149                /*Check point: */
    151150                femmodel->CheckPoint();
    152151

  • ../trunk-jpl/src/c/bamg/BamgQuadtree.cpp

     
    284284                                int xiv = b->v[k]->i.x;
    285285                                int yiv = b->v[k]->i.y;
    286286
    287 
    288287                                int h0 = Norm(xi2,xiv,yi2,yiv);
    289288
    290289                                /*is it smaller than previous value*/

  • ../trunk-jpl/src/c/bamg/Mesh.cpp

     
    500500                int* connectivitysize_1=NULL;
    501501                int  connectivitymax_1=0;
    502502                int verbose=0;
    503        
     503
    504504                /*Check needed pointer*/
    505505                _assert_(bamgmesh);
    506506
    … …  
    10331033                //      ----------------         !
    10341034                //   s0       tt2       s1
    10351035                //-------------------------------
    1036                
     1036
    10371037                /*Intermediaries*/
    10381038                Triangle* tt[3];       //the three triangles
    10391039                long long det3local[3];   //three determinants (integer)
    … …  
    11571157
    11581158                int verbose=0;
    11591159                double lminaniso = 1./ (Max(hminaniso*hminaniso,1e-100));
    1160        
     1160
    11611161                //Get options
    11621162                if(bamgopts) verbose=bamgopts->verbose;
    1163                        
     1163
    11641164                //display info
    11651165                if (verbose > 1)  _printf_("   BoundAnisotropy by " << anisomax << "\n");
    11661166
    … …  
    18531853
    18541854                /*Check pointer*/
    18551855                _assert_(bamgopts);
    1856                
     1856
    18571857                /*Recover options*/
    18581858                verbose=bamgopts->verbose;
    18591859                NbJacobi=bamgopts->nbjacobi;
    … …  
    27272727                long k0=this->RandomNumber(nbv);
    27282728                if (verbose>4) _printf_("      Prime Number = "<<PrimeNumber<<"\n");
    27292729                if (verbose>4) _printf_("      k0 = "<<k0<<"\n");
    2730                
     2730
    27312731                //Build orderedvertices
    27322732                for (i=0; i<nbv; i++){
    27332733                        orderedvertices[i]=&vertices[k0=(k0+PrimeNumber)%nbv];
    … …  
    29522952        /*}}}*/
    29532953        void  Mesh::MaxSubDivision(BamgOpts* bamgopts,double maxsubdiv) {/*{{{*/
    29542954                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/MaxSubDivision)*/
    2955                
     2955
    29562956                /*Intermediaries*/
    29572957                int verbose=0;
    29582958
    … …  
    37393739        /*}}}*/
    37403740        void Mesh::SmoothingVertex(BamgOpts* bamgopts,int nbiter,double omega ) { /*{{{*/
    37413741                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Mesh2.cpp/SmoothingVertex)*/
    3742                
     3742
    37433743                /*Intermediaries*/
    37443744                int verbose=0;
    3745        
     3745
    37463746                /*Get options*/
    37473747                if(bamgopts) verbose=bamgopts->verbose;
    37483748
    … …  
    37843784        /*}}}*/
    37853785        void Mesh::SmoothMetric(BamgOpts* bamgopts,double raisonmax) { /*{{{*/
    37863786                /*Original code from Frederic Hecht <hecht@ann.jussieu.fr> (BAMG v1.01, Metric.cpp/SmoothMetric)*/
    3787                
     3787
    37883788                /*Intermediaries*/
    37893789                int verbose=0;
    37903790
    … …  
    40864086                int verbose=0;
    40874087                int i;
    40884088                Metric M1(1);
    4089        
     4089
    40904090                /*Get options*/
    40914091                if(bamgopts) verbose=bamgopts->verbose;
    40924092

  • ../trunk-jpl/src/c/shared/Elements/Paterson.cpp

     
    1515        /*Switch to celsius from Kelvin: */
    1616        T=temperature-273.15;
    1717
    18 
    1918        if(T<=-45.0){
    2019                B=1.e+8*(-0.000292866376675*pow(T+50.,3)+ 0.011672640664130*pow(T+50.,2)  -0.325004442485481*(T+50.)+  6.524779401948101);
    2120        }

  • ../trunk-jpl/src/c/shared/Elements/ComputeMungsmTemperaturePrecipitation.cpp

     
    1515  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
    1616  IssmDouble tdiffh; 
    1717
    18 
    1918  for (int imonth = 0; imonth<12; imonth++){
    2019    tdiffh = TdiffTime*( TemperaturesLgm[imonth] - TemperaturesPresentday[imonth] );
    2120    monthlytemperaturestmp[imonth] = tdiffh + TemperaturesPresentday[imonth] ;
    … …  
    2827  }
    2928  // printf(" tempera %f\n",monthlytemperaturestmp[1]);
    3029}
    31  

  • ../trunk-jpl/src/c/shared/Elements/EstarComponents.cpp

     
    6565
    6666        /*Get norm of epsprime*/
    6767        epsprime_norm = sqrt(epsprime[0]*epsprime[0] + epsprime[1]*epsprime[1] + epsprime[2]*epsprime[2]);
    68        
     68
    6969        /*Assign output pointers*/
    7070        *pepsprime_norm=epsprime_norm;
    7171}/*}}}*/

  • ../trunk-jpl/src/c/shared/Elements/ComputeD18OTemperaturePrecipitationFromPD.cpp

     
    1010                        bool isPrecipScaled, IssmDouble f, IssmDouble* PrecipitationPresentday,IssmDouble* TemperaturePresentday,
    1111                        IssmDouble* PrecipitationReconstructed,IssmDouble* TemperatureReconstructed, IssmDouble* monthlytemperaturesout,
    1212                        IssmDouble* monthlyprecout){
    13  
     13
    1414  IssmDouble monthlytemperaturestmp[12],monthlyprectmp[12];
    1515  IssmDouble deltaTemp;
    16  
     16
    1717  /* Constants */
    1818  // dpermil = 2.4;/*degrees C per mil*/
    19  
     19
    2020  /*Create Delta Temp to be applied to monthly temps and used in precip scaling*/
    2121  deltaTemp = dpermil * (d018+34.83);   
    22    
     22
    2323  for (int imonth = 0; imonth<12; imonth++){
    24    
     24
    2525         if (isTemperatureScaled)monthlytemperaturestmp[imonth] = TemperaturePresentday[imonth] + deltaTemp;
    2626         else{
    2727                 monthlytemperaturestmp[imonth] = TemperatureReconstructed[imonth];
    … …  
    3030
    3131         if (isPrecipScaled)monthlyprectmp[imonth] = PrecipitationPresentday[imonth]*exp((f/dpermil)*deltaTemp);
    3232         else monthlyprectmp[imonth] = PrecipitationReconstructed[imonth];
    33          
     33
    3434    /*Assign output pointer*/
    3535    *(monthlytemperaturesout+imonth) = monthlytemperaturestmp[imonth];
    3636    *(monthlyprecout+imonth) = monthlyprectmp[imonth];

  • ../trunk-jpl/src/c/shared/Elements/DrainageFunctionWaterfraction.cpp

     
    1818        /*get drainage function value*/
    1919        if((w0==w1)||(w1==w2)||(w0==w2))
    2020                _error_("Error: equal ordinates in DrainageFunctionWaterfraction -> division by zero. Abort");
    21        
     21
    2222        if(waterfraction<=w0)
    2323                Dret=D0;
    2424        else if((waterfraction>w0) && (waterfraction<=w1))

  • ../trunk-jpl/src/c/shared/Elements/PddSurfaceMassBalance.cpp

     
    8282  // seasonal loop
    8383  for (iqj = 0; iqj < 12; iqj++){
    8484    imonth =  ismon[iqj];
    85    
     85
    8686    /*********compute lapse rate ****************/
    8787    st=(s-s0t)/1000.;
    8888    rtlaps=TdiffTime*rlapslgm + (1.-TdiffTime)*rlaps; // lapse rate
    89    
     89
    9090    /*********compute Surface temperature *******/
    9191    monthlytemperatures[imonth]=monthlytemperatures[imonth] - rtlaps *max(st,sealevTime*0.001);
    9292    tstar = monthlytemperatures[imonth];
    … …  
    9898      if (tstar >= -siglimc){ pd = pds[reCast<int,IssmDouble>(tstar/DT + siglim0c)];}}
    9999    else {
    100100      pd = 0.;}
    101    
     101
    102102    /******exp des/elev precip reduction*******/
    103103    sp=(s-s0p)/1000.-deselcut; // deselev effect is wrt chng in topo
    104104    if (sp>0.0){q = exp(-desfac*sp);}
    105105    else {q = 1.0;}
    106    
     106
    107107    qmt= qmt + monthlyprec[imonth]*sconv;  //*sconv to convert in m of ice equivalent per month 
    108108    qmpt= q*monthlyprec[imonth]*sconv;
    109109    qmp= qmp + qmpt;
    … …  
    113113    // ndd(month)=-(tstar-pdd(month)) since ndd+pdd gives expectation of
    114114    // gaussian=T_m, so ndd=-(Tsurf-pdd)
    115115    if (iqj>5 && iqj<9){ Tsum=Tsum+tstar;}
    116    
     116
    117117    if (tstar >= siglim) {pdd = pdd + tstar*deltm;}
    118118    else if (tstar> -siglim){
    119119      pddsig=pdds[reCast<int,IssmDouble>(tstar/DT + siglim0)];
    … …  
    120120      pdd = pdd + pddsig*deltm;
    121121      frzndd = frzndd - (tstar-pddsig)*deltm;}
    122122    else{frzndd = frzndd - tstar*deltm; }
    123    
     123
    124124    /*Assign output pointer*/
    125125    *(monthlytemperatures+imonth) = monthlytemperatures[imonth];
    126126    *(monthlyprec+imonth) = monthlyprec[imonth];     
    … …  
    149149            icefac=pddicefac;
    150150          }
    151151  }
    152  
     152
    153153  /***** determine PDD factors *****/
    154154  if(Tsum< -1.) {
    155155    snwmf=(2.65+snowfac-pddsnowfac0)*0.001;   //  ablation factor for snow per positive degree day.*0.001 to go from mm to m/ppd
    … …  
    165165  }
    166166  snwmf=0.95*snwmf;
    167167  smf=0.95*smf;
    168  
     168
    169169  /*****  compute PDD ablation and refreezing *****/
    170170  pddt = pdd *365;
    171171  snwm = snwmf*pddt;       // snow that could have been melted in a year
    172172  hmx2 = min(h,dfrz);   // refreeze active layer max depth: dfrz
    173  
     173
    174174  if(snwm < saccu) {
    175175    water=prect-saccu + snwm; //water=rain + snowmelt
    176176    //     l 2.2= capillary factor

  • ../trunk-jpl/src/c/shared/Numerics/legendre.cpp

     
    120120          v(1:m,2) = x;
    121121
    122122          for i = 2 : n
    123          
     123
    124124                v(:,i+1) = ( ( 2 * i - 1 ) * x .* v(:,i)   ...
    125125                                        -  (     i - 1 ) *    v(:,i-1) ) ...
    126126                                        /  (     i     );
    127          
     127
    128128          end
    129129          }}}  */
    130130
    … …  
    239239                                           Output, double P_POLYNOMIAL_VALUE[M*(N+1)], the values of the Legendre
    240240                                           polynomials of order 0 through N.
    241241        }}}*/
    242        
     242
    243243        int i,j;
    244244
    245245        if(n<0) return NULL;
    … …  
    253253
    254254        for ( i = 0; i < m; i++ ) {
    255255                for ( j = 2; j <= n; j++ ) {
    256                        
     256
    257257                        v[j+i*(n+1)] = ( ( IssmDouble ) ( 2 * j - 1 ) * x[i] * v[(j-1)+i*(n+1)]
    258258                                        - ( IssmDouble ) (     j - 1 ) *        v[(j-2)+i*(n+1)] )
    259259                                / ( IssmDouble ) (     j     );

  • ../trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp

     
    6767                /*Get current Gradient at xmin=0*/
    6868                _printf0_(" x = "<<setw(9)<<xmin<<" | ");
    6969                fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
    70                
     70
    7171                /*Get f(xmax)*/
    7272                _printf0_(" x = "<<setw(9)<<xmax<<" | ");
    7373                for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
    … …  
    243243                xDelete<IssmDouble>(G);
    244244                J[n]=fxbest;
    245245        }
    246        
     246
    247247        /*return*/
    248248        xDelete<IssmDouble>(X);
    249249        *pJ=J;

  • ../trunk-jpl/src/c/shared/Matrix/MatrixUtils.cpp

     
    6262        else{
    6363                dtemp = &dtemp_static[0];
    6464        }
    65        
     65
    6666        /*  perform the matrix triple product in the order that minimizes the
    6767                 number of multiplies and the temporary space used, noting that
    6868                 (a*b)*c requires ac(b+d) multiplies and ac IssmDoubles, and a*(b*c)
    … …  
    334334        /*Compute determinant*/
    335335        Matrix2x2Determinant(&det,A);
    336336        if (fabs(det) < DBL_EPSILON) _error_("Determinant smaller than machine epsilon");
    337        
     337
    338338        /*Multiplication is faster than divsion, so we multiply by the reciprocal*/
    339339        det_reciprocal = 1./det; 
    340340
    … …  
    426426        *pvx      = vx;
    427427        *pvy      = vy;
    428428
    429 
    430429}/*}}}*/
    431430
    432431void Matrix3x3Determinant(IssmDouble* Adet,IssmDouble* A){/*{{{*/
    … …  
    576575}/*}}}*/
    577576
    578577void newcell(IssmDouble** pcell, IssmDouble newvalue, bool top, int m){  /*{{{*/
    579    
     578
    580579    IssmDouble* cell=NULL;
    581580    IssmDouble* dummy=NULL;
    582    
     581
    583582    /*recover pointer: */
    584583    cell=*pcell;
    585    
     584
    586585    /*reallocate:*/
    587586    dummy=xNew<IssmDouble>(m+1);
    588    
     587
    589588        /*copy data:*/
    590589    if(top){
    591590        dummy[0]=newvalue;
    … …  
    595594        dummy[m]=newvalue;
    596595        for(int i=0;i<m;i++)dummy[i]=cell[i];
    597596    }
    598    
     597
    599598    /*delete and reassign: */
    600599    xDelete<IssmDouble>(cell); cell=dummy;
    601    
     600
    602601    /*assign output pointer:*/
    603602    *pcell=cell;
    604603} /*}}}*/
    … …  
    614613
    615614        /*input: */
    616615        IssmDouble* cell=*pcell;
    617        
     616
    618617        /*output: */
    619618        IssmDouble* newcell=xNew<IssmDouble>(nind);
    620619
    … …  
    621620        for(int i=0;i<nind;i++){
    622621                newcell[i]=cell[indices[i]];
    623622        }
    624        
     623
    625624        /*free allocation:*/
    626625        xDelete<IssmDouble>(cell);
    627626
    … …  
    632631
    633632        /*input: */
    634633        IssmDouble* cell=*pcell;
    635        
     634
    636635        /*output: */
    637636        IssmDouble* newcell=xNew<IssmDouble>(m+1);
    638637
    … …  
    640639        newcell[i]=scale*cell[i];
    641640        newcell[i+1]=scale* cell[i];
    642641        for(int j=i+2;j<m+1;j++)newcell[j]=cell[j-1];
    643        
     642
    644643        /*free allocation:*/
    645644        xDelete<IssmDouble>(cell);
    646645
    … …  
    661660                celllist[x]= va_arg ( arguments, IssmDouble*);
    662661        }
    663662        va_end ( arguments );                 
    664        
     663
    665664        _printf_("Echo of cell: \n");
    666665        for(int i=0;i<m;i++){
    667666                _printf_(i << ": ");

  • ../trunk-jpl/src/c/classes/Materials/Matpar.cpp

     
    5858        lithosphere_density       = 0;
    5959        mantle_shear_modulus      = 0;
    6060        mantle_density            = 0;
    61        
     61
    6262        earth_density             = 0;
    6363
    6464        int nnat,dummy;
    … …  
    337337        matpar->lithosphere_density=this->lithosphere_density;
    338338        matpar->mantle_shear_modulus=this->mantle_shear_modulus;
    339339        matpar->mantle_density=this->mantle_density;
    340        
     340
    341341        matpar->earth_density=this->earth_density;
    342342
    343343        return matpar;
    … …  
    438438        MARSHALLING(lithosphere_density);
    439439        MARSHALLING(mantle_shear_modulus);
    440440        MARSHALLING(mantle_density);
    441        
     441
    442442        //slr:
    443443        MARSHALLING(earth_density);
    444444}

  • ../trunk-jpl/src/c/classes/Loads/Friction.cpp

     
    195195        /*Check to prevent dividing by zero if vmag==0*/
    196196        if(vmag==0. && (s-1.)<0.) alpha_complement=0.;
    197197        else alpha_complement=pow(Neff,r)*pow(vmag,(s-1));
    198        
     198
    199199        _assert_(!xIsNan<IssmDouble>(alpha_complement));
    200200        _assert_(!xIsInf<IssmDouble>(alpha_complement));
    201201
    … …  
    270270        r=drag_q/drag_p;
    271271        s=1./drag_p;
    272272
    273 
    274273        /*Get effective pressure*/
    275274        IssmDouble Neff = EffectivePressure(gauss);
    276275

  • ../trunk-jpl/src/c/classes/FemModel.cpp

     
    242242        char *lockfilename   = NULL;
    243243        bool  waitonlock     = false;
    244244
    245 
    246245        /*Write lock file if requested: */
    247246        this->parameters->FindParam(&waitonlock,SettingsWaitonlockEnum);
    248247        this->parameters->FindParam(&lockfilename,LockFileNameEnum);
    … …  
    43924391        /*Free ressources:*/
    43934392        xDelete<IssmDouble>(RSLg_old);
    43944393
    4395 
    43964394}
    43974395/*}}}*/
    43984396void FemModel::SealevelriseElastic(Vector<IssmDouble>* pUp, Vector<IssmDouble>* pNorth, Vector<IssmDouble>* pEast, Vector<IssmDouble>* pRSLg, IssmDouble* latitude, IssmDouble* longitude, IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz,int loop,int horiz){/*{{{*/
    … …  
    48064804        char** poutputbuffer;
    48074805        size_t* poutputbuffersize;
    48084806
    4809 
    48104807        /*Before we delete the profiler, report statistics for this run: */
    48114808        profiler->Stop(TOTAL);  //final tagging
    48124809        _printf0_("\n");
    … …  
    48684865}/*}}}*/
    48694866#endif
    48704867
    4871 
    48724868#if defined(_HAVE_BAMG_) && !defined(_HAVE_ADOLC_)
    48734869void FemModel::ReMeshBamg(int* pnewnumberofvertices,int* pnewnumberofelements,IssmDouble** pnewx,IssmDouble** pnewy,IssmDouble** pnewz,int** pnewelementslist){/*{{{*/
    48744870
    … …  
    52455241      newx[i] = newxylist[2*i];
    52465242      newy[i] = newxylist[2*i+1];
    52475243   }
    5248        
     5244
    52495245        /*Assign the pointers*/
    52505246   (*pnewelementslist)  = newelementslist; //Matlab indexing
    52515247   (*pnewx)             = newx;

  • ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

     
    27352735                }
    27362736        }
    27372737
    2738 
    27392738        /*Clean-up*/
    27402739        xDelete<IssmDouble>(dbasis);
    27412740}/*}}}*/
Note: See TracBrowser for help on using the repository browser.