Skip to content

Commit 0f9dc3a

Browse files
author
Cristian Lussana
committed
fixed bug in allow_extrapolation
1 parent c104db0 commit 0f9dc3a

File tree

1 file changed

+37
-19
lines changed

1 file changed

+37
-19
lines changed

‎src/api/oi_ensi_lr.cpp‎

Lines changed: 37 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -157,6 +157,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
157157
throw std::invalid_argument("Background RIGTH and points size mismatch");
158158

159159
float hmax = 7 * dh;
160+
float default_min_std = 0.0013;
160161
/*
161162
BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax)t;
162163
if(which_structfun == 0) {
@@ -228,7 +229,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
228229
if(nValidEns == 0)
229230
return background_l;
230231

231-
// gZ_R(nY, nValidEns): used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observationsCompute Y
232+
// gZ_R(nY, nValidEns): used to compute ensemble-based background correlations i) between yth gridpoint and observations ii) among observations
232233
vec2 gZ_R = gridpp::init_vec2(nY, nValidEns); // useful to compute dynamical correlations
233234
for(int i = 0; i < nS; i++) {
234235
vec pbackgroundValid_R(nValidEns);
@@ -243,7 +244,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
243244
gZ_R[i][e] = 0;
244245
}
245246
else {
246-
if(std == 0) {
247+
if(std <= default_min_std) {
247248
for(int e = 0; e < nValidEns; e++)
248249
gZ_R[i][e] = 0;
249250
}
@@ -267,11 +268,13 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
267268
/* float localizationRadius = structure.localization_distance(p1); */
268269
float localizationRadius = 0;
269270
if(which_structfun == 0) {
270-
BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);
271+
/* BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);*/
272+
BarnesStructure structure = BarnesStructure( dh, dz, dw);
271273
localizationRadius = structure.localization_distance(p1);
272274
}
273275
else if(which_structfun == 1) {
274-
MixAStructure structure = MixAStructure( dh, dz, dw, hmax);
276+
/* MixAStructure structure = MixAStructure( dh, dz, dw, hmax); */
277+
MixAStructure structure = MixAStructure( dh, dz, dw);
275278
localizationRadius = structure.localization_distance(p1);
276279
}
277280

@@ -295,11 +298,13 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
295298
}
296299
vec rhos(lLocIndices0.size());
297300
if(which_structfun == 0) {
298-
BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);
301+
/* BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); */
302+
BarnesStructure structure = BarnesStructure( dh, dz, dw);
299303
rhos = structure.corr_background(p1, p2);
300304
}
301305
else if(which_structfun == 1) {
302-
MixAStructure structure = MixAStructure( dh, dz, dw, hmax);
306+
/* MixAStructure structure = MixAStructure( dh, dz, dw, hmax); */
307+
MixAStructure structure = MixAStructure( dh, dz, dw);
303308
rhos = structure.corr_background(p1, p2);
304309
}
305310
/* vec rhos = structure.corr_background(p1, p2); */
@@ -358,7 +363,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
358363
}
359364
float mean = gridpp::calc_statistic(backgroundValid_L, gridpp::Mean);
360365
float std = gridpp::calc_statistic(backgroundValid_L, gridpp::Std);
361-
if(gridpp::is_valid(mean) && gridpp::is_valid(std) && std != 0) {
366+
if(gridpp::is_valid(mean) && gridpp::is_valid(std) && std > default_min_std) {
362367
for(int e = 0; e < nValidEns; e++)
363368
lX_L(0,e) = 1 / sqrt(nValidEns-1) * (backgroundValid_L[e] - mean) / std;
364369
}
@@ -393,11 +398,13 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
393398
}
394399
vec corr(lS);
395400
if(which_structfun == 0) {
396-
BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);
401+
/* BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); */
402+
BarnesStructure structure = BarnesStructure( dh, dz, dw);
397403
corr = structure.corr_background(p1, p2);
398404
}
399405
else if(which_structfun == 1) {
400-
MixAStructure structure = MixAStructure( dh, dz, dw, hmax);
406+
/* MixAStructure structure = MixAStructure( dh, dz, dw, hmax); */
407+
MixAStructure structure = MixAStructure( dh, dz, dw);
401408
corr = structure.corr_background(p1, p2);
402409
}
403410
/* vec corr = structure.corr(p1, p2); */
@@ -429,13 +436,15 @@ vec2 gridpp::R_optimal_interpolation_ensi_lr(const gridpp::Points& bpoints,
429436
increment = maxInc;
430437
}
431438
else if(maxInc < 0 && increment > 0) {
432-
increment = maxInc;
439+
/* increment = maxInc; */
440+
increment = 0;
433441
}
434442
else if(minInc < 0 && increment < minInc) {
435443
increment = minInc;
436444
}
437445
else if(minInc > 0 && increment < 0) {
438-
increment = minInc;
446+
/* increment = minInc; */
447+
increment = 0;
439448
}
440449
dx[e] = increment;
441450
}
@@ -491,6 +500,7 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp
491500
throw std::invalid_argument("Background rigth and points size mismatch");
492501

