source: issm/oecreview/Archive/24307-24683/ISSM-24682-24683.diff

Last change on this file was 24684, checked in by Mathieu Morlighem, 5 years ago

CHG: added new review
File size: 35.6 KB

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

     
    714714        SmbDziniEnum,
    715715        SmbEAirEnum,
    716716        SmbECEnum,
     717        SmbECDtEnum,
    717718        SmbECiniEnum,
    718719        SmbElaEnum,
    719720        SmbEvaporationEnum,
    … …  
    733734        SmbMeanSHFEnum,
    734735        SmbMeanULWEnum,
    735736        SmbMeltEnum,
     737        SmbMInitnum,
    736738        SmbMonthlytemperaturesEnum,
    737739        SmbNetLWEnum,
    738740        SmbNetSWEnum,
    … …  
    771773        SmbVmeanEnum,
    772774        SmbVzEnum,
    773775        SmbWEnum,
     776        SmbWAddEnum,
    774777        SmbWiniEnum,
    775778        SmbZMaxEnum,
    776779        SmbZMinEnum,

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

     
    719719                case SmbDziniEnum : return "SmbDzini";
    720720                case SmbEAirEnum : return "SmbEAir";
    721721                case SmbECEnum : return "SmbEC";
     722                case SmbECDtEnum : return "SmbECDt";
    722723                case SmbECiniEnum : return "SmbECini";
    723724                case SmbElaEnum : return "SmbEla";
    724725                case SmbEvaporationEnum : return "SmbEvaporation";
    … …  
    776777                case SmbVmeanEnum : return "SmbVmean";
    777778                case SmbVzEnum : return "SmbVz";
    778779                case SmbWEnum : return "SmbW";
     780                case SmbWAddEnum : return "SmbWAdd";
    779781                case SmbWiniEnum : return "SmbWini";
    780782                case SmbZMaxEnum : return "SmbZMax";
    781783                case SmbZMinEnum : return "SmbZMin";

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

     
    734734              else if (strcmp(name,"SmbDzini")==0) return SmbDziniEnum;
    735735              else if (strcmp(name,"SmbEAir")==0) return SmbEAirEnum;
    736736              else if (strcmp(name,"SmbEC")==0) return SmbECEnum;
     737              else if (strcmp(name,"SmbECDt")==0) return SmbECDtEnum;
    737738              else if (strcmp(name,"SmbECini")==0) return SmbECiniEnum;
    738739              else if (strcmp(name,"SmbEla")==0) return SmbElaEnum;
    739740              else if (strcmp(name,"SmbEvaporation")==0) return SmbEvaporationEnum;
    … …  
    750751              else if (strcmp(name,"SmbMassBalanceSubstep")==0) return SmbMassBalanceSubstepEnum;
    751752              else if (strcmp(name,"SmbMassBalanceTransient")==0) return SmbMassBalanceTransientEnum;
    752753              else if (strcmp(name,"SmbMeanLHF")==0) return SmbMeanLHFEnum;
    753               else if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
    754754         else stage=7;
    755755   }
    756756   if(stage==7){
    757               if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
     757              if (strcmp(name,"SmbMeanSHF")==0) return SmbMeanSHFEnum;
     758              else if (strcmp(name,"SmbMeanULW")==0) return SmbMeanULWEnum;
    758759              else if (strcmp(name,"SmbMelt")==0) return SmbMeltEnum;
    759760              else if (strcmp(name,"SmbMonthlytemperatures")==0) return SmbMonthlytemperaturesEnum;
    760761              else if (strcmp(name,"SmbNetLW")==0) return SmbNetLWEnum;
    … …  
    794795              else if (strcmp(name,"SmbVmean")==0) return SmbVmeanEnum;
    795796              else if (strcmp(name,"SmbVz")==0) return SmbVzEnum;
    796797              else if (strcmp(name,"SmbW")==0) return SmbWEnum;
     798              else if (strcmp(name,"SmbWAdd")==0) return SmbWAddEnum;
    797799              else if (strcmp(name,"SmbWini")==0) return SmbWiniEnum;
    798800              else if (strcmp(name,"SmbZMax")==0) return SmbZMaxEnum;
    799801              else if (strcmp(name,"SmbZMin")==0) return SmbZMinEnum;
    … …  
    872874              else if (strcmp(name,"Outputdefinition16")==0) return Outputdefinition16Enum;
    873875              else if (strcmp(name,"Outputdefinition17")==0) return Outputdefinition17Enum;
    874876              else if (strcmp(name,"Outputdefinition18")==0) return Outputdefinition18Enum;
    875               else if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
    876               else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
    877877         else stage=8;
    878878   }
    879879   if(stage==8){
    880               if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
     880              if (strcmp(name,"Outputdefinition19")==0) return Outputdefinition19Enum;
     881              else if (strcmp(name,"Outputdefinition20")==0) return Outputdefinition20Enum;
     882              else if (strcmp(name,"Outputdefinition21")==0) return Outputdefinition21Enum;
    881883              else if (strcmp(name,"Outputdefinition22")==0) return Outputdefinition22Enum;
    882884              else if (strcmp(name,"Outputdefinition23")==0) return Outputdefinition23Enum;
    883885              else if (strcmp(name,"Outputdefinition24")==0) return Outputdefinition24Enum;
    … …  
    995997              else if (strcmp(name,"BoolInput")==0) return BoolInputEnum;
    996998              else if (strcmp(name,"BoolInput2")==0) return BoolInput2Enum;
    997999              else if (strcmp(name,"IntInput2")==0) return IntInput2Enum;
    998               else if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
    999               else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
    10001000         else stage=9;
    10011001   }
    10021002   if(stage==9){
    1003               if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
     1003              if (strcmp(name,"BoolParam")==0) return BoolParamEnum;
     1004              else if (strcmp(name,"Boundary")==0) return BoundaryEnum;
     1005              else if (strcmp(name,"BuddJacka")==0) return BuddJackaEnum;
    10041006              else if (strcmp(name,"CalvingDev2")==0) return CalvingDev2Enum;
    10051007              else if (strcmp(name,"CalvingHab")==0) return CalvingHabEnum;
    10061008              else if (strcmp(name,"CalvingLevermann")==0) return CalvingLevermannEnum;
    … …  
    11181120              else if (strcmp(name,"IceVolumeScaled")==0) return IceVolumeScaledEnum;
    11191121              else if (strcmp(name,"IcefrontMassFlux")==0) return IcefrontMassFluxEnum;
    11201122              else if (strcmp(name,"IcefrontMassFluxLevelset")==0) return IcefrontMassFluxLevelsetEnum;
    1121               else if (strcmp(name,"Incremental")==0) return IncrementalEnum;
    1122               else if (strcmp(name,"Indexed")==0) return IndexedEnum;
    11231123         else stage=10;
    11241124   }
    11251125   if(stage==10){
    1126               if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
     1126              if (strcmp(name,"Incremental")==0) return IncrementalEnum;
     1127              else if (strcmp(name,"Indexed")==0) return IndexedEnum;
     1128              else if (strcmp(name,"IntExternalResult")==0) return IntExternalResultEnum;
    11271129              else if (strcmp(name,"IntInput")==0) return IntInputEnum;
    11281130              else if (strcmp(name,"ElementInput2")==0) return ElementInput2Enum;
    11291131              else if (strcmp(name,"SegInput2")==0) return SegInput2Enum;
    … …  
    12411243              else if (strcmp(name,"ProfilingSolutionTime")==0) return ProfilingSolutionTimeEnum;
    12421244              else if (strcmp(name,"Regionaloutput")==0) return RegionaloutputEnum;
    12431245              else if (strcmp(name,"Regular")==0) return RegularEnum;
    1244               else if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
    1245               else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
    12461246         else stage=11;
    12471247   }
    12481248   if(stage==11){
    1249               if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
     1249              if (strcmp(name,"RecoveryAnalysis")==0) return RecoveryAnalysisEnum;
     1250              else if (strcmp(name,"Riftfront")==0) return RiftfrontEnum;
     1251              else if (strcmp(name,"SIAApproximation")==0) return SIAApproximationEnum;
    12501252              else if (strcmp(name,"SMBcomponents")==0) return SMBcomponentsEnum;
    12511253              else if (strcmp(name,"SMBd18opdd")==0) return SMBd18opddEnum;
    12521254              else if (strcmp(name,"SMBforcing")==0) return SMBforcingEnum;

  • ../trunk-jpl/src/c/shared/Enum/Enum.vim

     
    148148syn keyword cConstant FlowequationIsSSAEnum
    149149syn keyword cConstant FrictionCouplingEnum
    150150syn keyword cConstant FrictionDeltaEnum
     151syn keyword cConstant FrictionEffectivePressureLimitEnum
    151152syn keyword cConstant FrictionFEnum
    152153syn keyword cConstant FrictionGammaEnum
    153154syn keyword cConstant FrictionLawEnum
    … …  
    716717syn keyword cConstant SmbDziniEnum
    717718syn keyword cConstant SmbEAirEnum
    718719syn keyword cConstant SmbECEnum
     720syn keyword cConstant SmbECDtEnum
    719721syn keyword cConstant SmbECiniEnum
    720722syn keyword cConstant SmbElaEnum
    721723syn keyword cConstant SmbEvaporationEnum
    … …  
    773775syn keyword cConstant SmbVmeanEnum
    774776syn keyword cConstant SmbVzEnum
    775777syn keyword cConstant SmbWEnum
     778syn keyword cConstant SmbWAddEnum
    776779syn keyword cConstant SmbWiniEnum
    777780syn keyword cConstant SmbZMaxEnum
    778781syn keyword cConstant SmbZMinEnum
    … …  
    13341337syn keyword cType Cfsurfacelogvel
    13351338syn keyword cType Cfsurfacesquare
    13361339syn keyword cType Channel
    1337 syn keyword cType classes
    13381340syn keyword cType Constraint
    13391341syn keyword cType Constraints
    13401342syn keyword cType Contour
    … …  
    13411343syn keyword cType Contours
    13421344syn keyword cType ControlInput2
    13431345syn keyword cType Covertree
     1346syn keyword cType DataSetParam
    13441347syn keyword cType DatasetInput2
    1345 syn keyword cType DataSetParam
    13461348syn keyword cType Definition
    13471349syn keyword cType DependentObject
    13481350syn keyword cType DoubleMatArrayParam
    … …  
    13541356syn keyword cType ElementHook
    13551357syn keyword cType ElementInput2
    13561358syn keyword cType ElementMatrix
     1359syn keyword cType ElementVector
    13571360syn keyword cType Elements
    1358 syn keyword cType ElementVector
    13591361syn keyword cType ExponentialVariogram
    13601362syn keyword cType ExternalResult
    13611363syn keyword cType FemModel
    … …  
    13621364syn keyword cType FileParam
    13631365syn keyword cType Friction
    13641366syn keyword cType Gauss
    1365 syn keyword cType GaussianVariogram
    1366 syn keyword cType gaussobjects
    13671367syn keyword cType GaussPenta
    13681368syn keyword cType GaussSeg
    13691369syn keyword cType GaussTetra
    13701370syn keyword cType GaussTria
     1371syn keyword cType GaussianVariogram
    13711372syn keyword cType GenericExternalResult
    13721373syn keyword cType GenericOption
    13731374syn keyword cType GenericParam
    … …  
    13821383syn keyword cType IoModel
    13831384syn keyword cType IssmDirectApplicInterface
    13841385syn keyword cType IssmParallelDirectApplicInterface
    1385 syn keyword cType krigingobjects
    13861386syn keyword cType Load
    13871387syn keyword cType Loads
    13881388syn keyword cType Masscon
    … …  
    13931393syn keyword cType Matestar
    13941394syn keyword cType Matice
    13951395syn keyword cType Matlitho
    1396 syn keyword cType matrixobjects
    13971396syn keyword cType MatrixParam
    13981397syn keyword cType Misfit
    13991398syn keyword cType Moulin
    … …  
    14061405syn keyword cType Observation
    14071406syn keyword cType Observations
    14081407syn keyword cType Option
     1408syn keyword cType OptionUtilities
    14091409syn keyword cType Options
    1410 syn keyword cType OptionUtilities
    14111410syn keyword cType Param
    14121411syn keyword cType Parameters
    14131412syn keyword cType Pengrid
    … …  
    14211420syn keyword cType Radar
    14221421syn keyword cType Regionaloutput
    14231422syn keyword cType Results
     1423syn keyword cType RiftStruct
    14241424syn keyword cType Riftfront
    1425 syn keyword cType RiftStruct
    14261425syn keyword cType Seg
    14271426syn keyword cType SegInput2
     1427syn keyword cType SegRef
    14281428syn keyword cType Segment
    1429 syn keyword cType SegRef
    14301429syn keyword cType SpcDynamic
    14311430syn keyword cType SpcStatic
    14321431syn keyword cType SpcTransient
    … …  
    14451444syn keyword cType VectorParam
    14461445syn keyword cType Vertex
    14471446syn keyword cType Vertices
     1447syn keyword cType classes
     1448syn keyword cType gaussobjects
     1449syn keyword cType krigingobjects
     1450syn keyword cType matrixobjects
    14481451syn keyword cType AdjointBalancethickness2Analysis
    14491452syn keyword cType AdjointBalancethicknessAnalysis
    14501453syn keyword cType AdjointHorizAnalysis
    … …  
    14631466syn keyword cType ExtrudeFromTopAnalysis
    14641467syn keyword cType FreeSurfaceBaseAnalysis
    14651468syn keyword cType FreeSurfaceTopAnalysis
     1469syn keyword cType GLheightadvectionAnalysis
    14661470syn keyword cType GiaIvinsAnalysis
    1467 syn keyword cType GLheightadvectionAnalysis
    14681471syn keyword cType HydrologyDCEfficientAnalysis
    14691472syn keyword cType HydrologyDCInefficientAnalysis
    14701473syn keyword cType HydrologyGlaDSAnalysis

  • ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp

     
    55#include "./SurfaceMassBalancex.h"
    66#include "../../shared/shared.h"
    77#include "../../toolkits/toolkits.h"
     8#include "../modules.h"
     9#include "../../classes/Inputs2/TransientInput2.h"
    810
    911const double Pi = 3.141592653589793;
    1012const double CtoK = 273.15;             // Kelvin to Celcius conversion/ice melt. point T in K
    … …  
    3133
    3234void Gembx(FemModel* femmodel){  /*{{{*/
    3335
    34         for(int i=0;i<femmodel->elements->Size();i++){
    35         Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
    36         element->SmbGemb();
     36        int        count=0;
     37        IssmDouble time,dt,finaltime,starttime;
     38        IssmDouble timeclim=0.0;
     39        IssmDouble t,smb_dt;
     40   IssmDouble delta;
     41        bool       isclimatology=false;
     42
     43        femmodel->parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
     44   femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
     45        femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
     46        femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
     47   femmodel->parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
     48        femmodel->parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
     49
     50        //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
     51        //go back to time - deltaT:
     52        time-=dt;
     53
     54        IssmDouble timeinputs = time;
     55
     56        /*Start loop: */
     57        count=1;
     58        for (t=time;t<time+dt;t=t+smb_dt){
     59
     60                for(int i=0;i<femmodel->elements->Size();i++){
     61                        Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
     62
     63                        timeclim=time;
     64                        if (isclimatology){
     65                                //If this is a climatology, we need to repeat the forcing after the final time
     66                                TransientInput2* Ta_input_tr  = element->inputs2->GetTransientInput(SmbTaEnum);    _assert_(Ta_input_tr);
     67
     68                                /*Get temperature climatology value*/
     69                                int offsetend = Ta_input_tr->GetTimeInputOffset(finaltime);
     70                                IssmDouble time0     = Ta_input_tr->GetTimeByOffset(-1);
     71                                IssmDouble timeend   = Ta_input_tr->GetTimeByOffset(offsetend);
     72                                if (time>time0 & timeend>time0){
     73                                        delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
     74                                        timeclim=time0+delta;
     75                                }
     76                        }
     77                        timeinputs = t-time+timeclim;
     78                        element->SmbGemb(timeinputs,count);
     79                }
     80                count=count+1;
    3781        }
    3882
    3983} /*}}}*/
    … …  
    10441088
    10451089                        // SWs and SWss coefficients need to be better constranted. Greuell
    10461090                        // and Konzelmann 1994 used SWs = 0.36 and SWss = 0.64 as this the
    1047                         // the // of SW radiation with wavelengths > and < 800 nm
     1091                        // the % of SW radiation with wavelengths > and < 800 nm
    10481092                        // respectively.  This, however, may not account for the fact that
    10491093                        // the albedo of wavelengths > 800 nm has a much lower albedo.
    10501094
    … …  
    21062150        if(V< 0.01-Dtol)V=0.01;
    21072151
    21082152        // calculate the Bulk Richardson Number (Ri)
    2109         Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2));
     2153        Ri = (2.0*9.81* (Vz - z0) * (Ta - Ts)) / ((Ta + Ts)* pow(V,2.0));
    21102154
    2111         // calculate Monin-Obukhov stability factors 'coef_M' and 'coef_H'
     2155        // calculate Monin-Obukhov stability factors 'coefM' and 'coefH'
    21122156
    21132157        // do not allow Ri to exceed 0.19
    21142158        Ri = min(Ri, 0.19);

  • ../trunk-jpl/src/c/classes/Elements/Element.h

     
    171171                void               SetIntInput(Inputs2* inputs2,int enum_in,int value);
    172172                void               SmbSemic();
    173173                int                Sid();
    174                 void               SmbGemb();
     174                void               SmbGemb(IssmDouble timeinputs, int count);
    175175                void               StrainRateESA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input);
    176176                void               StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input,Input2* vz_input);
    177177                void               StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input2* vx_input,Input2* vy_input);

  • ../trunk-jpl/src/c/classes/Elements/Element.cpp

     
    561561        this->parameters->FindParam(&isTemperatureScaled,SmbIstemperaturescaledEnum);
    562562        this->parameters->FindParam(&isPrecipScaled,SmbIsprecipscaledEnum);
    563563
    564                 /*Recover present day temperature and precipitation*/
     564        /*Recover present day temperature and precipitation*/
    565565        DatasetInput2 *dinput3 = NULL;
    566566        DatasetInput2 *dinput4 = NULL;
    567567        int            offset_t,offset_p,N;
    … …  
    35873587
    35883588}
    35893589/*}}}*/
    3590 void       Element::SmbGemb(){/*{{{*/
     3590void       Element::SmbGemb(IssmDouble timeinputs, int count){/*{{{*/
    35913591
    35923592        /*only compute SMB at the surface: */
    35933593        if (!IsOnSurface()) return;
    … …  
    36043604        IssmDouble Vmean=0.0;
    36053605        IssmDouble C=0.0;
    36063606        IssmDouble Tz,Vz=0.0;
    3607         IssmDouble time,dt,starttime,finaltime;
    3608         IssmDouble timeclim=0.0;
    3609         IssmDouble t,smb_dt;
    36103607        IssmDouble yts;
    36113608        IssmDouble Ta=0.0;
    36123609        IssmDouble V=0.0;
    … …  
    36173614        IssmDouble pAir=0.0;
    36183615        IssmDouble teValue=1.0;
    36193616        IssmDouble aValue=0.0;
     3617        IssmDouble dt,time,smb_dt;
    36203618        int        aIdx=0;
    36213619        int        denIdx=0;
    36223620        int        dsnowIdx=0;
    … …  
    36263624        IssmDouble shf=0.0;
    36273625        IssmDouble dayEC=0.0;
    36283626        IssmDouble initMass=0.0;
    3629         IssmDouble sumR=0.0;
    3630         IssmDouble sumM=0.0;
    3631         IssmDouble sumEC=0.0;
    3632         IssmDouble sumP=0.0;
    3633         IssmDouble sumW=0.0;
    3634         IssmDouble sumMassAdd=0.0;
    3635         IssmDouble sumdz_add=0.0;
     3627   IssmDouble sumR=0.0;
     3628   IssmDouble sumM=0.0;
     3629   IssmDouble sumEC=0.0;
     3630   IssmDouble sumP=0.0;
     3631   IssmDouble sumW=0.0;
     3632   IssmDouble sumMassAdd=0.0;
    36363633        IssmDouble fac=0.0;
    36373634        IssmDouble sumMass=0.0;
    36383635        IssmDouble dMass=0.0;
    … …  
    36413638        IssmDouble init_scaling=0.0;
    36423639        IssmDouble thermo_scaling=1.0;
    36433640        IssmDouble adThresh=1023.0;
    3644         int offsetend=-1;
    3645         IssmDouble time0, timeend, delta;
    36463641        /*}}}*/
    36473642        /*Output variables:{{{ */
    36483643        IssmDouble* dz=NULL;
    … …  
    36753670        IssmDouble* aini = NULL;
    36763671        IssmDouble* Tini = NULL;
    36773672        int         m=0;
    3678         int         count=0;
    36793673        /*}}}*/
    36803674
    36813675        /*Retrieve material properties and parameters:{{{ */
    … …  
    36863680        parameters->FindParam(&time,TimeEnum);                        /*transient core time at which we run the smb core*/
    36873681        parameters->FindParam(&dt,TimesteppingTimeStepEnum);          /*transient core time step*/
    36883682        parameters->FindParam(&yts,ConstantsYtsEnum);
    3689         parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
    3690         parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
    36913683        parameters->FindParam(&smb_dt,SmbDtEnum);                     /*time period for the smb solution,  usually smaller than the glaciological dt*/
    36923684        parameters->FindParam(&aIdx,SmbAIdxEnum);
    36933685        parameters->FindParam(&denIdx,SmbDenIdxEnum);
    … …  
    36973689        parameters->FindParam(&t0wet,SmbT0wetEnum);
    36983690        parameters->FindParam(&t0dry,SmbT0dryEnum);
    36993691        parameters->FindParam(&K,SmbKEnum);
    3700         parameters->FindParam(&isclimatology,SmbIsclimatologyEnum);
    37013692        parameters->FindParam(&isgraingrowth,SmbIsgraingrowthEnum);
    37023693        parameters->FindParam(&isalbedo,SmbIsalbedoEnum);
    37033694        parameters->FindParam(&isshortwave,SmbIsshortwaveEnum);
    … …  
    37223713        Input2 *C_input             = this->GetInput2(SmbCEnum);            _assert_(C_input);
    37233714        Input2 *Tz_input            = this->GetInput2(SmbTzEnum);           _assert_(Tz_input);
    37243715        Input2 *Vz_input            = this->GetInput2(SmbVzEnum);           _assert_(Vz_input);
    3725         Input2 *teValue_input       = this->GetInput2(SmbTeValueEnum);      _assert_(teValue_input);
    3726         Input2 *aValue_input        = this->GetInput2(SmbAValueEnum);       _assert_(aValue_input);
    37273716        Input2 *EC_input            = NULL;
    37283717
    37293718        /*Retrieve input values:*/
    … …  
    37413730        C_input->GetInputValue(&C,gauss);
    37423731        Tz_input->GetInputValue(&Tz,gauss);
    37433732        Vz_input->GetInputValue(&Vz,gauss);
    3744         teValue_input->GetInputValue(&teValue,gauss);
    3745         aValue_input->GetInputValue(&aValue,gauss);
    37463733        /*}}}*/
    37473734
    37483735        /*First, check that the initial structures have been setup in GEMB. If not, initialize profile variables: layer thickness dz, * density d, temperature T, etc. {{{*/
    … …  
    37623749                EC_input = this->GetInput2(SmbECiniEnum);  _assert_(EC_input);
    37633750                EC_input->GetInputAverage(&EC);
    37643751
    3765                 /*Retrive the correct value of m (without the zeroes at the end)*/
     3752                /*Retrieve the correct value of m (without the zeroes at the end)*/
    37663753                this->GetInput2Value(&m,SmbSizeiniEnum);
    37673754
    37683755                if(m==2){ //Snow properties are initialized with default values. Vertical grid has to be initialized too
    … …  
    38153802                this->inputs2->GetArray(SmbWEnum,this->lid,&W,&m);
    38163803                this->inputs2->GetArray(SmbAEnum,this->lid,&a,&m);
    38173804                this->inputs2->GetArray(SmbTEnum,this->lid,&T,&m);
    3818                 EC_input = this->GetInput2(SmbECEnum);  _assert_(EC_input);
     3805                EC_input = this->GetInput2(SmbECDtEnum);  _assert_(EC_input);
    38193806                EC_input->GetInputAverage(&EC);
    3820                 EC=EC*dt*rho_ice;
    38213807
    38223808                //fixed lower temperature bounday condition - T is fixed
    38233809                _assert_(m>0);
    … …  
    38293815
    38303816        // initialize cumulative variables
    38313817        sumR = 0; sumM = 0; sumEC = 0; sumP = 0; sumMassAdd = 0;
    3832         sumdz_add=0;
    38333818
    38343819        //before starting loop, realize that the transient core runs this smb_core at time = time +deltaT.
    3835         //go back to time - deltaT:
    3836         time-=dt;
     3820   //go back to time - deltaT:
     3821   time-=dt;
    38373822
    3838         timeclim=time;
    3839         if (isclimatology){
    3840                 //If this is a climatology, we need to repeat the forcing after the final time
    3841                 TransientInput2* Ta_input_tr  = this->inputs2->GetTransientInput(SmbTaEnum);    _assert_(Ta_input_tr);
    3842                 offsetend = Ta_input_tr->GetTimeInputOffset(finaltime);
    3843                 time0     = Ta_input_tr->GetTimeByOffset(-1);
    3844                 timeend   = Ta_input_tr->GetTimeByOffset(offsetend);
    3845                 if (time>time0 & timeend>time0){
    3846                         delta=(time-time0) - (timeend-time0)*(reCast<int,IssmDouble>((time-time0)/(timeend-time0)));
    3847                         timeclim=time0+delta;
    3848                 }
    3849         }
     3823        if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << timeinputs/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");
    38503824
    3851         /*Start loop: */
    3852         count=1;
    3853         for (t=time;t<time+dt;t=t+smb_dt){
     3825        /*Get daily accumulated inputs {{{*/
     3826        if (count>1){
     3827                Input2 *sumEC_input         = this->GetInput2(SmbECEnum);  _assert_(sumEC_input);
     3828                Input2 *sumM_input          = this->GetInput2(SmbMeltEnum);  _assert_(sumM_input);
     3829                Input2 *sumR_input          = this->GetInput2(SmbRunoffEnum);  _assert_(sumR_input);
     3830                Input2 *sumP_input          = this->GetInput2(SmbPrecipitationEnum);  _assert_(sumP_input);
     3831                Input2 *ULW_input           = this->GetInput2(SmbMeanULWEnum);  _assert_(ULW_input);
     3832                Input2 *LW_input            = this->GetInput2(SmbNetLWEnum);  _assert_(LW_input);
     3833                Input2 *SW_input            = this->GetInput2(SmbNetSWEnum);  _assert_(SW_input);
     3834                Input2 *LHF_input           = this->GetInput2(SmbMeanLHFEnum);  _assert_(LHF_input);
     3835                Input2 *SHF_input           = this->GetInput2(SmbMeanSHFEnum);  _assert_(SHF_input);
     3836                Input2 *DzAdd_input         = this->GetInput2(SmbDzAddEnum);  _assert_(DzAdd_input);
     3837                Input2 *MassAdd_input       = this->GetInput2(SmbMAddEnum);  _assert_(MassAdd_input);
     3838                Input2 *InitMass_input      = this->GetInput2(SmbMInitnum);  _assert_(InitMass_input);
    38543839
    3855                 if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0)_printf0_("Time: t=" << setprecision(8) << t/365.0/24.0/3600.0 << " yr/" << (time+dt)/365.0/24.0/3600.0 << " yr" << setprecision(3) << " Step: " << count << "\n");
     3840                ULW_input->GetInputAverage(&meanULW);
     3841                LW_input->GetInputAverage(&netLW);
     3842                SW_input->GetInputAverage(&netSW);
     3843                LHF_input->GetInputAverage(&meanLHF);
     3844                SHF_input->GetInputAverage(&meanSHF);
     3845                DzAdd_input->GetInputAverage(&dz_add);
     3846                MassAdd_input->GetInputAverage(&sumMassAdd);
     3847                sumMassAdd=sumMassAdd*dt;
     3848                InitMass_input->GetInputAverage(&initMass);
     3849                sumEC_input->GetInputAverage(&sumEC);
     3850                sumEC=sumEC*dt*rho_ice;
     3851                sumM_input->GetInputAverage(&sumM);
     3852                sumM=sumM*dt*rho_ice;
     3853                sumR_input->GetInputAverage(&sumR);
     3854                sumR=sumR*dt*rho_ice;
     3855                sumP_input->GetInputAverage(&sumP);
     3856                sumP=sumP*dt*rho_ice;
     3857        }
     3858        /*}}}*/
    38563859
    3857                 Input2 *Ta_input  = this->GetInput2(SmbTaEnum,t-time+timeclim);    _assert_(Ta_input);
    3858                 Input2 *V_input   = this->GetInput2(SmbVEnum,t-time+timeclim);     _assert_(V_input);
    3859                 Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,t-time+timeclim); _assert_(Dlwr_input);
    3860                 Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,t-time+timeclim); _assert_(Dswr_input);
    3861                 Input2 *P_input   = this->GetInput2(SmbPEnum,t-time+timeclim);     _assert_(P_input);
    3862                 Input2 *eAir_input= this->GetInput2(SmbEAirEnum,t-time+timeclim);  _assert_(eAir_input);
    3863                 Input2 *pAir_input= this->GetInput2(SmbPAirEnum,t-time+timeclim);  _assert_(pAir_input);
     3860        // Get time forcing inputs
     3861        Input2 *Ta_input  = this->GetInput2(SmbTaEnum,timeinputs);    _assert_(Ta_input);
     3862        Input2 *V_input   = this->GetInput2(SmbVEnum,timeinputs);     _assert_(V_input);
     3863        Input2 *Dlwr_input= this->GetInput2(SmbDlwrfEnum,timeinputs); _assert_(Dlwr_input);
     3864        Input2 *Dswr_input= this->GetInput2(SmbDswrfEnum,timeinputs); _assert_(Dswr_input);
     3865        Input2 *P_input   = this->GetInput2(SmbPEnum,timeinputs);     _assert_(P_input);
     3866        Input2 *eAir_input= this->GetInput2(SmbEAirEnum,timeinputs);  _assert_(eAir_input);
     3867        Input2 *pAir_input= this->GetInput2(SmbPAirEnum,timeinputs);  _assert_(pAir_input);
     3868        Input2 *teValue_input= this->GetInput2(SmbTeValueEnum,timeinputs); _assert_(teValue_input);
     3869        Input2 *aValue_input= this->GetInput2(SmbAValueEnum,timeinputs); _assert_(aValue_input);
    38643870
    3865                 /*extract daily data:{{{*/
    3866                 Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
    3867                 V_input->GetInputValue(&V,gauss);  //wind speed [m s-1]
    3868                 Dlwr_input->GetInputValue(&dlw,gauss);   //downward longwave radiation flux [W m-2]
    3869                 Dswr_input->GetInputValue(&dsw,gauss);   //downward shortwave radiation flux [W m-2]
    3870                 P_input->GetInputValue(&P,gauss);        //precipitation [kg m-2]
    3871                 eAir_input->GetInputValue(&eAir,gauss);  //screen level vapor pressure [Pa]
    3872                 pAir_input->GetInputValue(&pAir,gauss);  // screen level air pressure [Pa]
    3873                 teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
    3874                 aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
    3875                 //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
    3876                 /*}}}*/
     3871        /*extract daily data:{{{*/
     3872        Ta_input->GetInputValue(&Ta,gauss);//screen level air temperature [K]
     3873        V_input->GetInputValue(&V,gauss);  //wind speed [m s-1]
     3874        Dlwr_input->GetInputValue(&dlw,gauss);   //downward longwave radiation flux [W m-2]
     3875        Dswr_input->GetInputValue(&dsw,gauss);   //downward shortwave radiation flux [W m-2]
     3876        P_input->GetInputValue(&P,gauss);        //precipitation [kg m-2]
     3877        eAir_input->GetInputValue(&eAir,gauss);  //screen level vapor pressure [Pa]
     3878        pAir_input->GetInputValue(&pAir,gauss);  // screen level air pressure [Pa]
     3879        teValue_input->GetInputValue(&teValue,gauss);  // Emissivity [0-1]
     3880        aValue_input->GetInputValue(&aValue,gauss);  // Albedo [0 1]
     3881        //_printf_("Time: " << t << " Ta: " << Ta << " V: " << V << " dlw: " << dlw << " dsw: " << dsw << " P: " << P << " eAir: " << eAir << " pAir: " << pAir << "\n");
     3882        /*}}}*/
    38773883
    3878                 /*Snow grain metamorphism:*/
    3879                 if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
     3884        /*Snow grain metamorphism:*/
     3885        if(isgraingrowth)grainGrowth(&re, &gdn, &gsp, T, dz, d, W, smb_dt, m, aIdx,this->Sid());
    38803886
    3881                 /*Snow, firn and ice albedo:*/
    3882                 if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
     3887        /*Snow, firn and ice albedo:*/
     3888        if(isalbedo)albedo(&a,aIdx,re,d,cldFrac,aIce,aSnow,aValue,adThresh,T,W,P,EC,t0wet,t0dry,K,smb_dt,rho_ice,m,this->Sid());
    38833889
    3884                 /*Distribution of absorbed short wave radation with depth:*/
    3885                 if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
     3890        /*Distribution of absorbed short wave radation with depth:*/
     3891        if(isshortwave)shortwave(&swf, swIdx, aIdx, dsw, a[0], d, dz, re,rho_ice,m,this->Sid());
    38863892
    3887                 /*Calculate net shortwave [W m-2]*/
    3888                 netSW = netSW + cellsum(swf,m)*smb_dt/dt;
     3893        /*Calculate net shortwave [W m-2]*/
     3894        netSW = netSW + cellsum(swf,m)*smb_dt/dt;
    38893895
    3890                 /*Thermal profile computation:*/
    3891                 if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
     3896        /*Thermal profile computation:*/
     3897        if(isthermal)thermo(&EC, &T, &ulw, dz, d, swf, dlw, Ta, V, eAir, pAir, teValue, W[0], smb_dt, m, Vz, Tz, thermo_scaling,rho_ice,this->Sid());
    38923898
    3893                 /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
    3894                 * need to fix this in case all or more of cell evaporates */
    3895                 dz[0] = dz[0] + EC / d[0];
     3899        /*Change in thickness of top cell due to evaporation/condensation  assuming same density as top cell.
     3900        * need to fix this in case all or more of cell evaporates */
     3901        dz[0] = dz[0] + EC / d[0];
    38963902
    3897                 /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
    3898                 if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
     3903        /*Add snow/rain to top grid cell adjusting cell depth, temperature and density*/
     3904        if(isaccumulation)accumulation(&T, &dz, &d, &W, &a, &re, &gdn, &gsp, &m, aIdx, dsnowIdx, Tmean, Ta, P, dzMin, aSnow, C, V, Vmean, rho_ice,this->Sid());
    38993905
    3900                 /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
    3901                 * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
    3902                 if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
     3906        /*Calculate water production, M [kg m-2] resulting from snow/ice temperature exceeding 273.15 deg K
     3907        * (> 0 deg C), runoff R [kg m-2] and resulting changes in density and determine wet compaction [m]*/
     3908        if(ismelt)melt(&M, &R, &mAdd, &dz_add, &T, &d, &dz, &W, &a, &re, &gdn, &gsp, &m, dzMin, zMax, zMin, zTop,rho_ice,this->Sid());
    39033909
    3904                 /*Allow non-melt densification and determine compaction [m]*/
    3905                 if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
     3910        /*Allow non-melt densification and determine compaction [m]*/
     3911        if(isdensification)densification(&d,&dz, T, re, denIdx, C, smb_dt, Tmean,rho_ice,m,this->Sid());
    39063912
    3907                 /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
    3908                 * sub-time step in thermo equations*/
    3909                 //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
     3913        /*Calculate upward longwave radiation flux [W m-2] not used in energy balance. Calculated for every
     3914        * sub-time step in thermo equations*/
     3915        //ulw = 5.67E-8 * pow(T[0],4.0) * teValue; // + deltatest here
    39103916
    3911                 /*Calculate net longwave [W m-2]*/
    3912                 meanULW = meanULW + ulw*smb_dt/dt;
    3913                 netLW = netLW + (dlw - ulw)*smb_dt/dt;
     3917        /*Calculate net longwave [W m-2]*/
     3918        meanULW = meanULW + ulw*smb_dt/dt;
     3919        netLW = netLW + (dlw - ulw)*smb_dt/dt;
    39143920
    3915                 /*Calculate turbulent heat fluxes [W m-2]*/
    3916                 if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
     3921        /*Calculate turbulent heat fluxes [W m-2]*/
     3922        if(isturbulentflux)turbulentFlux(&shf, &lhf, &dayEC, Ta, T[0], V, eAir, pAir, d[0], W[0], Vz, Tz,rho_ice,this->Sid());
    39173923
    3918                 /*Verbose some results in debug mode: {{{*/
    3919                 if(VerboseSmb() && 0){
    3920                         _printf_("smb log: count[" << count << "] m[" << m << "] "
    3921                                                 << setprecision(16)   << "T[" << cellsum(T,m)  << "] "
    3922                                                 << "d[" << cellsum(d,m)  << "] "
    3923                                                 << "dz[" << cellsum(dz,m)  << "] "
    3924                                                 << "a[" << cellsum(a,m)  << "] "
    3925                                                 << "W[" << cellsum(W,m)  << "] "
    3926                                                 << "re[" << cellsum(re,m)  << "] "
    3927                                                 << "gdn[" << cellsum(gdn,m)  << "] "
    3928                                                 << "gsp[" << cellsum(gsp,m)  << "] "
    3929                                                 << "swf[" << netSW << "] "
    3930                                                 << "lwf[" << netLW << "] "
    3931                                                 << "a[" << a << "] "
    3932                                                 << "te[" << teValue << "] "
    3933                                                 << "\n");
    3934                 } /*}}}*/
     3924        /*Verbose some results in debug mode: {{{*/
     3925        if(VerboseSmb() && 0){
     3926                _printf_("smb log: count[" << count << "] m[" << m << "] "
     3927                                        << setprecision(16)   << "T[" << cellsum(T,m)  << "] "
     3928                                        << "d[" << cellsum(d,m)  << "] "
     3929                                        << "dz[" << cellsum(dz,m)  << "] "
     3930                                        << "a[" << cellsum(a,m)  << "] "
     3931                                        << "W[" << cellsum(W,m)  << "] "
     3932                                        << "re[" << cellsum(re,m)  << "] "
     3933                                        << "gdn[" << cellsum(gdn,m)  << "] "
     3934                                        << "gsp[" << cellsum(gsp,m)  << "] "
     3935                                        << "swf[" << netSW << "] "
     3936                                        << "lwf[" << netLW << "] "
     3937                                        << "a[" << a << "] "
     3938                                        << "te[" << teValue << "] "
     3939                                        << "\n");
     3940        } /*}}}*/
    39353941
    3936                 meanLHF = meanLHF + lhf*smb_dt/dt;
    3937                 meanSHF = meanSHF + shf*smb_dt/dt;
     3942        meanLHF = meanLHF + lhf*smb_dt/dt;
     3943        meanSHF = meanSHF + shf*smb_dt/dt;
    39383944
    3939                 /*Sum component mass changes [kg m-2]*/
    3940                 sumMassAdd = mAdd + sumMassAdd;
    3941                 sumM = M + sumM;
    3942                 sumR = R + sumR;
    3943                 sumW = cellsum(W,m);
    3944                 sumP = P +  sumP;
    3945                 sumEC = sumEC + EC;  // evap (-)/cond(+)
     3945        /*Sum component mass changes [kg m-2]*/
     3946        sumMassAdd = mAdd + sumMassAdd;
     3947        sumM = M + sumM;
     3948        sumR = R + sumR;
     3949        sumW = cellsum(W,m);
     3950        sumP = P +  sumP;
     3951        sumEC = sumEC + EC;  // evap (-)/cond(+)
    39463952
    3947                 /*Calculate total system mass:*/
    3948                 sumMass=0;
    3949                 fac=0;
    3950                 for(int i=0;i<m;i++){
    3951                         sumMass += dz[i]*d[i];
    3952                         fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
    3953                 }
     3953        /*Calculate total system mass:*/
     3954        sumMass=0;
     3955        fac=0;
     3956        for(int i=0;i<m;i++){
     3957                sumMass += dz[i]*d[i];
     3958                fac += dz[i]*(rho_ice - fmin(d[i],rho_ice));
     3959        }
    39543960
    3955                 #if defined(_HAVE_AD_)
    3956                 /*we want to avoid the round operation at all cost. Not differentiable.*/
    3957                 _error_("not implemented yet");
    3958                 #else
    3959                 dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
    3960                 dMass = round(dMass * 100.0)/100.0;
     3961        #if defined(_HAVE_AD_)
     3962        /*we want to avoid the round operation at all cost. Not differentiable.*/
     3963        _error_("not implemented yet");
     3964        #else
     3965        dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd;
     3966        dMass = round(dMass * 100.0)/100.0;
    39613967
    3962                 /*Check mass conservation:*/
    3963                 if (dMass != 0.0) _printf_("total system mass not conserved in MB function \n");
    3964                 #endif
     3968        /*Check mass conservation:*/
     3969        if (dMass != 0.0){
     3970                _printf_("total system mass not conserved in MB function \n");
     3971        }
     3972        #endif
    39653973
    3966                 /*Check bottom grid cell T is unchanged:*/
    3967                 if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
    3968                         if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
    3969                 }
     3974        /*Check bottom grid cell T is unchanged:*/
     3975        if(VerboseSmb() && this->Sid()==0 && IssmComm::GetRank()==0){
     3976                if (T[m-1]!=T_bottom) _printf_("T(end)~=T_bottom" << "\n");
     3977        }
    39703978
    3971                 /*Free ressources: */
    3972                 xDelete<IssmDouble>(swf);
    3973 
    3974                 /*increase counter:*/
    3975                 count++;
    3976         } //for (t=time;t<time+dt;t=t+smb_dt)
    3977 
    39783979        /*Save generated inputs: */
    39793980        this->inputs2->SetArrayInput(SmbDzEnum,this->lid,dz,m);
    39803981        this->inputs2->SetArrayInput(SmbDEnum,this->lid,d,m);
    … …  
    39943995        this->SetElementInput(SmbNetSWEnum,netSW);
    39953996        this->SetElementInput(SmbMeanLHFEnum,meanLHF);
    39963997        this->SetElementInput(SmbMeanSHFEnum,meanSHF);
    3997         this->SetElementInput(SmbDzAddEnum,sumdz_add);
     3998        this->SetElementInput(SmbDzAddEnum,dz_add);
     3999        this->SetElementInput(SmbMInitnum,initMass);
    39984000        this->SetElementInput(SmbMAddEnum,sumMassAdd/dt);
     4001        this->SetElementInput(SmbWAddEnum,sumW/dt);
    39994002        this->SetElementInput(SmbFACEnum,fac/1000.); // output in meters
     4003        this->SetElementInput(SmbECDtEnum,EC);
    40004004
    40014005        /*Free allocations:{{{*/
    40024006        if(dz) xDelete<IssmDouble>(dz);
    … …  
    40154019        if(Wini) xDelete<IssmDouble>(Wini);
    40164020        if(aini) xDelete<IssmDouble>(aini);
    40174021        if(Tini) xDelete<IssmDouble>(Tini);
     4022        if(swf) xDelete<IssmDouble>(swf);
    40184023
    40194024        delete gauss;
    40204025        /*}}}*/
Note: See TracBrowser for help on using the repository browser.