Skip to content

Commit 947b397

Browse files
Merge branch 'develop' into feature/gdas-validation
2 parents 654dcdf + 50f5ebc commit 947b397

6 files changed

+67
-185
lines changed

CMakeLists.txt

+1-1
Original file line numberDiff line numberDiff line change
@@ -58,7 +58,7 @@ if(BUILD_GDASBUNDLE)
5858

5959
ecbuild_bundle( PROJECT eckit GIT "https://github.com/ecmwf/eckit.git" TAG 1.16.0 )
6060
ecbuild_bundle( PROJECT fckit GIT "https://github.com/ecmwf/fckit.git" TAG 0.9.2 )
61-
ecbuild_bundle( PROJECT atlas GIT "https://github.com/ecmwf/atlas.git" TAG 0.29.0 )
61+
ecbuild_bundle( PROJECT atlas GIT "https://github.com/ecmwf/atlas.git" TAG 0.35.0 )
6262

6363
# External (required) observation operators
6464
option("BUNDLE_SKIP_CRTM" "Don't build CRTM" "OFF") # Don't build crtm unless user passes -DBUNDLE_SKIP_CRTM=OFF

parm/soca/berror/soca_ensb.yaml

+3-1
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,9 @@ vertical geometry:
1616
ocn_filename: MOM.res.nc
1717
read_from_file: 3
1818

19-
soca increments:
19+
add recentering increment: false
20+
21+
soca increments: # Could also be states, but they are read as increments
2022
number of increments: ${ENS_SIZE}
2123
pattern: '%mem%'
2224
template:

scripts/exgdas_global_marine_analysis_bmat.sh

+12-92
Original file line numberDiff line numberDiff line change
@@ -36,36 +36,6 @@ pwd=$(pwd)
3636
# Utilities
3737
export NLN=${NLN:-"/bin/ln -sf"}
3838

