Skip to content

Commit 0b2e4e7

Browse files
NJSchlegelNJSchlegel
authored
CHG: Merge pull request #69 from NJSchlegel/main
CHG: add GEMB mapping option for precip downscaling and dlw downscaling with temperature when lapse rate is set to zero
2 parents d11622c + e127a1a commit 0b2e4e7

File tree

16 files changed

+360
-73
lines changed

16 files changed

+360
-73
lines changed

‎src/c/analyses/SmbAnalysis.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -429,6 +429,7 @@ void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
429429
parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
430430
parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdeltaLWup",SmbIsdeltaLWupEnum));
431431
parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismappedforcing",SmbIsmappedforcingEnum));
432+
parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipforcingremapped",SmbIsprecipforcingremappedEnum));
432433
parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
433434
parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
434435
parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));

‎src/c/classes/Elements/Element.cpp‎

Lines changed: 32 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -5578,6 +5578,9 @@ void Element::SmbGemb(IssmDouble timeinputs, int count, int steps){/*{{{*/
55785578
pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa]
55795579
//_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
55805580
} else {
5581+
IssmDouble Dtol = 1e-11;
5582+
IssmDouble gravity;
5583+
parameters->FindParam(&gravity,ConstantsGEnum);
55815584

55825585
int timestepping;
55835586
IssmDouble dt;
@@ -5587,38 +5590,61 @@ void Element::SmbGemb(IssmDouble timeinputs, int count, int steps){/*{{{*/
55875590
Input *currentsurface_input = this->GetInput(SurfaceEnum); _assert_(currentsurface_input);
55885591
currentsurface_input->GetInputAverage(&currentsurface);
55895592

5593+
bool isprecipmap=true;
5594+
parameters->FindParam(&isprecipmap,SmbIsprecipforcingremappedEnum);
5595+
55905596
parameters->FindParam(&tlapse,SmbLapseTaValueEnum);
55915597
parameters->FindParam(&dlwlapse,SmbLapsedlwrfValueEnum);
55925598

55935599
IssmDouble* elevation = NULL;
55945600
parameters->FindParam(&elevation,&N,SmbMappedforcingelevationEnum); _assert_(elevation);
55955601

55965602
//Variables for downscaling
5597-
IssmDouble taparam, dlwrfparam, rhparam, eaparam, pparam;
5598-
IssmDouble pscale = -8400;
5603+
IssmDouble taparam, dlwrfparam, rhparam, eaparam, pparam, prparam;
55995604
parameters->FindParam(&taparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbTaParamEnum);
56005605
parameters->FindParam(&dlwrfparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbDlwrfParamEnum);
56015606
parameters->FindParam(&eaparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbEAirParamEnum);
56025607
parameters->FindParam(&pparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbPAirParamEnum);
5608+
parameters->FindParam(&prparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbPParamEnum);
56035609

56045610
//Variables not downscaled
56055611
parameters->FindParam(&V, Mappedpoint-1, timeinputs, timestepping, dt, SmbVParamEnum);
56065612
parameters->FindParam(&dsw, Mappedpoint-1, timeinputs, timestepping, dt, SmbDswrfParamEnum);
56075613
parameters->FindParam(&dswdiff, Mappedpoint-1, timeinputs, timestepping, dt, SmbDswdiffrfParamEnum);
5608-
parameters->FindParam(&P, Mappedpoint-1, timeinputs, timestepping, dt, SmbPParamEnum);
56095614

56105615
Ta = taparam + (currentsurface - elevation[Mappedpoint-1])*tlapse;
5611-
dlw = fmax(dlwrfparam + (currentsurface - elevation[Mappedpoint-1])*dlwlapse,0.0);
5612-
if (dlwlapse!=0 || tlapse!=0) pAir=pparam*exp((currentsurface - elevation[Mappedpoint-1])/pscale);
5616+
if (fabs(dlwlapse) > Dtol) dlw = fmax(dlwrfparam + (currentsurface - elevation[Mappedpoint-1])*dlwlapse,0.0);
5617+
else{
5618+
//adjust downward longwave, holding emissivity equal (Glover et al, 1999)
5619+
IssmDouble SB = 5.67e-8; // Stefan-Boltzmann constant (W m-2 K-4)
5620+
IssmDouble effe = 1.;
5621+
effe = dlwrfparam/(SB * pow(taparam,4.0));
5622+
dlw = fmax(effe*SB*pow(Ta,4.0),0.0);
5623+
}
56135624

5614-
//Hold reltive humidity constant and calculte new saturation vapor pressure
5625+
if ( (fabs(dlwlapse) > Dtol) || (fabs(tlapse) > Dtol)){
5626+
IssmDouble Rg = 8.314; // gas constant (J mol-1 K-1)
5627+
IssmDouble dAir = 0.0;
5628+
// calculated air density [kg/m3]
5629+
// dAir = 0.029 * pAir /(R * Ta);
5630+
dAir=0.029 * pparam /(Rg * Ta);
5631+
pAir=pparam-gravity*dAir*(currentsurface - elevation[Mappedpoint-1]);
5632+
}
5633+
else pAir=pparam;
5634+
5635+
//Hold relative humidity constant and calculte new saturation vapor pressure
56155636
//https://cran.r-project.org/web/packages/humidity/vignettes/humidity-measures.html
56165637
//es over ice calculation
56175638
//Ding et al., 2019 after Bolton, 1980
56185639
//ea37 = rh37*100*6.112.*exp((17.67*(t237-273.15))./(t237-29.65));
56195640
rhparam=eaparam/6.112/exp((17.67*(taparam-273.15))/(taparam-29.65));
56205641
eAir=rhparam*6.112*exp((17.67*(Ta-273.15))/(Ta-29.65));
56215642

5643+
if (isprecipmap && (eaparam>0)){
5644+
P=prparam*eAir/eaparam;
5645+
}
5646+
else P=prparam;
5647+
56225648
xDelete<IssmDouble>(elevation);
56235649
}
56245650
/*}}}*/

‎src/c/shared/Enum/Enum.vim‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -607,6 +607,7 @@ syn keyword cConstant SmbIsgraingrowthEnum
607607
syn keyword cConstant SmbIsmappedforcingEnum
608608
syn keyword cConstant SmbIsmeltEnum
609609
syn keyword cConstant SmbIsmungsmEnum
610+
syn keyword cConstant SmbIsprecipforcingremappedEnum
610611
syn keyword cConstant SmbIsprecipscaledEnum
611612
syn keyword cConstant SmbIssetpddfacEnum
612613
syn keyword cConstant SmbIsshortwaveEnum
@@ -3794,16 +3795,15 @@ syn keyword cType Cfsurfacelogvel
37943795
syn keyword cType Cfsurfacesquare
37953796
syn keyword cType Cfsurfacesquaretransient
37963797
syn keyword cType Channel
3797-
syn keyword cType classes
37983798
syn keyword cType Constraint
37993799
syn keyword cType Constraints
38003800
syn keyword cType Contour
38013801
syn keyword cType Contours
38023802
syn keyword cType ControlInput
38033803
syn keyword cType ControlParam
38043804
syn keyword cType Covertree
3805-
syn keyword cType DatasetInput
38063805
syn keyword cType DataSetParam
3806+
syn keyword cType DatasetInput
38073807
syn keyword cType Definition
38083808
syn keyword cType DependentObject
38093809
syn keyword cType DoubleInput
@@ -3816,20 +3816,19 @@ syn keyword cType Element
38163816
syn keyword cType ElementHook
38173817
syn keyword cType ElementInput
38183818
syn keyword cType ElementMatrix
3819-
syn keyword cType Elements
38203819
syn keyword cType ElementVector
3820+
syn keyword cType Elements
38213821
syn keyword cType ExponentialVariogram
38223822
syn keyword cType ExternalResult
38233823
syn keyword cType FemModel
38243824
syn keyword cType FileParam
38253825
syn keyword cType Friction
38263826
syn keyword cType Gauss
3827-
syn keyword cType GaussianVariogram
3828-
syn keyword cType gaussobjects
38293827
syn keyword cType GaussPenta
38303828
syn keyword cType GaussSeg
38313829
syn keyword cType GaussTetra
38323830
syn keyword cType GaussTria
3831+
syn keyword cType GaussianVariogram
38333832
syn keyword cType GenericExternalResult
38343833
syn keyword cType GenericOption
38353834
syn keyword cType GenericParam
@@ -3846,7 +3845,6 @@ syn keyword cType IntVecParam
38463845
syn keyword cType IoModel
38473846
syn keyword cType IssmDirectApplicInterface
38483847
syn keyword cType IssmParallelDirectApplicInterface
3849-
syn keyword cType krigingobjects
38503848
syn keyword cType Load
38513849
syn keyword cType Loads
38523850
syn keyword cType Masscon
@@ -3857,7 +3855,6 @@ syn keyword cType Materials
38573855
syn keyword cType Matestar
38583856
syn keyword cType Matice
38593857
syn keyword cType Matlitho
3860-
syn keyword cType matrixobjects
38613858
syn keyword cType MatrixParam
38623859
syn keyword cType Misfit
38633860
syn keyword cType Moulin
@@ -3884,13 +3881,13 @@ syn keyword cType Quadtree
38843881
syn keyword cType Radar
38853882
syn keyword cType Regionaloutput
38863883
syn keyword cType Results
3887-
syn keyword cType Riftfront
38883884
syn keyword cType RiftStruct
3885+
syn keyword cType Riftfront
38893886
syn keyword cType SealevelGeometry
38903887
syn keyword cType Seg
38913888
syn keyword cType SegInput
3892-
syn keyword cType Segment
38933889
syn keyword cType SegRef
3890+
syn keyword cType Segment
38943891
syn keyword cType SpcDynamic
38953892
syn keyword cType SpcStatic
38963893
syn keyword cType SpcTransient
@@ -3911,6 +3908,10 @@ syn keyword cType Variogram
39113908
syn keyword cType VectorParam
39123909
syn keyword cType Vertex
39133910
syn keyword cType Vertices
3911+
syn keyword cType classes
3912+
syn keyword cType gaussobjects
3913+
syn keyword cType krigingobjects
3914+
syn keyword cType matrixobjects
39143915
syn keyword cType AdjointBalancethickness2Analysis
39153916
syn keyword cType AdjointBalancethicknessAnalysis
39163917
syn keyword cType AdjointHorizAnalysis

‎src/c/shared/Enum/EnumDefinitions.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -601,6 +601,7 @@ enum definitions{
601601
SmbIsmappedforcingEnum,
602602
SmbIsmeltEnum,
603603
SmbIsmungsmEnum,
604+
SmbIsprecipforcingremappedEnum,
604605
SmbIsprecipscaledEnum,
605606
SmbIssetpddfacEnum,
606607
SmbIsshortwaveEnum,

‎src/c/shared/Enum/EnumToStringx.cpp‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -609,6 +609,7 @@ const char* EnumToStringx(int en){
609609
case SmbIsmappedforcingEnum : return "SmbIsmappedforcing";
610610
case SmbIsmeltEnum : return "SmbIsmelt";
611611
case SmbIsmungsmEnum : return "SmbIsmungsm";
612+
case SmbIsprecipforcingremappedEnum : return "SmbIsprecipforcingremapped";
612613
case SmbIsprecipscaledEnum : return "SmbIsprecipscaled";
613614
case SmbIssetpddfacEnum : return "SmbIssetpddfac";
614615
case SmbIsshortwaveEnum : return "SmbIsshortwave";

‎src/c/shared/Enum/Enumjl.vim‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -600,6 +600,7 @@ syn keyword juliaConstC SmbIsgraingrowthEnum
600600
syn keyword juliaConstC SmbIsmappedforcingEnum
601601
syn keyword juliaConstC SmbIsmeltEnum
602602
syn keyword juliaConstC SmbIsmungsmEnum
603+
syn keyword juliaConstC SmbIsprecipforcingremappedEnum
603604
syn keyword juliaConstC SmbIsprecipscaledEnum
604605
syn keyword juliaConstC SmbIssetpddfacEnum
605606
syn keyword juliaConstC SmbIsshortwaveEnum

0 commit comments

Comments
 (0)