Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 17084) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 17085) @@ -304,10 +304,6 @@ ./modules/SpcNodesx/SpcNodesx.h\ ./modules/SpcNodesx/SpcNodesx.cpp\ - ./modules/PositiveDegreeDayx/PositiveDegreeDayx.h\ - ./modules/PositiveDegreeDayx/PositiveDegreeDayx.cpp\ - ./modules/Delta18oParameterizationx/Delta18oParameterizationx.h\ - ./modules/Delta18oParameterizationx/Delta18oParameterizationx.cpp\ - ./modules/SmbGradientsx/SmbGradientsx.h\ - ./modules/SmbGradientsx/SmbGradientsx.cpp\ + ./modules/SurfaceMassBalancex/SurfaceMassBalancex.h\ + ./modules/SurfaceMassBalancex/SurfaceMassBalancex.cpp\ ./modules/MeshPartitionx/MeshPartitionx.h\ ./modules/Reducevectorgtofx/Reducevectorgtofx.cpp\ Index: /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17084) +++ /issm/trunk-jpl/src/c/analyses/MasstransportAnalysis.cpp (revision 17085) @@ -109,4 +109,7 @@ iomodel->FetchDataToInput(elements,SurfaceforcingsBNegEnum); break; + case SMBhenningEnum: + iomodel->FetchDataToInput(elements,SurfaceforcingsMassBalanceEnum,0.); + break; default: _error_("Surface mass balance model "<parameters->FindParam(&save_results,SaveResultsEnum); - femmodel->parameters->FindParam(&smb_model,SurfaceforcingsEnum); femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum); @@ -41,25 +40,8 @@ } - switch(smb_model){ - case SMBEnum: - /*Nothing to be done*/ - break; - case SMBpddEnum: - femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); - if(isdelta18o){ - if(VerboseSolution()) _printf0_(" call Delta18oParametrization module\n"); - Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - } - if(VerboseSolution()) _printf0_(" call positive degree day module\n"); - PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - break; - case SMBgradientsEnum: - if(VerboseSolution())_printf_(" call smb gradients module\n\n"); - SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - break; - default: - _error_("Surface mass balance model "<parameters->FindParam(&smb_model,SurfaceforcingsEnum); + + /*branch to correct module*/ + switch(smb_model){ + case SMBEnum: + /*Nothing to be done*/ + break; + case SMBpddEnum: + femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + if(isdelta18o){ + if(VerboseSolution()) _printf0_(" call Delta18oParametrization module\n"); + Delta18oParameterizationx(femmodel); + } + if(VerboseSolution()) _printf0_(" call positive degree day module\n"); + PositiveDegreeDayx(femmodel); + break; + case SMBgradientsEnum: + if(VerboseSolution())_printf_(" call smb gradients module\n"); + SmbGradientsx(femmodel); + break; + case SMBhenningEnum: + if(VerboseSolution())_printf_(" call smb Henning module\n"); + SmbHenningx(femmodel); + default: + _error_("Surface mass balance model "<elements->Size();i++){ + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + element->SmbGradients(); + } + +}/*}}}*/ +void Delta18oParameterizationx(FemModel* femmodel){/*{{{*/ + + for(int i=0;ielements->Size();i++){ + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + element->Delta18oParameterization(); + } + +}/*}}}*/ +void PositiveDegreeDayx(FemModel* femmodel){/*{{{*/ + + // void PositiveDegreeDayx(hd,vTempsea,vPrec,agd,Tsurf,ni){ + // note "v" prefix means 12 monthly means, ie time dimension + // INPUT parameters: ni: working size of arrays + // INPUT: surface elevation (m): hd(NA) + // monthly mean surface sealevel temperature (degrees C): vTempsea(NA + // ,NTIME) + // monthly mean precip rate (m/yr water equivalent): vPrec(NA,NTIME) + // OUTPUT: mass-balance (m/yr ice): agd(NA) + // mean annual surface temperature (degrees C): Tsurf(NA) + + int i, it, jj, itm; + IssmDouble DT = 0.02, sigfac, snormfac; + IssmDouble signorm = 5.5; // signorm : sigma of the temperature distribution for a normal day + IssmDouble siglim; // sigma limit for the integration which is equal to 2.5 sigmanorm + IssmDouble signormc = signorm - 0.5; // sigma of the temperature distribution for cloudy day + IssmDouble siglimc, siglim0, siglim0c; + IssmDouble tstep, tsint, tint, tstepc; + int NPDMAX = 1504, NPDCMAX = 1454; + //IssmDouble pdds[NPDMAX]={0}; + //IssmDouble pds[NPDCMAX]={0}; + IssmDouble pddt, pd ; // pd : snow/precip fraction, precipitation falling as snow + IssmDouble PDup, PDCUT = 2.0; // PDcut: rain/snow cutoff temperature (C) + IssmDouble tstar; // monthly mean surface temp + + IssmDouble *pdds = NULL; + IssmDouble *pds = NULL; + Element *element = NULL; + + pdds=xNew(NPDMAX+1); + pds=xNew(NPDCMAX+1); + + /* initialize PDD (creation of a lookup table)*/ + tstep = 0.1; + tsint = tstep*0.5; + sigfac = -1.0/(2.0*pow(signorm,2)); + snormfac = 1.0/(signorm*sqrt(2.0*acos(-1.0))); + siglim = 2.5*signorm; + siglimc = 2.5*signormc; + siglim0 = siglim/DT + 0.5; + siglim0c = siglimc/DT + 0.5; + PDup = siglimc+PDCUT; + + itm = reCast((2*siglim/DT + 1.5)); + + if(itm >= NPDMAX) _error_("increase NPDMAX in massBalance.cpp"); + for(it = 0; it < itm; it++){ + // tstar = REAL(it)*DT-siglim; + tstar = it*DT-siglim; + tint = tsint; + pddt = 0.; + for ( jj = 0; jj < 600; jj++){ + if (tint > (tstar+siglim)){break;} + pddt = pddt + tint*exp(sigfac*(pow((tint-tstar),2)))*tstep; + tint = tint+tstep; + } + pdds[it] = pddt*snormfac; + } + pdds[itm+1] = siglim + DT; + + //*********compute PD(T) : snow/precip fraction. precipitation falling as snow + tstepc = 0.1; + tsint = PDCUT-tstepc*0.5; + signormc = signorm - 0.5; + sigfac = -1.0/(2.0*pow(signormc,2)); + snormfac = 1.0/(signormc*sqrt(2.0*acos(-1.0))); + siglimc = 2.5*signormc ; + itm = reCast((PDCUT+2.*siglimc)/DT + 1.5); + if(itm >= NPDCMAX) _error_("increase NPDCMAX in p35com"); + for(it = 0; it < itm; it++ ){ + tstar = it*DT-siglimc; + // tstar = REAL(it)*DT-siglimc; + tint = tsint; // start against upper bound + pd = 0.; + for (jj = 0; jj < 600; jj++){ + if (tint<(tstar-siglimc)) {break;} + pd = pd + exp(sigfac*(pow((tint-tstar),2)))*tstepc; + tint = tint-tstepc; + } + pds[it] = pd*snormfac; // gaussian integral lookup table for snow fraction + } + pds[itm+1] = 0.; + // *******END initialize PDD + + for(i=0;ielements->Size();i++){ + element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + element->PositiveDegreeDay(pdds,pds,signorm); + } + + /*free ressouces: */ + xDelete(pdds); + xDelete(pds); +}/*}}}*/ +void SmbHenningx(FemModel* femmodel){/*{{{*/ + + // void SurfaceMassBalancex(hd,agd,ni){ + // INPUT parameters: ni: working size of arrays + // INPUT: surface elevation (m): hd(NA) + // OUTPUT: mass-balance (m/yr ice): agd(NA) + + for(int i=0;ielements->Size();i++){ + Element* element=dynamic_cast(femmodel->elements->GetObjectByOffset(i)); + } + +}/*}}}*/ Index: /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h =================================================================== --- /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 17085) +++ /issm/trunk-jpl/src/c/modules/SurfaceMassBalancex/SurfaceMassBalancex.h (revision 17085) @@ -0,0 +1,17 @@ +/*!\file: SurfaceMassBalancex.h + * \brief header file for SMB + */ + +#ifndef _SurfaceMassBalancex_H +#define _SurfaceMassBalancex_H + +#include "../../classes/classes.h" + +/* local prototypes: */ +void SurfaceMassBalancex(FemModel* femmodel); +void SmbGradientsx(FemModel* femmodel); +void Delta18oParameterizationx(FemModel* femmodel); +void PositiveDegreeDayx(FemModel* femmodel); +void SmbHenningx(FemModel* femmodel); + +#endif /* _SurfaceMassBalancex_H*/ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 17084) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 17085) @@ -24,5 +24,4 @@ #include "./CreateNodalConstraintsx/CreateNodalConstraintsx.h" #include "./CreateJacobianMatrixx/CreateJacobianMatrixx.h" -#include "./Delta18oParameterizationx/Delta18oParameterizationx.h" #include "./DragCoefficientAbsGradientx/DragCoefficientAbsGradientx.h" #include "./ElementConnectivityx/ElementConnectivityx.h" @@ -75,5 +74,4 @@ #include "./ConstraintsStatex/ConstraintsStatex.h" #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h" -#include "./PositiveDegreeDayx/PositiveDegreeDayx.h" #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h" #include "./Reduceloadx/Reduceloadx.h" @@ -85,5 +83,5 @@ #include "./Scotchx/Scotchx.h" #include "./Shp2Kmlx/Shp2Kmlx.h" -#include "./SmbGradientsx/SmbGradientsx.h" +#include "./SurfaceMassBalancex/SurfaceMassBalancex.h" #include "./Solverx/Solverx.h" #include "./SystemMatricesx/SystemMatricesx.h" Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17084) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 17085) @@ -297,5 +297,5 @@ SurfaceforcingsBPosEnum, SurfaceforcingsBNegEnum, - SMBlapserate, + SMBhenningEnum, /*TO BE REMOVED*/ SurfaceforcingsIspddEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17084) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 17085) @@ -301,4 +301,5 @@ case SurfaceforcingsBPosEnum : return "SurfaceforcingsBPos"; case SurfaceforcingsBNegEnum : return "SurfaceforcingsBNeg"; + case SMBhenningEnum : return "SMBhenning"; case SurfaceforcingsIspddEnum : return "SurfaceforcingsIspdd"; case SurfaceforcingsIssmbgradientsEnum : return "SurfaceforcingsIssmbgradients"; Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17084) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 17085) @@ -307,4 +307,5 @@ else if (strcmp(name,"SurfaceforcingsBPos")==0) return SurfaceforcingsBPosEnum; else if (strcmp(name,"SurfaceforcingsBNeg")==0) return SurfaceforcingsBNegEnum; + else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; else if (strcmp(name,"SurfaceforcingsIspdd")==0) return SurfaceforcingsIspddEnum; else if (strcmp(name,"SurfaceforcingsIssmbgradients")==0) return SurfaceforcingsIssmbgradientsEnum; @@ -382,9 +383,9 @@ else if (strcmp(name,"AdolcParam")==0) return AdolcParamEnum; else if (strcmp(name,"BoolInput")==0) return BoolInputEnum; - else if (strcmp(name,"BoolParam")==0) return BoolParamEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"Contour")==0) return ContourEnum; + if (strcmp(name,"BoolParam")==0) return BoolParamEnum; + else if (strcmp(name,"Contour")==0) return ContourEnum; else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"DatasetInput")==0) return DatasetInputEnum; @@ -505,9 +506,9 @@ else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else if (strcmp(name,"ThicknessAbsGradient")==0) return ThicknessAbsGradientEnum; - else if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; + if (strcmp(name,"ThicknessAlongGradient")==0) return ThicknessAlongGradientEnum; + else if (strcmp(name,"ThicknessAcrossGradient")==0) return ThicknessAcrossGradientEnum; else if (strcmp(name,"IntMatParam")==0) return IntMatParamEnum; else if (strcmp(name,"RheologyBbarAbsGradient")==0) return RheologyBbarAbsGradientEnum; @@ -628,9 +629,9 @@ else if (strcmp(name,"TriangleInterp")==0) return TriangleInterpEnum; else if (strcmp(name,"BilinearInterp")==0) return BilinearInterpEnum; - else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"XY")==0) return XYEnum; + if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; + else if (strcmp(name,"XY")==0) return XYEnum; else if (strcmp(name,"XYZ")==0) return XYZEnum; else if (strcmp(name,"Dense")==0) return DenseEnum;