39-
40-
function bump_vars()
41-
{
42-
tmp=$(ls -d ${1}_* )
43-
lov=()
44-
for v in $tmp; do
45-
lov[${#lov[@]}]=${v#*_}
46-
done
47-
echo "$lov"
48-
}
49-
50-
function concatenate_bump()
51-
{
52-
bumpdim=$1
53-
# Concatenate the bump files
54-
vars=$(bump_vars $bumpdim)
55-
n=$(wc -w <<< "$vars")
56-
echo "concatenating $n variables: $vars"
57-
lof=$(ls ${bumpdim}_${vars[0]})
58-
echo $lof
59-
for f in $lof; do
60-
bumpbase=${f#*_}
61-
output=bump/${bumpdim}_$bumpbase
62-
lob=$(ls ${bumpdim}_*/*$bumpbase)
63-
for b in $lob; do
64-
ncks -A $b $output
65-
done
66-
done
67-
}
68-
6939
function clean_yaml()
7040
{
7141
mv $1 tmp_yaml;
@@ -86,23 +56,6 @@ else
8656
fi
8757
fi
8858

89-
################################################################################
90-
# Prepare the diagonal of the parametric B
91-
shopt -s nullglob
92-
files=(ocn.bkgerr_stddev.incr.*.nc)
93-
echo $files
94-
if [ ${#files[@]} -gt 0 ]; then
95-
echo "Diag of B already staged, skipping the parametric diag B"
96-
else
97-
# TODO: this step should be replaced by a post processing step of the ens. derived std. dev.
98-
$APRUN_OCNANAL $JEDI_BIN/soca_convertincrement.x parametric_stddev_b.yaml > parametric_stddev_b.out 2>&1
99-
export err=$?; err_chk
100-
if [ $err -gt 0 ]; then
101-
exit $err
102-
fi
103-
fi
104-
shopt -u nullglob
105-
10659
################################################################################
10760
# Write ensemble weights for the hybrid envar
10861
$APRUN_OCNANAL $JEDI_BIN/gdas_socahybridweights.x soca_ensweights.yaml
@@ -112,16 +65,18 @@ if [ $err -gt 0 ]; then
11265
fi
11366

11467
################################################################################
115-
# Compute the ens std. dev, ignore ssh variance
116-
# TODO (G): Implement what's below into one single oops application
117-
# 0 - Compute moments of original ensemble, used at the diag of the first
118-
# component of the hybrid B
119-
# 1 - Ensemble perturbations:
120-
# o Vertically Interpolate to the deterministic layers
121-
# o Recenter around 0 to create an ensemble of perurbations
122-
# 2 - Filter members and apply the linear steric height balance to each members
123-
# 3 - Copy h from deterministic to unbalanced perturbations
124-
# 4 - Compute moments of converted ensemble perturbations
68+
# Ensemble perturbations for the EnVAR and diagonal of static B
69+
# Static B:
70+
# - Compute moments of original ensemble with the balanced ssh removed
71+
# from the statistics
72+
#
73+
# Ensemble of perturbations:
74+
# - apply the linear steric height balance to each members, using the
75+
# deterministic for trajectory.
76+
# - add the unblanced ssh to the steric ssh field
77+
# Diagnostics:
78+
# - variance explained by the steric heigh
79+
# - moments
12580

12681
# Process static ensemble
12782
clean_yaml soca_clim_ens_moments.yaml
@@ -146,14 +101,6 @@ if [ ${#files[@]} -gt 0 ]; then
146101
shopt -u nullglob
147102
fi
148103

149-
################################################################################
150-
# Set decorrelation scales for bump C
151-
$APRUN_OCNANAL $JEDI_BIN/soca_setcorscales.x soca_setcorscales.yaml
152-
export err=$?; err_chk
153-
if [ $err -gt 0 ]; then
154-
exit $err
155-
fi
156-
157104
################################################################################
158105
# Set localization scales for the hybrid en. var.
159106
$APRUN_OCNANAL $JEDI_BIN/soca_setcorscales.x soca_setlocscales.yaml
@@ -162,33 +109,6 @@ if [ $err -gt 0 ]; then
162109
exit $err
163110
fi
164111

165-
################################################################################
166-
# Compute convolution coefs for C
167-
# TODO (G, C, R, ...): problem with ' character when reading yaml, removing from file for now
168-
# 2D C from bump
169-
if false; then
170-
# TODO: resurect this section when making use of bump 3D in the static B, skip for now
171-
yaml_bump2d=soca_bump2d.yaml
172-
clean_yaml $yaml_bump2d
173-
$APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x $yaml_bump2d 2>$yaml_bump2d.err
174-
export err=$?; err_chk
175-
if [ $err -gt 0 ]; then
176-
exit $err
177-
fi
178-
179-
# 3D C from bump
180-
yaml_list=`ls soca_bump3d_*.yaml`
181-
for yaml in $yaml_list; do
182-
clean_yaml $yaml
183-
$APRUN_OCNANAL $JEDI_BIN/soca_error_covariance_toolbox.x $yaml 2>$yaml.err
184-
export err=$?; err_chk
185-
if [ $err -gt 0 ]; then
186-
exit $err
187-
fi
188-
done
189-
concatenate_bump 'bump3d'
190-
fi
191-
192112
################################################################################
193113
# Compute convolution coefs for L
194114
clean_yaml soca_bump_loc.yaml

scripts/exgdas_global_marine_analysis_prep.py

+1-61
Original file line numberDiff line numberDiff line change
@@ -308,16 +308,6 @@ def find_clim_ens(input_date):
308308
logging.info(f"---------------- Stage static files")
309309
ufsda.stage.soca_fix(stage_cfg)
310310

311-
################################################################################
312-
# stage background error files
313-
314-
logging.info(f"---------------- Stage static files")
315-
bkgerr_list = []
316-
for domain in ['ocn', 'ice']:
317-
fname_stddev = find_bkgerr(pytz.utc.localize(window_begin, is_dst=None), domain=domain)
318-
fname_out = domain+'.bkgerr_stddev.incr.'+window_begin_iso+'.nc'
319-
bkgerr_list.append([fname_stddev, fname_out])
320-
FileHandler({'copy': bkgerr_list}).sync()
321311

322312
################################################################################
323313
# stage ensemble members
@@ -389,14 +379,6 @@ def find_clim_ens(input_date):
389379
config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get)
390380
config.save(berr_yaml)
391381

392-
################################################################################
393-
# copy yaml for decorrelation length scales
394-
395-
logging.info(f"---------------- generate soca_setcorscales.yaml")
396-
corscales_yaml_src = os.path.join(gdas_home, 'parm', 'soca', 'berror', 'soca_setcorscales.yaml')
397-
corscales_yaml_dst = os.path.join(stage_cfg['stage_dir'], 'soca_setcorscales.yaml')
398-
FileHandler({'copy': [[corscales_yaml_src, corscales_yaml_dst]]}).sync()
399-
400382
################################################################################
401383
# copy yaml for localization length scales
402384

@@ -408,49 +390,7 @@ def find_clim_ens(input_date):
408390
################################################################################
409391
# generate yaml for bump/nicas (used for correlation and/or localization)
410392

411-
logging.info(f"---------------- generate BUMP/NICAS yamls")
412-
# TODO (Guillaume): move the possible vars somewhere else
413-
vars3d = ['tocn', 'socn', 'uocn', 'vocn', 'chl', 'biop']
414-
vars2d = ['ssh', 'cicen', 'hicen', 'hsnon', 'swh',
415-
'sw', 'lw', 'lw_rad', 'lhf', 'shf', 'us']
416-
417-
# 2d bump yaml (all 2d vars at once)
418-
bumpdir = 'bump'
419-
ufsda.disk_utils.mkdir(os.path.join(anl_dir, bumpdir))
420-
bump_yaml = os.path.join(anl_dir, 'soca_bump2d.yaml')
421-
bump_yaml_template = os.path.join(gdas_home,
422-
'parm',
423-
'soca',
424-
'berror',
425-
'soca_bump2d.yaml')
426-
config = YAMLFile(path=bump_yaml_template)
427-
config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get)
428-
config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get)
429-
config.save(bump_yaml)
430-
431-
# 3d bump yaml, 1 yaml per variable
432-
soca_vars = ['tocn', 'socn', 'uocn', 'vocn']
433-
for v in soca_vars:
434-
logging.info(f"creating the yaml to initialize bump for {v}")
435-
if v in vars2d:
436-
continue
437-
else:
438-
dim = '3d'
439-
bump_yaml = os.path.join(anl_dir, 'soca_bump'+dim+'_'+v+'.yaml')
440-
bump_yaml_template = os.path.join(gdas_home,
441-
'parm',
442-
'soca',
443-
'berror',
444-
'soca_bump_split.yaml')
445-
bumpdir = 'bump'+dim+'_'+v
446-
os.environ['BUMPDIR'] = bumpdir
447-
ufsda.disk_utils.mkdir(os.path.join(anl_dir, bumpdir))
448-
os.environ['CVAR'] = v
449-
config = YAMLFile(path=bump_yaml_template)
450-
config = Template.substitute_structure(config, TemplateConstants.DOUBLE_CURLY_BRACES, envconfig.get)
451-
config = Template.substitute_structure(config, TemplateConstants.DOLLAR_PARENTHESES, envconfig.get)
452-
config.save(bump_yaml)
453-
393+
logging.info(f"---------------- generate BUMP/NICAS localization yamls")
454394
# localization bump yaml
455395
bumpdir = 'bump'
456396
ufsda.disk_utils.mkdir(os.path.join(anl_dir, bumpdir))

utils/soca/gdas_ens_handler.h

+38-20
Original file line numberDiff line numberDiff line change
@@ -97,6 +97,12 @@ namespace gdasapp {
9797
ensMembers.push_back(postProcIncr.read(i));
9898
}
9999

100+
// Check if we need to recenter the increment around the deterministic
101+
bool addRecenterIncr(false);
102+
if ( fullConfig.has("add recentering increment") ) {
103+
fullConfig.get("add recentering increment", addRecenterIncr);
104+
}
105+
100106
// Compute ensemble moments
101107
soca::Increment ensMean(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
102108
soca::Increment ensStd(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
@@ -120,24 +126,25 @@ namespace gdasapp {
120126
fullConfig.get("steric height", stericVarChangeConfig);
121127
oops::Log::info() << "steric config 0000: " << stericVarChangeConfig << std::endl;
122128

123-
// Initialize trajectories
129+
// Initialize the trajectories used in the linear variable changes
124130
const eckit::LocalConfiguration trajConfig(fullConfig, "trajectory");
125-
soca::State cycleTraj(geom, trajConfig); // trajectory of the cycle
126-
soca::State meanTraj(cycleTraj); // trajectory defined as the ens. mean
127-
meanTraj.zero();
128-
meanTraj += ensMean;
129-
130-
// Compute the error between the ensemble mean and the deterministic
131-
soca::Increment deterministicError(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
132-
deterministicError.diff(cycleTraj, meanTraj);
131+
soca::State determTraj(geom, trajConfig); // deterministic trajectory
132+
soca::State ensMeanTraj(determTraj); // trajectory defined as the ens. mean
133+
ensMeanTraj.zero();
134+
ensMeanTraj += ensMean;
135+
136+
// Compute the recentering increment as the difference between
137+
// the ensemble mean and the deterministic
138+
soca::Increment recenteringIncr(geom, postProcIncr.socaIncrVar_, postProcIncr.dt_);
139+
recenteringIncr.diff(determTraj, ensMeanTraj);
133140
eckit::LocalConfiguration sshRecErrOutputConfig(fullConfig, "ssh output.recentering error");
134-
deterministicError = postProcIncr.setToZero(deterministicError);
141+
recenteringIncr = postProcIncr.setToZero(recenteringIncr);
135142
oops::Log::info() << "steric config : " << stericVarChangeConfig << std::endl;
136-
postProcIncr.applyLinVarChange(deterministicError, stericVarChangeConfig, cycleTraj);
137-
oops::Log::info() << "ensemble mean: " << meanTraj << std::endl;
138-
oops::Log::info() << "deterministic: " << cycleTraj << std::endl;
139-
oops::Log::info() << "error: " << deterministicError << std::endl;
140-
deterministicError.write(sshRecErrOutputConfig);
143+
postProcIncr.applyLinVarChange(recenteringIncr, stericVarChangeConfig, determTraj);
144+
oops::Log::info() << "ensemble mean: " << ensMeanTraj << std::endl;
145+
oops::Log::info() << "deterministic: " << determTraj << std::endl;
146+
oops::Log::info() << "error: " << recenteringIncr << std::endl;
147+
recenteringIncr.write(sshRecErrOutputConfig);
141148

142149
// Re-process the ensemble of perturbations
143150
int result = 0;
@@ -165,7 +172,7 @@ namespace gdasapp {
165172

166173
// Compute the original steric height perturbation from T and S
167174
eckit::LocalConfiguration stericConfig(fullConfig, "steric height");
168-
postProcIncr.applyLinVarChange(incr, stericConfig, meanTraj);
175+
postProcIncr.applyLinVarChange(incr, stericConfig, ensMeanTraj);
169176
ssh_tmp = incr;
170177
sshSteric.push_back(ssh_tmp);
171178

@@ -183,10 +190,10 @@ namespace gdasapp {
183190
incr = postProcIncr.setToZero(incr);
184191

185192
// Filter ensemble member and recompute steric ssh, recentering around
186-
// the cycle's trajectory
193+
// the deterministic trajectory
187194
if ( fullConfig.has("linear variable change") ) {
188195
eckit::LocalConfiguration lvcConfig(fullConfig, "linear variable change");
189-
postProcIncr.applyLinVarChange(incr, lvcConfig, cycleTraj);
196+
postProcIncr.applyLinVarChange(incr, lvcConfig, determTraj);
190197
}
191198

192199
// Add the unbalanced ssh to the recentered perturbation
@@ -200,12 +207,23 @@ namespace gdasapp {
200207
incr.fromFieldSet(incrFs);
201208
oops::Log::debug() << "&&&&& after adding ssh_u " << incr << std::endl;
202209

203-
// Save final perturbation, used in the offline EnVAR
210+
// Save final perturbation, used in the EnVAR or to initialize the ensemble forecast
211+
if (addRecenterIncr) {
212+
// Add the recentering increment to the increment ensemble members.
213+
// This is generally used to recenter the ensemble fcst around the
214+
// deterministic
215+
incr += recenteringIncr;
216+
}
204217
result = postProcIncr.save(incr, i+1);
205218

206219
// Update ensemble
207220
ensMembers[i] = incr;
208-
}
221+
}
222+
223+
// Exit early if only recentering around the deterministic
224+
if (addRecenterIncr) {
225+
return 0;
226+
}
209227

210228
// Compute ensemble moments for total ssh
211229
soca::Increment sshMean(geom, socaSshVar, postProcIncr.dt_);

utils/test/CMakeLists.txt

+12-10
Original file line numberDiff line numberDiff line change
@@ -62,11 +62,12 @@ ecbuild_add_test( TARGET test_gdasapp_util_ghrsst2ioda
6262
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
6363

6464
# Test the SMAP to IODA converter
65-
ecbuild_add_test( TARGET test_gdasapp_util_smap2ioda
66-
COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x
67-
ARGS "../testinput/gdas_smap2ioda.yaml"
68-
LIBS gdas-utils
69-
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
65+
# TODO(Mindo): Turn back on when date is fixed
66+
#ecbuild_add_test( TARGET test_gdasapp_util_smap2ioda
67+
# COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x
68+
# ARGS "../testinput/gdas_smap2ioda.yaml"
69+
# LIBS gdas-utils
70+
# WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
7071

7172
# Test the SMOS to IODA converter
7273
ecbuild_add_test( TARGET test_gdasapp_util_smos2ioda
@@ -76,8 +77,9 @@ ecbuild_add_test( TARGET test_gdasapp_util_smos2ioda
7677
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
7778

7879
# Test the AMSR2 to IODA converter
79-
ecbuild_add_test( TARGET test_gdasapp_util_icecamsr2ioda
80-
COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x
81-
ARGS "../testinput/gdas_icecamsr2ioda.yaml"
82-
LIBS gdas-utils
83-
WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)
80+
# TODO(Mindo): Turn back on when date is fixed
81+
#ecbuild_add_test( TARGET test_gdasapp_util_icecamsr2ioda
82+
# COMMAND ${CMAKE_BINARY_DIR}/bin/gdas_obsprovider2ioda.x
83+
# ARGS "../testinput/gdas_icecamsr2ioda.yaml"
84+
# LIBS gdas-utils
85+
# WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/obsproc)

0 commit comments

Comments
 (0)