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

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

Setting up the new STATION_ASF with sea-ice support

  • Property svn:executable set to *
File size: 8.5 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,num2date
20import matplotlib as mpl
21mpl.use('Agg')
22import matplotlib.pyplot as plt
23import matplotlib.dates as mdates
24
25l_t_shift = False  ; # because time interp. is set to FALSE into "&namsbc_blk" of NEMO...
26#                    #  ==> so time array is shifted by 30 minutes but fluxes are the same (persitence of input fields ???)
27
28l_alg = [ 'ECMWF' , 'NCAR' , 'COARE3p0', 'COARE3p6', 'ANDREAS' ]
29nb_alg     =    len(l_alg)
30
31
32# Variables to check:
33l_var_rf = [  'Qsen'   ,   'Qlat'   ,  'Qlw'   ,   'Tau',   'Cd'  ,   'Ce'   ] ; # In forcing file "IDEALIZED/input_output_VALIDATION_IDEALIZED.nc"
34l_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)
35nb_var  = len(l_var_rf)
36
37dir_figs='.'
38size_fig=(13,8)
39fig_ext='png'
40
41
42rDPI=100.
43
44
45class fclrs:
46    OKGR = '\033[92m'
47    FAIL = '\033[91m'
48    ENDC = '\033[0m'
49
50
51
52# Getting arguments:
53narg = len(sys.argv)
54if not narg in [3,4]:
55    print('Usage: '+sys.argv[0]+' <forcing_+_validation_file> <NEMO-STATION_ASF_output_directory> (<m> for more/debug)'); sys.exit(0)
56cf_rf = sys.argv[1]
57cdir_out = sys.argv[2]
58l_more = False
59if narg==4:
60    l_more = ( sys.argv[3] in ['m','M'] )
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.