Changeset 22185


Ignore:
Timestamp:
10/23/17 15:05:08 (7 years ago)
Author:
aleahsommers
Message:

NEW: added degree of channelization as output; CHG: removed lower cap on water pressure to allow negative pressure, turned off geothermal flux and frictional heat in melt rate for SHMIP simulations

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

Legend:

Unmodified
Added
Removed

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

    r21526 r22185  
    313313                PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
    314314
    315                 meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
     315                meltrate = 1/latentheat*(rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); //Taking out geothermal and frictional heat for SHMIP
     316   //   meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
    316317                _assert_(meltrate>0.);
    317318       
    … …  
    380381
    381382                /*Make sure that negative pressure is not allowed*/
    382       if(values[i]<bed[i]){
    383                         values[i] = bed[i];
    384                 }
     383  //    if(values[i]<bed[i]){
     384        //              values[i] = bed[i];
     385        //      }
    385386
    386387                /*Under-relaxation*/
    … …  
    466467   IssmDouble  dpressure_water[2],dbed[2],PMPheat;
    467468        IssmDouble q = 0.;
     469   IssmDouble channelization = 0.;
    468470
    469471        /*Retrieve all inputs and parameters*/
    … …  
    542544                dpressure_water[1] = rho_water*g*(dh[1] - dbed[1]);
    543545                PMPheat=-CT*CW*conductivity*(dh[0]*dpressure_water[0]+dh[1]*dpressure_water[1]);
    544        
    545                 meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
     546
     547                meltrate = 1/latentheat*(rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat); // Taking out geothermal and frictional heat for SHMIP
     548//              meltrate = 1/latentheat*(G+frictionheat+rho_water*g*conductivity*(dh[0]*dh[0]+dh[1]*dh[1])-PMPheat);
    546549                _assert_(meltrate>0.);
    547550
    … …  
    555558                /* Compute basal water flux */
    556559      q += gauss->weight*Jdet*(conductivity*sqrt(dh[0]*dh[0]+dh[1]*dh[1]));
     560
     561                /* Compute "degree of channelization" (ratio of melt opening to opening by sliding) */
     562                channelization += gauss->weight*Jdet*(meltrate/rho_ice/(beta*sqrt(vx*vx+vy*vy)));
    557563        }
    558564
    … …  
    565571        if(newgap>thickness)
    566572         newgap = thickness;
    567          
     573
    568574        /*Add new gap as an input*/
    569575        element->AddInput(HydrologyGapHeightEnum,&newgap,P0Enum);
    … …  
    573579        element->AddInput(HydrologyBasalFluxEnum,&q,P0Enum);
    574580
     581        /* Divide by connectivity, add degree of channelization as an input */
     582        channelization = channelization/totalweights;
     583        element->AddInput(DegreeOfChannelizationEnum,&channelization,P0Enum);
     584
    575585        /*Clean up and return*/
    576586        xDelete<IssmDouble>(xyz_list);

  • issm/trunk-jpl/src/c/cores/hydrology_core.cpp

    r21463 r22185  
    8989                if(save_results){
    9090                        if(VerboseSolution()) _printf0_("   saving results \n");
    91                         int outputs[4] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum,HydrologyBasalFluxEnum};
    92                         femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4);
     91                        int outputs[5] = {HydrologyHeadEnum,HydrologyGapHeightEnum,EffectivePressureEnum,HydrologyBasalFluxEnum,DegreeOfChannelizationEnum};
     92                        femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],5);
    9393                       
    9494                        /*unload results*/

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

    r22181 r22185  
    185185        HydrologyBasalFluxEnum,
    186186        HydrologyStorageEnum,
     187        DegreeOfChannelizationEnum,
    187188        InversionControlParametersEnum,
    188189        InversionControlScalingFactorsEnum,

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

    r22181 r22185  
    191191                case HydrologyBasalFluxEnum : return "HydrologyBasalFlux";
    192192                case HydrologyStorageEnum : return "HydrologyStorage";
     193                case DegreeOfChannelizationEnum : return "DegreeOfChannelization";
    193194                case InversionControlParametersEnum : return "InversionControlParameters";
    194195                case InversionControlScalingFactorsEnum : return "InversionControlScalingFactors";

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

    r22181 r22185  
    194194              else if (strcmp(name,"HydrologyBasalFlux")==0) return HydrologyBasalFluxEnum;
    195195              else if (strcmp(name,"HydrologyStorage")==0) return HydrologyStorageEnum;
     196              else if (strcmp(name,"DegreeOfChannelization")==0) return DegreeOfChannelizationEnum;
    196197              else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
    197198              else if (strcmp(name,"InversionControlScalingFactors")==0) return InversionControlScalingFactorsEnum;
    … …  
    259260              else if (strcmp(name,"DamageEvolutionRequestedOutputs")==0) return DamageEvolutionRequestedOutputsEnum;
    260261              else if (strcmp(name,"Damage")==0) return DamageEnum;
    261               else if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
    262262         else stage=3;
    263263   }
    264264   if(stage==3){
    265               if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
     265              if (strcmp(name,"NewDamage")==0) return NewDamageEnum;
     266              else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum;
    266267              else if (strcmp(name,"CalvingLaw")==0) return CalvingLawEnum;
    267268              else if (strcmp(name,"CalvingCalvingrate")==0) return CalvingCalvingrateEnum;
    … …  
    270271              else if (strcmp(name,"CalvingDev")==0) return CalvingDevEnum;
    271272              else if (strcmp(name,"CalvingMinthickness")==0) return CalvingMinthicknessEnum;
    272          else stage=3;
    273    }
    274    if(stage==3){
    275               if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
     273              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    276274              else if (strcmp(name,"CalvingCrevasseDepth")==0) return CalvingCrevasseDepthEnum;
    277275              else if (strcmp(name,"SurfaceCrevasse")==0) return SurfaceCrevasseEnum;
    … …  
    385383              else if (strcmp(name,"TransientIsesa")==0) return TransientIsesaEnum;
    386384              else if (strcmp(name,"TransientIsdamageevolution")==0) return TransientIsdamageevolutionEnum;
    387               else if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
     385         else stage=4;
     386   }
     387   if(stage==4){
     388              if (strcmp(name,"TransientIshydrology")==0) return TransientIshydrologyEnum;
    388389              else if (strcmp(name,"TransientIsmovingfront")==0) return TransientIsmovingfrontEnum;
    389390              else if (strcmp(name,"TransientIsslr")==0) return TransientIsslrEnum;
    … …  
    393394              else if (strcmp(name,"BalancethicknessDiffusionCoefficient")==0) return BalancethicknessDiffusionCoefficientEnum;
    394395              else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum;
    395          else stage=4;
    396    }
    397    if(stage==4){
    398               if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
     396              else if (strcmp(name,"BalancethicknessD0")==0) return BalancethicknessD0Enum;
    399397              else if (strcmp(name,"Smb")==0) return SmbEnum;
    400398              else if (strcmp(name,"SmbAnalysis")==0) return SmbAnalysisEnum;
    … …  
    508506              else if (strcmp(name,"Adjointx")==0) return AdjointxEnum;
    509507              else if (strcmp(name,"Adjointy")==0) return AdjointyEnum;
    510               else if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
     508         else stage=5;
     509   }
     510   if(stage==5){
     511              if (strcmp(name,"Adjointz")==0) return AdjointzEnum;
    511512              else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum;
    512513              else if (strcmp(name,"BedSlopeX")==0) return BedSlopeXEnum;
    … …  
    516517              else if (strcmp(name,"Internal")==0) return InternalEnum;
    517518              else if (strcmp(name,"MassFlux")==0) return MassFluxEnum;
    518          else stage=5;
    519    }
    520    if(stage==5){
    521               if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
     519              else if (strcmp(name,"MeltingOffset")==0) return MeltingOffsetEnum;
    522520              else if (strcmp(name,"Misfit")==0) return MisfitEnum;
    523521              else if (strcmp(name,"Pressure")==0) return PressureEnum;
    … …  
    631629              else if (strcmp(name,"Outputdefinition12")==0) return Outputdefinition12Enum;
    632630              else if (strcmp(name,"Outputdefinition13")==0) return Outputdefinition13Enum;
    633               else if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
     631         else stage=6;
     632   }
     633   if(stage==6){
     634              if (strcmp(name,"Outputdefinition14")==0) return Outputdefinition14Enum;
    634635              else if (strcmp(name,"Outputdefinition15")==0) return Outputdefinition15Enum;
    635636              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    … …  
    639640              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    640641              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    641          else stage=6;
    642    }
    643    if(stage==6){
    644               if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
     642              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    645643              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    646644              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    … …  
    754752              else if (strcmp(name,"Gset")==0) return GsetEnum;
    755753              else if (strcmp(name,"Index")==0) return IndexEnum;
    756               else if (strcmp(name,"Indexed")==0) return IndexedEnum;
     754         else stage=7;
     755   }
     756   if(stage==7){
     757              if (strcmp(name,"Indexed")==0) return IndexedEnum;
    757758              else if (strcmp(name,"Intersect")==0) return IntersectEnum;
    758759              else if (strcmp(name,"Nodal")==0) return NodalEnum;
    … …  
    762763              else if (strcmp(name,"OutputFilePointer")==0) return OutputFilePointerEnum;
    763764              else if (strcmp(name,"ToolkitsFileName")==0) return ToolkitsFileNameEnum;
    764          else stage=7;
    765    }
    766    if(stage==7){
    767               if (strcmp(name,"RootPath")==0) return RootPathEnum;
     765              else if (strcmp(name,"RootPath")==0) return RootPathEnum;
    768766              else if (strcmp(name,"OutputFileName")==0) return OutputFileNameEnum;
    769          else stage=7;
    770    }
    771    if(stage==7){
    772               if (strcmp(name,"InputFileName")==0) return InputFileNameEnum;
     767              else if (strcmp(name,"InputFileName")==0) return InputFileNameEnum;
    773768              else if (strcmp(name,"LockFileName")==0) return LockFileNameEnum;
    774769              else if (strcmp(name,"RestartFileName")==0) return RestartFileNameEnum;
    … …  
    880875              else if (strcmp(name,"AmrThicknessErrorResolution")==0) return AmrThicknessErrorResolutionEnum;
    881876              else if (strcmp(name,"AmrThicknessErrorThreshold")==0) return AmrThicknessErrorThresholdEnum;
    882               else if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum;
     877         else stage=8;
     878   }
     879   if(stage==8){
     880              if (strcmp(name,"AmrDeviatoricErrorResolution")==0) return AmrDeviatoricErrorResolutionEnum;
    883881              else if (strcmp(name,"AmrDeviatoricErrorThreshold")==0) return AmrDeviatoricErrorThresholdEnum;
    884882              else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum;
    … …  
    888886              else if (strcmp(name,"XYZ")==0) return XYZEnum;
    889887              else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
    890          else stage=8;
    891    }
    892    if(stage==8){
    893               if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
     888              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    894889              else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    895890              else if (strcmp(name,"Contour")==0) return ContourEnum;
    … …  
    932927              else if (strcmp(name,"NodeSId")==0) return NodeSIdEnum;
    933928              else if (strcmp(name,"ElementSId")==0) return ElementSIdEnum;
    934          else stage=8;
    935    }
    936    if(stage==8){
    937               if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
     929              else if (strcmp(name,"VectorParam")==0) return VectorParamEnum;
    938930              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    939931              else if (strcmp(name,"Segment")==0) return SegmentEnum;
    … …  
    1006998              else if (strcmp(name,"EsaAnalysis")==0) return EsaAnalysisEnum;
    1007999              else if (strcmp(name,"LoveSolution")==0) return LoveSolutionEnum;
    1008               else if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum;
     1000         else stage=9;
     1001   }
     1002   if(stage==9){
     1003              if (strcmp(name,"LoveAnalysis")==0) return LoveAnalysisEnum;
    10091004              else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum;
    10101005              else if (strcmp(name,"ExtrapolationAnalysis")==0) return ExtrapolationAnalysisEnum;
    … …  
    10141009              else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum;
    10151010              else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum;
    1016          else stage=9;
    1017    }
    1018    if(stage==9){
    1019               if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
     1011              else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum;
    10201012              else if (strcmp(name,"L1L2Approximation")==0) return L1L2ApproximationEnum;
    10211013              else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum;
    … …  
    10571049              else if (strcmp(name,"P1bubblecondensed")==0) return P1bubblecondensedEnum;
    10581050              else if (strcmp(name,"P2")==0) return P2Enum;
    1059          else stage=9;
    1060    }
    1061    if(stage==9){
    1062               if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
     1051              else if (strcmp(name,"P2bubble")==0) return P2bubbleEnum;
    10631052              else if (strcmp(name,"P2bubblecondensed")==0) return P2bubblecondensedEnum;
    10641053              else if (strcmp(name,"P2xP1")==0) return P2xP1Enum;
Note: See TracChangeset for help on using the changeset viewer.