[13676] | 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 |
---|
[13682] | 19 | from netCDF4 import Dataset |
---|
[13676] | 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'] ) |
---|
[13682] | 57 | import matplotlib as mpl |
---|
| 58 | mpl.use('Agg') |
---|
| 59 | import matplotlib.pyplot as plt |
---|
[13676] | 60 | |
---|
| 61 | |
---|
[13682] | 62 | |
---|
[13676] | 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] |
---|
[13686] | 137 | print('\n\n ==== Checking variable '+cv_nemo+' against '+cv_rf_m+' !') |
---|
[13676] | 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: |
---|
[13733] | 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') |
---|
[13676] | 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() |
---|