493502
float hmax = 7 * dh;
503+
float default_min_std = 0.0013;
494504

495505
int nS = points.size();
496506
if(nS == 0)
@@ -560,11 +570,13 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp
560570
/* float localizationRadius = structure.localization_distance(p1); */
561571
float localizationRadius = 0;
562572
if(which_structfun == 0) {
563-
BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);
573+
/* BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); */
574+
BarnesStructure structure = BarnesStructure( dh, dz, dw);
564575
localizationRadius = structure.localization_distance(p1);
565576
}
566577
else if(which_structfun == 1) {
567-
MixAStructure structure = MixAStructure( dh, dz, dw, hmax);
578+
/* MixAStructure structure = MixAStructure( dh, dz, dw, hmax); */
579+
MixAStructure structure = MixAStructure( dh, dz, dw);
568580
localizationRadius = structure.localization_distance(p1);
569581
}
570582

@@ -588,11 +600,13 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp
588600
}
589601
vec rhos(lLocIndices0.size());
590602
if(which_structfun == 0) {
591-
BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);
603+
/* BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); */
604+
BarnesStructure structure = BarnesStructure( dh, dz, dw);
592605
rhos = structure.corr_background(p1, p2);
593606
}
594607
else if(which_structfun == 1) {
595-
MixAStructure structure = MixAStructure( dh, dz, dw, hmax);
608+
/* MixAStructure structure = MixAStructure( dh, dz, dw, hmax); */
609+
MixAStructure structure = MixAStructure( dh, dz, dw);
596610
rhos = structure.corr_background(p1, p2);
597611
}
598612
for(int i = 0; i < lLocIndices0.size(); i++) {
@@ -669,11 +683,13 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp
669683
}
670684
vec corr(lS);
671685
if(which_structfun == 0) {
672-
BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax);
686+
/* BarnesStructure structure = BarnesStructure( dh, dz, dw, hmax); */
687+
BarnesStructure structure = BarnesStructure( dh, dz, dw);
673688
corr = structure.corr_background(p1, p2);
674689
}
675690
else if(which_structfun == 1) {
676-
MixAStructure structure = MixAStructure( dh, dz, dw, hmax);
691+
/* MixAStructure structure = MixAStructure( dh, dz, dw, hmax); */
692+
MixAStructure structure = MixAStructure( dh, dz, dw);
677693
corr = structure.corr_background(p1, p2);
678694
}
679695
/* vec corr = structure.corr(p1, p2); */
@@ -699,13 +715,15 @@ vec2 gridpp::R_optimal_interpolation_ensi_staticcorr_lr(const gridpp::Points& bp
699715
increment = maxInc;
700716
}
701717
else if(maxInc < 0 && increment > 0) {
702-
increment = maxInc;
718+
/* increment = maxInc; */
719+
increment = 0;
703720
}
704721
else if(minInc < 0 && increment < minInc) {
705722
increment = minInc;
706723
}
707724
else if(minInc > 0 && increment < 0) {
708-
increment = minInc;
725+
/* increment = minInc; */
726+
increment = 0;
709727
}
710728
dx[e] = increment;
711729
}

0 commit comments

Comments
 (0)