1 | #!/usr/bin/env python3 |
---|
2 | # -*- Mode: Python; coding: utf-8; indent-tabs-mode: nil; tab-width: 4 -*- |
---|
3 | # |
---|
4 | ####################################################################################### |
---|
5 | # This script analyzes the output of STATION_ASF test-case with IDEALIZED forcing |
---|
6 | # in order to test the validity of computed fluxes and bulk transfer coefficient. |
---|
7 | # Beside an explicit standard output message, a result file is spawned: |
---|
8 | # * SBCBLK.success => the test passed ! |
---|
9 | # * SBCBLK.fail => the test failed ! |
---|
10 | # |
---|
11 | # Brodeau, 2020 |
---|
12 | # |
---|
13 | ######################################################################################## |
---|
14 | |
---|
15 | import sys |
---|
16 | from os import path as path |
---|
17 | import math |
---|
18 | import numpy as nmp |
---|
19 | from netCDF4 import Dataset |
---|
20 | |
---|
21 | l_t_shift = False ; # because time interp. is set to FALSE into "&namsbc_blk" of NEMO... |
---|
22 | # # ==> so time array is shifted by 30 minutes but fluxes are the same (persitence of input fields ???) |
---|
23 | |
---|
24 | l_alg = [ 'ECMWF' , 'NCAR' , 'COARE3p0', 'COARE3p6', 'ANDREAS' ] |
---|
25 | nb_alg = len(l_alg) |
---|
26 | |
---|
27 | |
---|
28 | # Variables to check: |
---|
29 | l_var_rf = [ 'Qsen' , 'Qlat' , 'Qlw' , 'Tau', 'Cd' , 'Ce' ] ; # In forcing file "IDEALIZED/input_output_VALIDATION_IDEALIZED.nc" |
---|
30 | l_var_ot = [ 'qsb_oce' , 'qla_oce' , 'qlw_oce', 'taum', 'Cd_oce', 'Ce_oce' ] ; # names in the NEMO/STATION_ASF output file (check file_def_oce.xml) |
---|
31 | nb_var = len(l_var_rf) |
---|
32 | |
---|
33 | dir_figs='.' |
---|
34 | size_fig=(13,8) |
---|
35 | fig_ext='png' |
---|
36 | |
---|
37 | |
---|
38 | rDPI=100. |
---|
39 | |
---|
40 | |
---|
41 | class fclrs: |
---|
42 | OKGR = '\033[92m' |
---|
43 | FAIL = '\033[91m' |
---|
44 | ENDC = '\033[0m' |
---|
45 | |
---|
46 | |
---|
47 | |
---|
48 | # Getting arguments: |
---|
49 | narg = len(sys.argv) |
---|
50 | if not narg in [3,4]: |
---|
51 | print('Usage: '+sys.argv[0]+' <forcing_+_validation_file> <NEMO-STATION_ASF_output_directory> (<m> for more/debug)'); sys.exit(0) |
---|
52 | cf_rf = sys.argv[1] |
---|
53 | cdir_out = sys.argv[2] |
---|
54 | l_more = False |
---|
55 | if narg==4: |
---|
56 | l_more = ( sys.argv[3] in ['m','M'] ) |
---|
57 | import matplotlib as mpl |
---|
58 | mpl.use('Agg') |
---|
59 | import matplotlib.pyplot as plt |
---|
60 | |
---|
61 | |
---|
62 | |
---|
63 | |
---|
64 | # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
65 | # Populating and checking existence of files to be read |
---|
66 | # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
67 | def chck4f(cf): |
---|
68 | cmesg = 'ERROR: File '+cf+' does not exist !!!' |
---|
69 | if not path.exists(cf): print(cmesg) ; sys.exit(0) |
---|
70 | |
---|
71 | print('\n') |
---|
72 | |
---|
73 | # Input forcing/valid file: |
---|
74 | chck4f(cf_rf) |
---|
75 | id_rf = Dataset(cf_rf) |
---|
76 | vt = id_rf.variables['time'][:] |
---|
77 | cunit_t = id_rf.variables['time'].units |
---|
78 | id_rf.close() |
---|
79 | Nt_rf = len(vt) |
---|
80 | vtime_rf = nmp.zeros(Nt_rf); vtime_rf[:] = vt[:] |
---|
81 | del vt |
---|
82 | print(' *** in forcing/valid file, "time" is in "'+cunit_t+'", Nt = '+str(Nt_rf)+'\n') |
---|
83 | |
---|
84 | |
---|
85 | # STATION_ASF output files (1 file per algorithm): |
---|
86 | cf_nemo = [] |
---|
87 | for ja in range(nb_alg): |
---|
88 | cfi = cdir_out+'/STATION_ASF-'+l_alg[ja]+'_IDEALIZED_1h_20200101_20200105_gridT.nc' |
---|
89 | chck4f(cfi) |
---|
90 | cf_nemo.append(cfi) |
---|
91 | print('\n *** NEMO/STATION_ASF output files we are goin to check:') |
---|
92 | for ja in range(nb_alg): print(cf_nemo[ja]) |
---|
93 | print('\n') |
---|
94 | #----------------------------------------------------------------- |
---|
95 | |
---|
96 | # Getting time array from the first file: |
---|
97 | id_nm = Dataset(cf_nemo[0]) |
---|
98 | vt = id_nm.variables['time_counter'][:] |
---|
99 | cunit_t = id_nm.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"') |
---|
100 | id_nm.close() |
---|
101 | Nt = len(vt) |
---|
102 | vtime_nm = nmp.zeros(Nt); vtime_nm[:] = vt[:] |
---|
103 | del vt |
---|
104 | |
---|
105 | if Nt != Nt_rf-1: print('ERROR: the two files do not agree in terms of record lengrth: '+Nt_rf-1+' vs '+Nt) |
---|
106 | |
---|
107 | print('\n *** Excellent! We are going to look at surface fluxes under '+str(Nt)+' different scenarios of air-sea stability/wind-speed conditions...\n') |
---|
108 | |
---|
109 | # 30 minute shift, just like NEMO |
---|
110 | vtime = nmp.zeros(Nt) |
---|
111 | if l_t_shift: |
---|
112 | vtime[:] = 0.5*(vtime_rf[:-1] + vtime_rf[1:]) |
---|
113 | else: |
---|
114 | vtime[:] = vtime_rf[:-1] |
---|
115 | # Debug: |
---|
116 | #for jt in range(3): |
---|
117 | # print(' Ref. , Nemo "', vtime_rf[jt], vtime_nm[jt], vtime[jt]) |
---|
118 | #sys.exit(0) |
---|
119 | |
---|
120 | |
---|
121 | |
---|
122 | |
---|
123 | ## |
---|
124 | IREPORT = nmp.zeros((nb_alg,nb_var), dtype=int) |
---|
125 | |
---|
126 | |
---|
127 | |
---|
128 | |
---|
129 | # Loop on the fields to control... |
---|
130 | ################################### |
---|
131 | |
---|
132 | jv=0 |
---|
133 | for cv in l_var_rf: |
---|
134 | cv_rf_m = cv+'_mean' |
---|
135 | cv_rf_t = cv+'_tol' |
---|
136 | cv_nemo = l_var_ot[jv] |
---|
137 | print('\n\n ==== Checking variable '+cv_nemo+' against '+cv_rf_m+' !') |
---|
138 | |
---|
139 | F_rf_m = nmp.zeros( Nt ) |
---|
140 | F_rf_t = nmp.zeros( Nt ) |
---|
141 | F_nemo = nmp.zeros((Nt,nb_alg)) |
---|
142 | |
---|
143 | wnd_rf = nmp.zeros( Nt ) ; #DEBUG |
---|
144 | tzt_rf = nmp.zeros( Nt ) ; #DEBUG |
---|
145 | qzt_rf = nmp.zeros( Nt ) ; #DEBUG |
---|
146 | sst_rf = nmp.zeros( Nt ) ; #DEBUG |
---|
147 | |
---|
148 | id_rf = Dataset(cf_rf) |
---|
149 | trf_m = id_rf.variables[cv_rf_m][:] ; # Nt+1 |
---|
150 | trf_t = id_rf.variables[cv_rf_t][:] ; # Nt+1 |
---|
151 | trf_wnd = id_rf.variables['wndspd'][:,1,1] ; # DEBUG |
---|
152 | trf_tzt = id_rf.variables['t_air'][:,1,1] ; # DEBUG |
---|
153 | trf_qzt = id_rf.variables['rh_air'][:,1,1] ; # DEBUG |
---|
154 | trf_sst = id_rf.variables['sst'][:,1,1] ; # DEBUG |
---|
155 | id_rf.close() |
---|
156 | |
---|
157 | if l_t_shift: |
---|
158 | # 30 minute shift, just like NEMO |
---|
159 | F_rf_m[:] = 0.5 * (trf_m[:-1] + trf_m[1:]) |
---|
160 | F_rf_t[:] = 0.5 * (trf_t[:-1] + trf_t[1:]) |
---|
161 | else: |
---|
162 | F_rf_m[:] = trf_m[:-1] |
---|
163 | F_rf_t[:] = trf_t[:-1] |
---|
164 | |
---|
165 | wnd_rf[:] = trf_wnd[:-1] |
---|
166 | tzt_rf[:] = trf_tzt[:-1] |
---|
167 | qzt_rf[:] = trf_qzt[:-1] |
---|
168 | sst_rf[:] = trf_sst[:-1] |
---|
169 | |
---|
170 | |
---|
171 | for ja in range(nb_alg): |
---|
172 | |
---|
173 | calgo = l_alg[ja] |
---|
174 | |
---|
175 | print(' *** '+calgo+' => '+cf_nemo[ja]) |
---|
176 | |
---|
177 | id_nemo = Dataset(cf_nemo[ja]) |
---|
178 | F_nemo[:,ja] = id_nemo.variables[cv_nemo][:,1,1] ; # it's 3x3 spatial domain, taking middle point ! |
---|
179 | id_nemo.close() |
---|
180 | |
---|
181 | |
---|
182 | if l_more: |
---|
183 | #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
184 | # How does all this look on on a figure? |
---|
185 | cfig = l_var_rf[jv]+'_'+calgo+'.'+fig_ext |
---|
186 | print(' *** will plot '+cfig) |
---|
187 | fig = plt.figure(num=1, figsize=size_fig, facecolor='w', edgecolor='k') |
---|
188 | ax1 = plt.axes([0.08, 0.25, 0.9, 0.7]) |
---|
189 | # |
---|
190 | plt.plot(vtime, F_nemo[:,ja], label='NEMO['+calgo+']', zorder=1) |
---|
191 | plt.plot(vtime, F_rf_m[:], color='k', label='MEAN REF!', zorder=10) |
---|
192 | # +- rtol enveloppe: |
---|
193 | plt.fill_between(vtime, F_rf_m[:]-F_rf_t[:], F_rf_m[:]+F_rf_t[:], alpha=0.2) |
---|
194 | # |
---|
195 | ax1.grid(color='k', linestyle='-', linewidth=0.3) |
---|
196 | plt.legend(loc='best', ncol=1, shadow=True, fancybox=True) |
---|
197 | plt.savefig(cfig, dpi=int(rDPI), transparent=False) |
---|
198 | plt.close(1) |
---|
199 | print('') |
---|
200 | #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
201 | |
---|
202 | |
---|
203 | # Does the field look okay with respect to reference +- tolerance ? |
---|
204 | |
---|
205 | l_overshoot_a = nmp.any( F_nemo[:,ja] > F_rf_m[:]+F_rf_t[:] ) |
---|
206 | l_overshoot_b = nmp.any( F_nemo[:,ja] < F_rf_m[:]-F_rf_t[:] ) |
---|
207 | |
---|
208 | if l_overshoot_a: print(fclrs.FAIL+'\n ***** BAD overshoot + for '+calgo+' for variable '+cv+' !!!\n'+fclrs.ENDC ) |
---|
209 | |
---|
210 | if l_overshoot_b: print(fclrs.FAIL+'\n ***** BAD overshoot - for '+calgo+' for variable '+cv+' !!!'+fclrs.ENDC ) |
---|
211 | |
---|
212 | if l_overshoot_a or l_overshoot_b: |
---|
213 | print(fclrs.FAIL+'\n ***** TEST NOT PASSED FOR '+calgo+' for variable '+cv+' !!!\n'+fclrs.ENDC ) |
---|
214 | else: |
---|
215 | # We're all good ! |
---|
216 | IREPORT[ja,jv] = 1 |
---|
217 | if l_more: print(fclrs.OKGR+'\n ***** TEST PASSED FOR '+calgo+' for variable '+cv+' :D !!!\n'+fclrs.ENDC ) |
---|
218 | |
---|
219 | jv=jv+1 |
---|
220 | |
---|
221 | |
---|
222 | |
---|
223 | l_ok = nmp.sum(IREPORT[:,:]) == nb_var*nb_alg |
---|
224 | |
---|
225 | |
---|
226 | |
---|
227 | if l_ok: |
---|
228 | ctxt = 'PASSED' |
---|
229 | ccol = fclrs.OKGR |
---|
230 | cf_report = 'SBCBLK.success' |
---|
231 | else: |
---|
232 | ctxt = 'FAILED' |
---|
233 | ccol = fclrs.FAIL |
---|
234 | cf_report = 'SBCBLK.fail' |
---|
235 | |
---|
236 | |
---|
237 | print(ccol+'\n\n ############ FINAL REPORT ############\n'+fclrs.ENDC) |
---|
238 | |
---|
239 | f = open(cf_report, 'w') |
---|
240 | |
---|
241 | f.write("### Sanity-check report for SBCBLK generated via 'STATION_ASF/EXP00/sbcblk_sanity_check.sh'\n\n") |
---|
242 | |
---|
243 | for ja in range(nb_alg): |
---|
244 | calgo = l_alg[ja] |
---|
245 | |
---|
246 | if nmp.sum(IREPORT[ja,:]) == nb_var: |
---|
247 | # Success for this algo |
---|
248 | cbla = ' ***** Algorithm "'+calgo+'" PASSED sanity check !!!\n\n' |
---|
249 | print(fclrs.OKGR+cbla+fclrs.ENDC ) ; f.write(cbla) |
---|
250 | else: |
---|
251 | # Algo FAILS! |
---|
252 | cbla = ' ***** Algorithm "'+calgo+'" FAILED sanity check !!!\n' |
---|
253 | print(fclrs.FAIL+cbla+fclrs.ENDC ) ; f.write(cbla) |
---|
254 | (idx_fail,) = nmp.where(IREPORT[ja,:]==0) |
---|
255 | for jv in idx_fail: |
---|
256 | cbla = ' ==> on variable '+l_var_rf[jv]+' !\n\n' |
---|
257 | print(fclrs.FAIL+cbla+fclrs.ENDC ) ; f.write(cbla) |
---|
258 | |
---|
259 | # Conclusion: |
---|
260 | clist='' |
---|
261 | for cc in l_var_ot: clist=clist+cc+', ' |
---|
262 | cbla = ' Test performed on the following NEMO prognostic variables:\n ==> '+clist[:-2]+'\n' |
---|
263 | print(ccol+cbla+fclrs.ENDC ) ; f.write(cbla+'\n') |
---|
264 | |
---|
265 | cbla = ' ####################################\n' +\ |
---|
266 | ' ### TEST '+ctxt+' FOR SBCBLK ! ###\n'+\ |
---|
267 | ' ####################################\n\n' |
---|
268 | |
---|
269 | print(ccol+cbla+fclrs.ENDC ) ; f.write(cbla) |
---|
270 | |
---|
271 | f.close() |
---|