-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathrun_rtm_enkf_steps.py
223 lines (160 loc) · 7.62 KB
/
run_rtm_enkf_steps.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
# First load modules!
# source ./modules.sh
'''
Assimilation of TBs to correct TOPAZ5 forecasts
Standalone analysis using the EnKF-C software
Model data: TOPAZ5 outputs
Observations: sea-ice concentration (SIT) and Tbs (19v, 19h, 37v, 37h) from AMSR2
Four data assimilation experiments are defined. Experiments are to be performed in this order:
exp_sic: synchronous assimilation of SIC
exp_tb: synchronous assimilation of Tbs
exp_sic_asyn: asynchronous assimilation of SIC
exp_tb_asyn: asynchronous assimilation of Tbs
This notebook is composed of five steps:
Computation of DAL and 2D-plane coefficients
RTM Tbs simulation
Plotting of RTM Tbs
Preparation of background ensemble, observations and mask for EnKF
Run EnKF
Plotting
'''
import importlib
from main_imports import cmd
config_filename = f'config_{config.date[0:4]}.py' # Name of you configuration file
cmd(f'cp {config_filename} config.py')
import config; importlib.reload(config)
# Model mask file to be used in the RTM simulations and DA analysis
cmd(f'rm {config.exps_dir}/conf/mask_topaz5.nc')
cmd(f'ln -s {config.exps_dir}/conf/{config.date[0:4]}/topaz5_grid_{config.date[0:4]}.nc {config.exps_dir}/conf/mask_topaz5.nc')
print("DA experiment: ", config.assim)
run_all_steps = False
if run_all_steps: compute_coeffs = run_rtm = make_plots_rtm = prepare_enkf = run_enkf = make_plots_da = True
else :
compute_coeffs = False
run_rtm = False
make_plots_rtm = False
prepare_enkf = False
run_enkf = False
make_plots_da = True
'''
1. Computation of DAL and 2D-plane coefficients
'''
if compute_coeffs :
from rtm_dal.main_fcts_rtm_tbs import *
print('Computation of DAL (Distance Along the Line) from TPD and TPA files...')
dal_norm, _, list_tpd_files = compute_dal()
print('\nComputation of 2D-plane coefficients. This 2D-plane is defined by the relation between Emissivity, DAL and T2m...')
compute_coeffs(dal_norm, list_tpd_files)
print('...computation finished.')
'''
2. RTM Tb simulation
'''
if run_rtm :
from rtm_dal.main_fcts_rtm_tbs import *
if 'exp_tb_asyn' in config.assim :
print('\nComputation of RTM Tbs (swaths)...')
import warnings
warnings.filterwarnings("ignore") # Ignore warning related to division in eq35
res = run_rtm_swaths(version = 1)
warnings.resetwarnings()
print('\nCreate directory where RTM Tbs will be saved')
cmd('mkdir -p ' + f"{config.rtm_tbs_dir}"); cmd('mkdir -p ' + f"{config.rtm_tbs_dir}/passes/")
print('\nCreate netCDF files containing RTM Tbs and saved them in the previously defined directory')
save_rtm_tbs(res[0], f"{config.rtm_tbs_dir}/passes/", swaths = True)
elif 'exp_sic' in config.assim :
print('\nComputation of RTM Tbs (daily means)')
import warnings
warnings.filterwarnings("ignore") # Ignore warning related to division in eq35
res = run_rtm(version = 1)
warnings.resetwarnings()
print('\nCreate directory where RTM Tbs will be saved')
cmd('mkdir -p ' + f"{config.rtm_tbs_dir}"); cmd('mkdir -p ' + f"{config.rtm_tbs_dir}/means/")
print('\nCreate netCDF files containing RTM Tbs and saved them in the previously defined directory')
save_rtm_tbs(res[0], f"{config.rtm_tbs_dir}/means/")
print('...RTM simulation finished.')
'''
3. Plotting from RTM Tb simulation (comparison to AMSR2 Tbs). Plot scripts might need to be adapted depending on your data format and dimensions
'''
if make_plots_rtm :
# Create a RTM Tbs map and difference to AMSR2 (two rows)
from plotting_codes import plotting_rtm
#from plotting_codes.plotting_rtm import *
importlib.reload(plotting_rtm)
plot_rtm = plotting_rtm.subplots_rtm_tbs #metrics
plot_rtm()
plot_hist = plotting_rtm.plot_histograms
plot_hist()
plot_diags = plotting_rtm.plot_diagrams
plot_diags()
print('...plotting finished.')
'''
4. Preparation of background ensemble, observations and mask for EnKF
'''
if prepare_enkf :
print('Preparation of background ensemble, observations and mask for EnKF-C')
cmd('mkdir -p ' + config.storage_dir)
cmd('mkdir -p ' + config.storage_dir + 'ensb');
cmd('rm ' + config.storage_dir + 'ensb/*nc')
if config.assim == 'exp_sic' :
from steps_da import prepare_ens; importlib.reload(prepare_ens); prep_ens = prepare_ens.prep_ensemble
from steps_da import prepare_obs; importlib.reload(prepare_obs); prep_topaz = prepare_obs.prep_topaz
from steps_da import model_mask; importlib.reload(model_mask); prep_mask = model_mask.generate_mask
print('Ensemble preparation...'); prep_ens()
print('Prepare observations (means)...'); prep_topaz()
print('Generation of TOPAZ mask...'); prep_mask()
elif config.assim == 'exp_tb' :
storage_dir_tbs = f"{config.exps_dir}/exps_2021/exp_sic/{config.date}/"
cmd('ln -s ' + storage_dir_tbs + 'ensb/* ' + config.storage_dir + 'ensb/') # Link to background ensemble
elif config.assim == 'exp_sic_asyn' :
storage_dir_asyn = f"{config.exps_dir}/exps_2021/exp_sic/{config.date}/"
cmd('ln -s ' + storage_dir_asyn + 'ensb/* ' + config.storage_dir + 'ensb/') # Link to background ensemble
from steps_da import prepare_ens; importlib.reload(prepare_ens); prep_ens = prepare_ens.prep_ensemble_asyn
from steps_da import prepare_obs; importlib.reload(prepare_obs); prep_topaz = prepare_obs.prep_topaz_passes
from steps_da import change_date; importlib.reload(change_date); update = change_date.update_enkf_prm
print('Ensemble preparation for ASYN experiment...'); prep_ens()
print('Prepare observations (passes)...'); prep_topaz()
print('Update the EnKF date...'); update()
elif config.assim == 'exp_tb_asyn' :
storage_dir_asyn = f"{config.exps_dir}/exps_2021/exp_sic/{config.date}/"
cmd('ln -s ' + storage_dir_asyn + 'ensb/* ' + config.storage_dir + 'ensb/') # Link to background ensemble
from steps_da import prepare_ens; importlib.reload(prepare_ens); prep_ens = prepare_ens.prep_ensemble_asyn
print('Ensemble preparation for ASYN experiment...'); prep_ens()
print('...preparation of model and observation data finished.')
'''
5. Run EnKF
'''
if run_enkf :
cmd('module use /modules/MET/rhel8/user-modules/')
cmd('module load enkfc/2.9.9')
print('Running EnKF-C...')
from steps_da.run_da import run_enkf
run_enkf()
print('...run finished.')
'''
6. Plotting from DA results
vari = 0 is aice
vari = 1 is hi
vari = 8 is iage
'''
vari = 0
metrics = True
if make_plots_da :
print(f'Making figures for variable {config.varnames[vari]}...')
from plotting_codes import plotting_da
importlib.reload(plotting_da)
from plotting_codes.plotting_da import *
plot_metrics = plotting_da.plot_metrics
if config.assim == 'exp_sic' :
if vari == 0 :
background_maps()
background_maps(var = 'tb', var_tb = 0)
analysis_maps(vari)
elif vari > 0 : no_obs_maps(vari)
else :
if vari == 0 :
analysis_maps(vari)
elif vari > 0 : no_obs_maps(vari)
if metrics :
print('\nPlotting EnKF metrics (DFS and SRF)...')
plot_metrics()
print('...plotting finished.')