source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/STATION_ASF/EXPREF/analyze_output.py @ 13682

Last change on this file since 13682 was 13682, checked in by laurent, 7 months ago

Still improving the oce+ice capable STATION_ASF…

  • Property svn:executable set to *
File size: 8.4 KB
Line 
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
15import sys
16from os import path as path
17import math
18import numpy as nmp
19from netCDF4 import Dataset
20
21l_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
24l_alg = [ 'ECMWF' , 'NCAR' , 'COARE3p0', 'COARE3p6', 'ANDREAS' ]
25nb_alg     =    len(l_alg)
26
27
28# Variables to check:
29l_var_rf = [  'Qsen'   ,   'Qlat'   ,  'Qlw'   ,   'Tau',   'Cd'  ,   'Ce'   ] ; # In forcing file "IDEALIZED/input_output_VALIDATION_IDEALIZED.nc"
30l_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)
31nb_var  = len(l_var_rf)
32
33dir_figs='.'
34size_fig=(13,8)
35fig_ext='png'
36
37
38rDPI=100.
39
40
41class fclrs:
42    OKGR = '\033[92m'
43    FAIL = '\033[91m'
44    ENDC = '\033[0m'
45
46
47
48# Getting arguments:
49narg = len(sys.argv)
50if 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)
52cf_rf = sys.argv[1]
53cdir_out = sys.argv[2]
54l_more = False
55if 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# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
67def chck4f(cf):
68    cmesg = 'ERROR: File '+cf+' does not exist !!!'
69    if not path.exists(cf): print(cmesg) ; sys.exit(0)
70
71print('\n')
72
73# Input forcing/valid file:
74chck4f(cf_rf)
75id_rf = Dataset(cf_rf)
76vt = id_rf.variables['time'][:]
77cunit_t = id_rf.variables['time'].units
78id_rf.close()
79Nt_rf = len(vt)
80vtime_rf = nmp.zeros(Nt_rf); vtime_rf[:]  = vt[:]
81del vt
82print(' *** 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):
86cf_nemo = []
87for 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)
91print('\n *** NEMO/STATION_ASF output files we are goin to check:')
92for ja in range(nb_alg): print(cf_nemo[ja])
93print('\n')
94#-----------------------------------------------------------------
95
96# Getting time array from the first file:
97id_nm = Dataset(cf_nemo[0])
98vt = id_nm.variables['time_counter'][:]
99cunit_t = id_nm.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
100id_nm.close()
101Nt = len(vt)
102vtime_nm = nmp.zeros(Nt); vtime_nm[:]  = vt[:]
103del vt
104
105if Nt != Nt_rf-1: print('ERROR: the two files do not agree in terms of record lengrth: '+Nt_rf-1+' vs '+Nt)
106
107print('\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
110vtime = nmp.zeros(Nt)
111if l_t_shift:
112    vtime[:] = 0.5*(vtime_rf[:-1] + vtime_rf[1:])
113else:
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##
124IREPORT = nmp.zeros((nb_alg,nb_var), dtype=int)
125
126
127
128
129# Loop on the fields to control...
130###################################
131
132jv=0
133for 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 flux '+cv_rf_m+' ('+cv_nemo+') !')
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
223l_ok =  nmp.sum(IREPORT[:,:]) == nb_var*nb_alg
224
225
226
227if l_ok:
228    ctxt      = 'PASSED'
229    ccol      = fclrs.OKGR
230    cf_report = 'SBCBLK.success'
231else:
232    ctxt      = 'FAILED'
233    ccol      = fclrs.FAIL
234    cf_report = 'SBCBLK.fail'
235   
236
237print(ccol+'\n\n ############ FINAL REPORT ############\n'+fclrs.ENDC)
238
239f = open(cf_report, 'w')
240
241f.write("### Sanity-check report for SBCBLK generated via 'STATION_ASF/EXP00/sbcblk_sanity_check.sh'\n\n")
242
243for 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
261cbla = '    ####################################\n'  +\
262       '    ###   TEST '+ctxt+' FOR SBCBLK !   ###\n'+\
263       '    ####################################\n\n'
264
265print(ccol+cbla+fclrs.ENDC ) ; f.write(cbla)
266
267f.close()
Note: See TracBrowser for help on using the repository browser.