Skip to content

CHG: add GEMB mapping option for precip downscaling and dlw downscaling with temperature when lapse rate is set to zero #69

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Jun 17, 2025
Merged

CHG: add GEMB mapping option for precip downscaling and dlw downscaling with temperature when lapse rate is set to zero #69

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
009ab62
CHG: add GEMB mapping option for precip downscaling and dlw downscali…
NJSchlegel Jun 17, 2025
07d5cea
CHG: update class description for the dlw down lapse rate for mapping…
NJSchlegel Jun 17, 2025
e127a1a
CHG: Merge pull request #7 from NJSchlegel/GEMB_mapping
NJSchlegel Jun 17, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/c/analyses/SmbAnalysis.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -429,6 +429,7 @@ void SmbAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int s
parameters->AddObject(iomodel->CopyConstantObject("md.smb.isconstrainsurfaceT",SmbIsconstrainsurfaceTEnum));
parameters->AddObject(iomodel->CopyConstantObject("md.smb.isdeltaLWup",SmbIsdeltaLWupEnum));
parameters->AddObject(iomodel->CopyConstantObject("md.smb.ismappedforcing",SmbIsmappedforcingEnum));
parameters->AddObject(iomodel->CopyConstantObject("md.smb.isprecipforcingremapped",SmbIsprecipforcingremappedEnum));
parameters->AddObject(iomodel->CopyConstantObject("md.smb.InitDensityScaling",SmbInitDensityScalingEnum));
parameters->AddObject(iomodel->CopyConstantObject("md.smb.ThermoDeltaTScaling",SmbThermoDeltaTScalingEnum));
parameters->AddObject(iomodel->CopyConstantObject("md.smb.adThresh",SmbAdThreshEnum));
Expand Down
38 changes: 32 additions & 6 deletions src/c/classes/Elements/Element.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -5578,6 +5578,9 @@ void Element::SmbGemb(IssmDouble timeinputs, int count, int steps){/*{{{*/
pAir_input->GetInputValue(&pAir,gauss); // screen level air pressure [Pa]
//_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
} else {
IssmDouble Dtol = 1e-11;
IssmDouble gravity;
parameters->FindParam(&gravity,ConstantsGEnum);

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

bool isprecipmap=true;
parameters->FindParam(&isprecipmap,SmbIsprecipforcingremappedEnum);

parameters->FindParam(&tlapse,SmbLapseTaValueEnum);
parameters->FindParam(&dlwlapse,SmbLapsedlwrfValueEnum);

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

//Variables for downscaling
IssmDouble taparam, dlwrfparam, rhparam, eaparam, pparam;
IssmDouble pscale = -8400;
IssmDouble taparam, dlwrfparam, rhparam, eaparam, pparam, prparam;
parameters->FindParam(&taparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbTaParamEnum);
parameters->FindParam(&dlwrfparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbDlwrfParamEnum);
parameters->FindParam(&eaparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbEAirParamEnum);
parameters->FindParam(&pparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbPAirParamEnum);
parameters->FindParam(&prparam, Mappedpoint-1, timeinputs, timestepping, dt, SmbPParamEnum);

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

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

//Hold reltive humidity constant and calculte new saturation vapor pressure
if ( (fabs(dlwlapse) > Dtol) || (fabs(tlapse) > Dtol)){
IssmDouble Rg = 8.314; // gas constant (J mol-1 K-1)
IssmDouble dAir = 0.0;
// calculated air density [kg/m3]
// dAir = 0.029 * pAir /(R * Ta);
dAir=0.029 * pparam /(Rg * Ta);
pAir=pparam-gravity*dAir*(currentsurface - elevation[Mappedpoint-1]);
}
else pAir=pparam;

//Hold relative humidity constant and calculte new saturation vapor pressure
//https://cran.r-project.org/web/packages/humidity/vignettes/humidity-measures.html
//es over ice calculation
//Ding et al., 2019 after Bolton, 1980
//ea37 = rh37*100*6.112.*exp((17.67*(t237-273.15))./(t237-29.65));
rhparam=eaparam/6.112/exp((17.67*(taparam-273.15))/(taparam-29.65));
eAir=rhparam*6.112*exp((17.67*(Ta-273.15))/(Ta-29.65));

if (isprecipmap && (eaparam>0)){
P=prparam*eAir/eaparam;
}
else P=prparam;

xDelete<IssmDouble>(elevation);
}
/*}}}*/
Expand Down
19 changes: 10 additions & 9 deletions src/c/shared/Enum/Enum.vim
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -607,6 +607,7 @@ syn keyword cConstant SmbIsgraingrowthEnum
syn keyword cConstant SmbIsmappedforcingEnum
syn keyword cConstant SmbIsmeltEnum
syn keyword cConstant SmbIsmungsmEnum
syn keyword cConstant SmbIsprecipforcingremappedEnum
syn keyword cConstant SmbIsprecipscaledEnum
syn keyword cConstant SmbIssetpddfacEnum
syn keyword cConstant SmbIsshortwaveEnum
Expand Down Expand Up @@ -3794,16 +3795,15 @@ syn keyword cType Cfsurfacelogvel
syn keyword cType Cfsurfacesquare
syn keyword cType Cfsurfacesquaretransient
syn keyword cType Channel
syn keyword cType classes
syn keyword cType Constraint
syn keyword cType Constraints
syn keyword cType Contour
syn keyword cType Contours
syn keyword cType ControlInput
syn keyword cType ControlParam
syn keyword cType Covertree
syn keyword cType DatasetInput
syn keyword cType DataSetParam
syn keyword cType DatasetInput
syn keyword cType Definition
syn keyword cType DependentObject
syn keyword cType DoubleInput
Expand All @@ -3816,20 +3816,19 @@ syn keyword cType Element
syn keyword cType ElementHook
syn keyword cType ElementInput
syn keyword cType ElementMatrix
syn keyword cType Elements
syn keyword cType ElementVector
syn keyword cType Elements
syn keyword cType ExponentialVariogram
syn keyword cType ExternalResult
syn keyword cType FemModel
syn keyword cType FileParam
syn keyword cType Friction
syn keyword cType Gauss
syn keyword cType GaussianVariogram
syn keyword cType gaussobjects
syn keyword cType GaussPenta
syn keyword cType GaussSeg
syn keyword cType GaussTetra
syn keyword cType GaussTria
syn keyword cType GaussianVariogram
syn keyword cType GenericExternalResult
syn keyword cType GenericOption
syn keyword cType GenericParam
Expand All @@ -3846,7 +3845,6 @@ syn keyword cType IntVecParam
syn keyword cType IoModel
syn keyword cType IssmDirectApplicInterface
syn keyword cType IssmParallelDirectApplicInterface
syn keyword cType krigingobjects
syn keyword cType Load
syn keyword cType Loads
syn keyword cType Masscon
Expand All @@ -3857,7 +3855,6 @@ syn keyword cType Materials
syn keyword cType Matestar
syn keyword cType Matice
syn keyword cType Matlitho
syn keyword cType matrixobjects
syn keyword cType MatrixParam
syn keyword cType Misfit
syn keyword cType Moulin
Expand All @@ -3884,13 +3881,13 @@ syn keyword cType Quadtree
syn keyword cType Radar
syn keyword cType Regionaloutput
syn keyword cType Results
syn keyword cType Riftfront
syn keyword cType RiftStruct
syn keyword cType Riftfront
syn keyword cType SealevelGeometry
syn keyword cType Seg
syn keyword cType SegInput
syn keyword cType Segment
syn keyword cType SegRef
syn keyword cType Segment
syn keyword cType SpcDynamic
syn keyword cType SpcStatic
syn keyword cType SpcTransient
Expand All @@ -3911,6 +3908,10 @@ syn keyword cType Variogram
syn keyword cType VectorParam
syn keyword cType Vertex
syn keyword cType Vertices
syn keyword cType classes
syn keyword cType gaussobjects
syn keyword cType krigingobjects
syn keyword cType matrixobjects
syn keyword cType AdjointBalancethickness2Analysis
syn keyword cType AdjointBalancethicknessAnalysis
syn keyword cType AdjointHorizAnalysis
Expand Down
1 change: 1 addition & 0 deletions src/c/shared/Enum/EnumDefinitions.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -601,6 +601,7 @@ enum definitions{
SmbIsmappedforcingEnum,
SmbIsmeltEnum,
SmbIsmungsmEnum,
SmbIsprecipforcingremappedEnum,
SmbIsprecipscaledEnum,
SmbIssetpddfacEnum,
SmbIsshortwaveEnum,
Expand Down
1 change: 1 addition & 0 deletions src/c/shared/Enum/EnumToStringx.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -609,6 +609,7 @@ const char* EnumToStringx(int en){
case SmbIsmappedforcingEnum : return "SmbIsmappedforcing";
case SmbIsmeltEnum : return "SmbIsmelt";
case SmbIsmungsmEnum : return "SmbIsmungsm";
case SmbIsprecipforcingremappedEnum : return "SmbIsprecipforcingremapped";
case SmbIsprecipscaledEnum : return "SmbIsprecipscaled";
case SmbIssetpddfacEnum : return "SmbIssetpddfac";
case SmbIsshortwaveEnum : return "SmbIsshortwave";
Expand Down
1 change: 1 addition & 0 deletions src/c/shared/Enum/Enumjl.vim
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -600,6 +600,7 @@ syn keyword juliaConstC SmbIsgraingrowthEnum
syn keyword juliaConstC SmbIsmappedforcingEnum
syn keyword juliaConstC SmbIsmeltEnum
syn keyword juliaConstC SmbIsmungsmEnum
syn keyword juliaConstC SmbIsprecipforcingremappedEnum
syn keyword juliaConstC SmbIsprecipscaledEnum
syn keyword juliaConstC SmbIssetpddfacEnum
syn keyword juliaConstC SmbIsshortwaveEnum
Expand Down
Loading