Changeset 21317


Ignore:
Timestamp:
10/26/16 17:09:55 (8 years ago)
Author:
Eric.Larour
Message:

CHG: some fixes on EarthMassBalance.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/branches/trunk-larour-NatGeoScience2016/src/c/cores/sealevelrise_core.cpp

    r21313 r21317  
    8080                femmodel->parameters->FindParam(&earthid,EarthIdEnum);
    8181                femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
    82                 if((modelid==EarthIdEnum) & !ismasstransport) EarthMassTransport(femmodel);
     82                if((modelid==earthid) && (ismasstransport==0)) EarthMassTransport(femmodel);
    8383        }
    8484
    … …  
    9696                InputUpdateFromVectorx(femmodel,Sg,SealevelEnum,VertexSIdEnum);
    9797
    98                 /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */
    99                 /*Initialize:*/
    100                 U_radial = new Vector<IssmDouble>(gsize);
    101                 U_north = new Vector<IssmDouble>(gsize);
    102                 U_east = new Vector<IssmDouble>(gsize);
    103                 Sg_absolute = new Vector<IssmDouble>(gsize);
    104                
    105                 /*call the geodetic main modlule:*/
    106                 femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz);
    107 
    108                 /*compute: absolute sea level change = relative sea level change + vertical motion*/
    109                 Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1);
    110                
    111                 /*get results into elements:*/
    112                 InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum);    // radial displacement
    113                 InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum);     // north motion
    114                 InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum);              // east motion
    115                 InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum);
    116                
     98                int geodetic=0;
     99                if (geodetic){
     100
     101                        /*compute other geodetic signatures, such as absolute sea level chagne, components of 3-D crustal motion: */
     102                        /*Initialize:*/
     103                        U_radial = new Vector<IssmDouble>(gsize);
     104                        U_north = new Vector<IssmDouble>(gsize);
     105                        U_east = new Vector<IssmDouble>(gsize);
     106                        Sg_absolute = new Vector<IssmDouble>(gsize);
     107
     108                        /*call the geodetic main modlule:*/
     109                        femmodel->SealevelriseGeodetic(U_radial,U_north,U_east,Sg,latitude,longitude,radius,xx,yy,zz);
     110
     111                        /*compute: absolute sea level change = relative sea level change + vertical motion*/
     112                        Sg->Copy(Sg_absolute); Sg_absolute->AXPY(U_radial,1);
     113
     114                        /*get results into elements:*/
     115                        InputUpdateFromVectorx(femmodel,U_radial,SealevelUmotionEnum,VertexSIdEnum);    // radial displacement
     116                        InputUpdateFromVectorx(femmodel,U_north,SealevelNmotionEnum,VertexSIdEnum);     // north motion
     117                        InputUpdateFromVectorx(femmodel,U_east,SealevelEmotionEnum,VertexSIdEnum);              // east motion
     118                        InputUpdateFromVectorx(femmodel,Sg_absolute,SealevelAbsoluteEnum,VertexSIdEnum);
     119                        delete U_radial;
     120                        delete U_north;
     121                        delete U_east;
     122                        delete Sg_absolute;
     123                }       
     124
    117125                if(save_results){
    118126                        if(VerboseSolution()) _printf0_("   saving results\n");
    … …  
    120128                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
    121129                }
    122                        
     130
    123131                if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
    124132
    … …  
    126134                delete Sg;
    127135                delete Sg_eustatic;
    128                 delete U_radial;
    129                 delete U_north;
    130                 delete U_east;
    131                 delete Sg_absolute;
    132136                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
    133137        }
    … …  
    209213                 *First, build the global delta thickness vector in the earth model: */
    210214                nv=femmodel->vertices->NumberOfVertices();
    211                 forcingglobal= new Vector<IssmDouble>(nv);
     215                GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
    212216
    213217                /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
    … …  
    377381        Vector<IssmDouble> *deltathickness    = NULL;
    378382        int nv;
     383       
     384        if(VerboseSolution()) _printf0_("      computing earth mass transport\n");
    379385
    380386        /*No mass transport module was called, so we are just going to retrieve the geometry thickness
    … …  
    390396        deltathickness = new Vector<IssmDouble>(nv);
    391397        newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1);
    392                
     398
    393399        /*plug into elements:*/
    394400        InputUpdateFromVectorx(femmodel,deltathickness,SealevelriseDeltathicknessEnum,VertexSIdEnum);
Note: See TracChangeset for help on using the changeset viewer.