Index: ../trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h =================================================================== --- ../trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 23241) +++ ../trunk-jpl/src/c/classes/ExternalResults/GenericExternalResult.h (revision 23242) @@ -621,7 +621,7 @@ } /*}}}*/ /*Specific instantiations for IssmDouble*: */ -#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization. +#if defined(_HAVE_AD_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization. template <> inline GenericExternalResult::~GenericExternalResult(){ /*{{{*/ xDelete(result_name); xDelete(value); @@ -708,7 +708,7 @@ template <> inline Object* GenericExternalResult*>::copy(void){ /*{{{*/ return new GenericExternalResult*>(this->id,StringToEnumx(this->result_name),this->value,this->step,this->time); } /*}}}*/ -#if defined(_HAVE_ADOLC_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization. +#if defined(_HAVE_AD_) && !defined(_WRAPPERS_) //We hook off this specific specialization when not running ADOLC, otherwise we get a redeclaration with the next specialization. template <> inline void GenericExternalResult*>::WriteData(FILE* fid,bool io_gather){ /*{{{*/ char *name = NULL; Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23241) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 23242) @@ -880,7 +880,7 @@ results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentMemEnum, solution_memory)); results->AddObject(new GenericExternalResult(results->Size()+1, ProfilingCurrentFlopsEnum, solution_flops)); - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ solution_time = profiler->TotalTime(ADCORE); solution_flops = profiler->TotalFlops(ADCORE); solution_memory = profiler->Memory(ADCORE); @@ -1844,7 +1844,7 @@ parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); if(isautodiff){ - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum); if(num_dependents){ Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23241) +++ ../trunk-jpl/src/c/classes/Elements/Element.cpp (revision 23242) @@ -3132,7 +3132,10 @@ fac += dz[i]*(rho_ice - fmin(d[i],rho_ice)); } - #ifndef _HAVE_ADOLC_ //we want to avoid the round operation at all cost. Not differentiable. + #if defined(_HAVE_AD_) + /*we want to avoid the round operation at all cost. Not differentiable.*/ + _error_("not implemented yet"); + #else dMass = sumMass + sumR + sumW - sumP - sumEC - initMass - sumMassAdd; dMass = round(dMass * 100.0)/100.0; Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23241) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 23242) @@ -52,6 +52,7 @@ AutodiffNumDependentsEnum, AutodiffNumIndependentsEnum, AutodiffObufsizeEnum, + AutodiffTapeAllocEnum, AutodiffTbufsizeEnum, AutodiffXpEnum, BalancethicknessStabilizationEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23241) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 23242) @@ -60,6 +60,7 @@ case AutodiffNumDependentsEnum : return "AutodiffNumDependents"; case AutodiffNumIndependentsEnum : return "AutodiffNumIndependents"; case AutodiffObufsizeEnum : return "AutodiffObufsize"; + case AutodiffTapeAllocEnum : return "AutodiffTapeAlloc"; case AutodiffTbufsizeEnum : return "AutodiffTbufsize"; case AutodiffXpEnum : return "AutodiffXp"; case BalancethicknessStabilizationEnum : return "BalancethicknessStabilization"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23241) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 23242) @@ -60,6 +60,7 @@ else if (strcmp(name,"AutodiffNumDependents")==0) return AutodiffNumDependentsEnum; else if (strcmp(name,"AutodiffNumIndependents")==0) return AutodiffNumIndependentsEnum; else if (strcmp(name,"AutodiffObufsize")==0) return AutodiffObufsizeEnum; + else if (strcmp(name,"AutodiffTapeAlloc")==0) return AutodiffTapeAllocEnum; else if (strcmp(name,"AutodiffTbufsize")==0) return AutodiffTbufsizeEnum; else if (strcmp(name,"AutodiffXp")==0) return AutodiffXpEnum; else if (strcmp(name,"BalancethicknessStabilization")==0) return BalancethicknessStabilizationEnum; @@ -135,11 +136,11 @@ else if (strcmp(name,"FrictionGamma")==0) return FrictionGammaEnum; else if (strcmp(name,"FrictionLaw")==0) return FrictionLawEnum; else if (strcmp(name,"FrictionPseudoplasticityExponent")==0) return FrictionPseudoplasticityExponentEnum; - else if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum; else stage=2; } if(stage==2){ - if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum; + if (strcmp(name,"FrictionThresholdSpeed")==0) return FrictionThresholdSpeedEnum; + else if (strcmp(name,"FrictionDelta")==0) return FrictionDeltaEnum; else if (strcmp(name,"FrictionVoidRatio")==0) return FrictionVoidRatioEnum; else if (strcmp(name,"GiaCrossSectionShape")==0) return GiaCrossSectionShapeEnum; else if (strcmp(name,"GroundinglineMigration")==0) return GroundinglineMigrationEnum; @@ -258,11 +259,11 @@ else if (strcmp(name,"SealevelriseEquatorialMoi")==0) return SealevelriseEquatorialMoiEnum; else if (strcmp(name,"SealevelriseFluidLove")==0) return SealevelriseFluidLoveEnum; else if (strcmp(name,"SealevelriseGElastic")==0) return SealevelriseGElasticEnum; - else if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum; else stage=3; } if(stage==3){ - if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum; + if (strcmp(name,"SealevelriseGeodetic")==0) return SealevelriseGeodeticEnum; + else if (strcmp(name,"SealevelriseGeodeticRunFrequency")==0) return SealevelriseGeodeticRunFrequencyEnum; else if (strcmp(name,"SealevelriseHElastic")==0) return SealevelriseHElasticEnum; else if (strcmp(name,"SealevelriseHoriz")==0) return SealevelriseHorizEnum; else if (strcmp(name,"SealevelriseLoopIncrement")==0) return SealevelriseLoopIncrementEnum; @@ -381,11 +382,11 @@ else if (strcmp(name,"TransientNumRequestedOutputs")==0) return TransientNumRequestedOutputsEnum; else if (strcmp(name,"TransientRequestedOutputs")==0) return TransientRequestedOutputsEnum; else if (strcmp(name,"Velocity")==0) return VelocityEnum; - else if (strcmp(name,"WorldComm")==0) return WorldCommEnum; else stage=4; } if(stage==4){ - if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; + if (strcmp(name,"WorldComm")==0) return WorldCommEnum; + else if (strcmp(name,"ParametersEND")==0) return ParametersENDEnum; else if (strcmp(name,"InputsSTART")==0) return InputsSTARTEnum; else if (strcmp(name,"Adjoint")==0) return AdjointEnum; else if (strcmp(name,"Adjointx")==0) return AdjointxEnum; @@ -504,11 +505,11 @@ else if (strcmp(name,"InversionVelObs")==0) return InversionVelObsEnum; else if (strcmp(name,"InversionVxObs")==0) return InversionVxObsEnum; else if (strcmp(name,"InversionVyObs")==0) return InversionVyObsEnum; - else if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; else stage=5; } if(stage==5){ - if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; + if (strcmp(name,"LevelsetfunctionSlopeX")==0) return LevelsetfunctionSlopeXEnum; + else if (strcmp(name,"LevelsetfunctionSlopeY")==0) return LevelsetfunctionSlopeYEnum; else if (strcmp(name,"LoadingforceX")==0) return LoadingforceXEnum; else if (strcmp(name,"LoadingforceY")==0) return LoadingforceYEnum; else if (strcmp(name,"LoadingforceZ")==0) return LoadingforceZEnum; @@ -627,11 +628,11 @@ else if (strcmp(name,"StrainRateperpendicular")==0) return StrainRateperpendicularEnum; else if (strcmp(name,"StrainRatexx")==0) return StrainRatexxEnum; else if (strcmp(name,"StrainRatexy")==0) return StrainRatexyEnum; - else if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; else stage=6; } if(stage==6){ - if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; + if (strcmp(name,"StrainRatexz")==0) return StrainRatexzEnum; + else if (strcmp(name,"StrainRateyy")==0) return StrainRateyyEnum; else if (strcmp(name,"StrainRateyz")==0) return StrainRateyzEnum; else if (strcmp(name,"StrainRatezz")==0) return StrainRatezzEnum; else if (strcmp(name,"StressMaxPrincipal")==0) return StressMaxPrincipalEnum; @@ -750,11 +751,11 @@ else if (strcmp(name,"DeviatoricStressErrorEstimator")==0) return DeviatoricStressErrorEstimatorEnum; else if (strcmp(name,"Divergence")==0) return DivergenceEnum; else if (strcmp(name,"Domain3Dsurface")==0) return Domain3DsurfaceEnum; - else if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; else stage=7; } if(stage==7){ - if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; + if (strcmp(name,"DoubleArrayInput")==0) return DoubleArrayInputEnum; + else if (strcmp(name,"DoubleExternalResult")==0) return DoubleExternalResultEnum; else if (strcmp(name,"DoubleInput")==0) return DoubleInputEnum; else if (strcmp(name,"DoubleMatArrayParam")==0) return DoubleMatArrayParamEnum; else if (strcmp(name,"DoubleMatExternalResult")==0) return DoubleMatExternalResultEnum; @@ -873,11 +874,11 @@ else if (strcmp(name,"LoveHi")==0) return LoveHiEnum; else if (strcmp(name,"LoveHr")==0) return LoveHrEnum; else if (strcmp(name,"LoveKernelsImag")==0) return LoveKernelsImagEnum; - else if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"LoveKi")==0) return LoveKiEnum; + if (strcmp(name,"LoveKernelsReal")==0) return LoveKernelsRealEnum; + else if (strcmp(name,"LoveKi")==0) return LoveKiEnum; else if (strcmp(name,"LoveKr")==0) return LoveKrEnum; else if (strcmp(name,"LoveLi")==0) return LoveLiEnum; else if (strcmp(name,"LoveLr")==0) return LoveLrEnum; @@ -996,11 +997,11 @@ else if (strcmp(name,"Outputdefinition45")==0) return Outputdefinition45Enum; else if (strcmp(name,"Outputdefinition46")==0) return Outputdefinition46Enum; else if (strcmp(name,"Outputdefinition47")==0) return Outputdefinition47Enum; - else if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; + if (strcmp(name,"Outputdefinition48")==0) return Outputdefinition48Enum; + else if (strcmp(name,"Outputdefinition49")==0) return Outputdefinition49Enum; else if (strcmp(name,"Outputdefinition4")==0) return Outputdefinition4Enum; else if (strcmp(name,"Outputdefinition50")==0) return Outputdefinition50Enum; else if (strcmp(name,"Outputdefinition51")==0) return Outputdefinition51Enum; @@ -1119,11 +1120,11 @@ else if (strcmp(name,"SMBgradientsela")==0) return SMBgradientselaEnum; else if (strcmp(name,"SMBgradients")==0) return SMBgradientsEnum; else if (strcmp(name,"SMBhenning")==0) return SMBhenningEnum; - else if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; + if (strcmp(name,"SmbMAdd")==0) return SmbMAddEnum; + else if (strcmp(name,"SMBmeltcomponents")==0) return SMBmeltcomponentsEnum; else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SmbRlaps")==0) return SmbRlapsEnum; else if (strcmp(name,"SmbRlapslgm")==0) return SmbRlapslgmEnum; Index: ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23241) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 23242) @@ -373,7 +373,7 @@ /*Now, deal with toolkits options, which need to be put into the parameters dataset: */ ParseToolkitsOptionsx(parameters,toolkitsoptionsfid); - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ if(VerboseMProcessor()) _printf0_(" starting autodiff parameters \n"); CreateParametersAutodiff(parameters,iomodel); if(VerboseMProcessor()) _printf0_(" ending autodiff parameters \n"); Index: ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp =================================================================== --- ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23241) +++ ../trunk-jpl/src/c/modules/ModelProcessorx/Autodiff/CreateParametersAutodiff.cpp (revision 23242) @@ -72,8 +72,36 @@ /*Free ressources: */ xDelete(options); /*}}}*/ + #elif _HAVE_CODIPACK_ + //fprintf(stderr, "*** Codipack CreateParametersAutodiff()\n"); + /*initialize a placeholder to store solver pointers: {{{*/ + /*Solver pointers depend on what type of solver we are implementing: */ + options=OptionsFromAnalysis(parameters,DefaultAnalysisEnum); //options database is not filled in yet, use default. + ToolkitOptions::Init(options); + switch(IssmSolverTypeFromToolkitOptions()){ + case MumpsEnum:{ + #ifndef _HAVE_MUMPS_ + _error_("CoDiPack: requesting mumps solver without MUMPS being compiled in!"); + #endif + break; + } + case GslEnum: { + #ifndef _HAVE_GSL_ + _error_("CoDiPack: requesting GSL solver without GSL being compiled in!"); + #endif + break; + } + default: + _error_("solver type not supported yet!"); + } + /*Free ressources: */ + xDelete(options); + #endif + #if defined(_HAVE_AD_) + if(isautodiff){ +#if _HAVE_ADOLC_ /*Copy some parameters from IoModel to parameters dataset: {{{*/ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.obufsize",AutodiffObufsizeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.cbufsize",AutodiffCbufsizeEnum)); @@ -82,11 +110,16 @@ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerRatio",AutodiffGcTriggerRatioEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.gcTriggerMaxSize",AutodiffGcTriggerMaxSizeEnum)); /*}}}*/ +#endif /*retrieve driver: {{{*/ iomodel->FindConstant(&autodiff_driver,"md.autodiff.driver"); parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.driver",AutodiffDriverEnum)); if(strcmp(autodiff_driver,"fos_forward")==0){ +#if _HAVE_CODIPACK_ + // FIXME codi support Foward Mode (scalar) + _error_("Foward Mode (scalar) not supported yet!"); +#endif parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_forward_index",AutodiffFosForwardIndexEnum)); } else if(strcmp(autodiff_driver,"fos_reverse")==0){ @@ -93,6 +126,10 @@ parameters->AddObject(iomodel->CopyConstantObject("md.autodiff.fos_reverse_index",AutodiffFosReverseIndexEnum)); } else if(strcmp(autodiff_driver,"fov_forward")==0){ +#if _HAVE_CODIPACK_ + // FIXME codi support Foward Mode (vector) + _error_("Foward Mode (vector) not supported yet!"); +#endif /*Retrieve list of indices: */ iomodel->FetchData(&indices,&num_indices,&dummy,"md.autodiff.fov_forward_indices"); parameters->AddObject(new IntMatParam(AutodiffFovForwardIndicesEnum,indices,num_indices,1)); Index: ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp =================================================================== --- ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23241) +++ ../trunk-jpl/src/c/modules/SurfaceMassBalancex/Gembx.cpp (revision 23242) @@ -61,7 +61,7 @@ //check to see if the top grid cell structure length (dzTop) goes evenly //into specified top structure depth (zTop). Also make sure top grid cell //structure length (dzTop) is greater than 5 cm - #ifndef _HAVE_ADOLC_ //avoid the round operation check! + #ifndef _HAVE_AD_ //avoid the round operation check! if (dgpTop != round(dgpTop)){ _error_("top grid cell structure length does not go evenly into specified top structure depth, adjust dzTop or zTop"); } @@ -424,11 +424,11 @@ // spectral range: // 0.3 - 0.8um - IssmDouble a0 = fmin(0.98, 1 - 1.58 *pow(gsz,0.5)); + IssmDouble a0 = min(0.98, 1 - 1.58 *pow(gsz,0.5)); // 0.8 - 1.5um - IssmDouble a1 = fmax(0, 0.95 - 15.4 *pow(gsz,0.5)); + IssmDouble a1 = max(0, 0.95 - 15.4 *pow(gsz,0.5)); // 1.5 - 2.8um - IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5)); + IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz - 32.31*pow(gsz,0.5)); // broadband surface albedo a[0] = sF[0]*a0 + sF[1]*a1 + sF[2]*a2; @@ -666,7 +666,7 @@ // determine minimum acceptable delta t (diffusion number > 1/2) [s] // NS: 2.16.18 divided dt by scaling factor, default set to 1/11 for stability dt=1e12; - for(int i=0;i 1/2 // must go evenly into one hour or the data frequency if it is smaller @@ -751,7 +751,7 @@ // calculated. The estimated enegy balance & melt are significanly // less when Ts is taken as the mean of the x top grid cells. Ts = (T[0] + T[1])/2.0; - Ts = fmin(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC) + Ts = min(CtoK,Ts); // don't allow Ts to exceed 273.15 K (0 degC) //TURBULENT HEAT FLUX @@ -766,7 +766,7 @@ // calculate Monin-Obukhov stability factors 'coefM' and 'coefH' // do not allow Ri to exceed 0.19 - Ri = fmin(Ri, 0.19); + Ri = min(Ri, 0.19); // calculate momentum 'coefM' stability factor if (Ri > 0.0+Ttol){ @@ -946,11 +946,11 @@ // spectral albedos: // 0.3 - 0.8um - IssmDouble a0 = fmin(0.98, 1.0 - 1.58 *pow(gsz[0],0.5)); + IssmDouble a0 = min(0.98, 1.0 - 1.58 *pow(gsz[0],0.5)); // 0.8 - 1.5um - IssmDouble a1 = fmax(0.0, 0.95 - 15.4 *pow(gsz[0],0.5)); + IssmDouble a1 = max(0.0, 0.95 - 15.4 *pow(gsz[0],0.5)); // 1.5 - 2.8um - IssmDouble a2 = fmax(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5)); + IssmDouble a2 = max(0.127, 0.88 + 346.3*gsz[0] - 32.31*pow(gsz[0],0.5)); // separate net shortwave radiative flux into spectral ranges IssmDouble swfS[3]; @@ -1203,7 +1203,7 @@ mass_diff = mass - massinit - P; - #ifndef _HAVE_ADOLC_ //avoid round operation. only check in forward mode. + #ifndef _HAVE_AD_ //avoid round operation. only check in forward mode. mass_diff = round(mass_diff * 100.0)/100.0; if (mass_diff > 0) _error_("mass not conserved in accumulation function"); #endif @@ -1322,11 +1322,11 @@ // calculate temperature excess above 0 deg C exsT=xNewZeroInit(n); - for(int i=0;i0.0+Wtol){ if(VerboseSmb() && sid==0 && IssmComm::GetRank()==0)_printf0_(" pore water refreeze\n"); // calculate maximum freeze amount, maxF [kg] - for(int i=0;i(n); for(int i=0;i(n); for(int i=0;i(n); for(int i=0;i 0.0 + Ttol ){ // _printf_("T Surplus"); @@ -1379,10 +1379,10 @@ // use surplus energy to increase the temperature of lower cell T[i+1] = surpE[i]/m[i+1]/CI + T[i+1]; - exsT[i+1] = fmax(0.0, T[i+1] - CtoK) + exsT[i+1]; - T[i+1] = fmin(CtoK, T[i+1]); + exsT[i+1] = max(0.0, T[i+1] - CtoK) + exsT[i+1]; + T[i+1] = min(CtoK, T[i+1]); - surpT[i+1] = fmax(0.0, exsT[i+1] - LF/CI); + surpT[i+1] = max(0.0, exsT[i+1] - LF/CI); surpE[i+1] = surpT[i+1] * CI * m[i+1]; // adjust current cell properties (again 159.1342 is the max T) @@ -1393,11 +1393,11 @@ } // convert temperature excess to melt [kg] - for(int i=0;i(n); @@ -1434,8 +1434,8 @@ m[i] = m[i] - M[i]; // mass after melt Wi = (dIce-d[i]) * Swi * (m[i]/d[i]); // irreducible water - dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]); // change in pore water - R[i] = fmax(0.0, inM - dW[i]); // runoff + dW[i] = max(min(inM, Wi - W[i]),-1*W[i]); // change in pore water + R[i] = max(0.0, inM - dW[i]); // runoff F[i] = 0.0; } // check if no energy to refreeze meltwater @@ -1446,8 +1446,8 @@ m[i] = m[i] - M[i]; // mass after melt Wi = (dIce-d[i]) * Swi * (m[i]/d[i]); // irreducible water - dW[i] = fmax(fmin(inM, Wi - W[i]),-1*W[i]); // change in pore water - flxDn[i+1] = fmax(0.0, inM - dW[i]); // meltwater out + dW[i] = max(min(inM, Wi - W[i]),-1*W[i]); // change in pore water + flxDn[i+1] = max(0.0, inM - dW[i]); // meltwater out R[i] = 0.0; F[i] = 0.0; // no freeze } @@ -1459,13 +1459,13 @@ m[i] = m[i] - M[i]; IssmDouble dz_0 = m[i]/d[i]; IssmDouble dMax = (dIce - d[i])*dz_0; // d max = dIce - IssmDouble F1 = fmin(fmin(inM,dMax),maxF[i]); // maximum refreeze + IssmDouble F1 = min(min(inM,dMax),maxF[i]); // maximum refreeze m[i] = m[i] + F1; // mass after refreeze d[i] = m[i]/dz_0; //-----------------------pore water----------------------------- Wi = (dIce-d[i])* Swi * dz_0; // irreducible water - dW[i] = fmin(inM - F1, Wi-W[i]); // change in pore water + dW[i] = min(inM - F1, Wi-W[i]); // change in pore water if (dW[i] < 0.0-Wtol && -1*dW[i]>W[i]-Wtol ){ dW[i]= -1*W[i]; } @@ -1473,8 +1473,8 @@ if (dW[i] < 0.0-Wtol){ // excess pore water dMax = (dIce - d[i])*dz_0; // maximum refreeze - IssmDouble maxF2 = fmin(dMax, maxF[i]-F1); // maximum refreeze - F2 = fmin(-1*dW[i], maxF2); // pore water refreeze + IssmDouble maxF2 = min(dMax, maxF[i]-F1); // maximum refreeze + F2 = min(-1*dW[i], maxF2); // pore water refreeze m[i] = m[i] + F2; // mass after refreeze d[i] = m[i]/dz_0; dW[i] = dW[i] - F2; @@ -1752,7 +1752,7 @@ for(int i=0;i 0.0+Ttol){ Index: ../trunk-jpl/src/c/cores/controlvalidation_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 23241) +++ ../trunk-jpl/src/c/cores/controlvalidation_core.cpp (revision 23242) @@ -87,7 +87,7 @@ } /*output*/ - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ IssmPDouble* J_passive=xNew(2*num); for(int i=0;i<2*num;i++) J_passive[i]=reCast(output[i]); femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,J_passive,num,2,0,0)); Index: ../trunk-jpl/src/c/cores/control_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/control_core.cpp (revision 23241) +++ ../trunk-jpl/src/c/cores/control_core.cpp (revision 23242) @@ -96,7 +96,7 @@ if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run! femmodel->OutputControlsx(&femmodel->results); - #ifdef _HAVE_ADOLC_ + #ifdef _HAVE_AD_ IssmPDouble* J_passive=xNew(nsteps); for(int i=0;i(J[i]); femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,0,0));