New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
plot_station_asf_ICE.py in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/tests/STATION_ASF/EXPREF – NEMO

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

Last change on this file since 13981 was 13981, checked in by laurent, 3 years ago

Canceled changes in dynatf.F90, NEW "iom_put-able" "[uv]tau_oce" and "[uv]tau_ice" at T-points!

  • Property svn:executable set to *
File size: 11.4 KB
Line 
1#!/usr/bin/env python3
2# -*- Mode: Python; coding: utf-8; indent-tabs-mode: nil; tab-width: 4 -*-
3#
4##########################################################################################
5# Post-diagnostic of STATION_ASF for air-ice fluxes (over sea-ice)
6#
7#  L. Brodeau, 2020
8##########################################################################################
9
10import sys
11from os import path, listdir
12import argparse as ap
13from math import floor, ceil, copysign, log
14import numpy as nmp
15from netCDF4 import Dataset,num2date
16import matplotlib as mpl
17mpl.use('Agg')
18import matplotlib.pyplot as plt
19import matplotlib.dates as mdates
20
21dir_figs='.'
22size_fig=(13,8.5)
23size_fig0=(12,10)
24
25clr_red = '#AD0000'
26clr_mod = '#008ab8'
27
28rDPI=100.
29
30#ffed00: yellow ON
31#E8A727: ornage
32
33l_color = [ '0.3' , '#E8A727', '0.1'  , '#008ab8'  ] ; # colors to differentiate algos on the plot
34l_width = [   2   ,    2     ,  1.5   ,    2       ] ; # line-width to differentiate algos on the plot
35l_style = [  '-'  ,    '-'   , '--'   ,   '-'      ] ; # line-style
36
37
38# Variables to compare between algorithms
39############################################
40crealm = 'sea-ice'
41L_VNEM = [   'Cd_ice',   'Ch_ice',   'qla_ice' ,   'qsb_ice'   ,  'qt_ice'    ,  'qlw_ice'  ,  'qsr_ice'  , 'taum_ice'  ]
42L_VARO = [     'Cd'  ,     'Ch'  ,   'Qlat'    ,    'Qsen'     ,   'Qnet'     ,   'Qlw'     ,    'Qsw'    ,  'Tau'      ]
43L_VARL = [ r'$C_{D}$', r'$C_{H}$', r'$Q_{lat}$', r'$Q_{sens}$' , r'$Q_{net}$' , r'$Q_{lw}$' , r'$Q_{sw}$' , r'$|\tau|$' ]
44L_VUNT = [     ''    ,     ''    , r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$W/m^2$'  , r'$N/m^2$'  ]
45L_BASE = [    0.0005 ,    0.0005 ,      5.     ,     5.        ,     5.       ,    5.       ,     5.      ,    0.05     ]
46L_PREC = [     3     ,     3     ,      0     ,      0        ,      0      ,      0        ,     0       ,     2       ]
47L_ANOM = [   False   ,   False   ,   True      ,    True       ,    True      ,    True     ,    True     ,   True      ]
48L_MAXT = [  10000.   ,     10000.,  10000.   ,  10000.      ,  10000.    ,  10000.      , 10000.     ,      1.5    ]
49L_MINT = [    0.001  ,     0.001 , -10000.   , -10000.      , -10000.    , -10000.      ,-10000.     ,   -10000.   ]
50
51
52# About STATION_ASF output files to read:
53cpref = 'STATION_ASF-'          ; np = len(cpref)
54csuff = '_icemod.nc'            ; ns = len(csuff)
55cclnd = '_1h_YYYY0101_YYYY1231' ; nc = len(cclnd)
56
57
58################## ARGUMENT PARSING / USAGE ################################################################################################
59parser = ap.ArgumentParser(description='Generate pixel maps of a given scalar.')
60#
61requiredNamed = parser.add_argument_group('required arguments')
62requiredNamed.add_argument('-d', '--dirout' , required=True,                 help='Path to (production) directory where STATION_ASF was run')
63requiredNamed.add_argument('-f', '--forcing', required=True, default="PAPA", help='Name of forcing (ex: PAPA, ERA5_arctic')
64#
65parser.add_argument('-C', '--conf',   default="STATION_ASF",  help='specify NEMO config (ex: STATION_ASF)')
66parser.add_argument('-s', '--ystart', default="2018",         help='specify first year of experiment (ex: 2018)')
67parser.add_argument('-e', '--yend',   default="2018",         help='specify last  year of experiment (ex: 2018)')
68#
69parser.add_argument('-t', '--itype',   default="png",         help='specify the type of image you want to create (ex: png, svg, etc.)')
70#parser.add_argument('-l', '--lev' , type=int, default=0,    help='specify the level to use if 3D field (default: 0 => 2D)')
71#parser.add_argument('-I', '--ice' , action='store_true',    help='draw sea-ice concentration layer onto the field')
72#
73args = parser.parse_args()
74#
75cdir_data = args.dirout
76cforcing  = args.forcing
77#
78CONF      = args.conf
79cy1       = args.ystart
80cy2       = args.yend
81#
82fig_ext   = args.itype
83#jk    = args.lev
84#lshow_ice = args.ice
85#
86#print(''); print(' *** cdir_data = ', cdir_data); print(' *** cforcing  = ', cforcing)
87#print(' *** CONF = ', CONF); print(' *** cy1 = ', cy1); print(' *** cy2 = ', cy2)
88###############################################################################################################################################
89
90
91
92# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
93# Populating and checking existence of files to be read
94# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
95
96dir_out = cdir_data+'/output'
97ldir = listdir(dir_out)
98cf_in    = []
99list_exp = []
100list_ext = [] ; # for title
101list_frc = []
102for fn in ldir:
103    fpn = dir_out+'/'+fn
104    if path.isfile(fpn):
105        if fn[:np]==cpref and fn[-ns:]==csuff and cforcing in fn:
106            print('\n file: '+fn)
107            clab = fn[np:-nc-ns]
108            [ cexp, cfrc ] = str.split(clab, '_', 1)
109            cexp = cexp.split('-')[1] ; # removing air-sea algo name...
110            cext = cexp
111            if   cexp == "LG15":
112                cext ="Lüpkes & Gryanik, 2015"
113            elif cexp == "LU12":
114                cext ="Lüpkes et al., 2012"
115            elif cexp == "AN05":
116                cext ="Andreas et al., 2005"
117            elif cexp == "CSTC":
118                cext ="Constant coefficients"
119            print('  ===> Experiment = '+cexp+' ('+cext+'), Forcing = '+cfrc)
120            list_exp.append(cexp)
121            list_ext.append(cext)
122            list_frc.append(cfrc)
123            cf_in.append(fpn)
124nbf = len( set(list_frc) )
125
126
127
128if nbf == 0:
129    print('PROBLEM: we found no files corresponding to a forcing!')
130    sys.exit(0)
131print("list_frc =>",list_frc)
132if not nbf == 1:
133    print('PROBLEM: we found files for more that one forcing: ', set(list_frc))
134    sys.exit(0)
135
136nb_exp = len(list_exp)
137
138
139print('\n\nThere are '+str(nb_exp)+' experiments to compare:')
140for ja in range(nb_exp): print('  * '+list_exp[ja]+'\n'+'     ==> '+cf_in[ja]+'\n')
141
142if nb_exp > len(l_color):
143    print('PROBLEM: the max number of experiments for comparison is '+str(len(l_color))+' for now...')
144    sys.exit(0)
145
146
147#-----------------------------------------------------------------
148
149
150def round_bounds( x1, x2,  base=5, prec=3 ):
151    rmin =  base * round( floor(float(x1)/base), prec )
152    rmax =  base * round(  ceil(float(x2)/base), prec )
153    return rmin, rmax
154
155
156# Getting time array from the first file:
157id_in = Dataset(cf_in[0])
158vt = id_in.variables['time_counter'][:]
159cunit_t = id_in.variables['time_counter'].units ; print(' "time_counter" is in "'+cunit_t+'"')
160id_in.close()
161Nt = len(vt)
162
163vtime = num2date(vt, units=cunit_t) ; # something human!
164vtime = vtime.astype(dtype='datetime64[D]')
165
166ii=Nt/300
167ib=max(ii-ii%10,1)
168xticks_d=int(30*ib)
169
170rat = 100./float(rDPI)
171params = { 'font.family':'Open Sans',
172           'font.size':       int(15.*rat),
173           'legend.fontsize': int(15.*rat),
174           'xtick.labelsize': int(15.*rat),
175           'ytick.labelsize': int(15.*rat),
176           'axes.labelsize':  int(16.*rat)
177}
178mpl.rcParams.update(params)
179font_inf = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':18.*rat }
180font_x   = { 'fontname':'Open Sans', 'fontweight':'normal', 'fontsize':15.*rat }
181
182
183# Now we compare output variables from bulk algorithms between them:
184
185nb_var = len(L_VNEM)
186
187xF  = nmp.zeros((Nt,nb_exp))
188xFa = nmp.zeros((Nt,nb_exp))
189
190
191for jv in range(nb_var):
192    print('\n *** Treating variable: '+L_VARO[jv]+' !')
193
194    for ja in range(nb_exp):
195        #
196        id_in = Dataset(cf_in[ja])
197        xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1]  ; # only the center point of the 3x3 spatial domain!
198        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name
199        id_in.close()
200        #
201        id_toolarge, = nmp.where( xF[:,ja] > L_MAXT[jv] ) #
202        xF[id_toolarge,ja] = L_MAXT[jv]
203        id_toosmall, = nmp.where( xF[:,ja] < L_MINT[jv] ) ; #print("id_toosmall =", id_toosmall)
204        xF[id_toosmall,ja] = L_MINT[jv]
205
206    idx_okay = nmp.where( nmp.abs(xF) < 1.e+10 )
207
208    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
209    fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k')
210    ax1 = plt.axes([0.083, 0.23, 0.9, 0.7])
211    ax1.set_xticks(vtime[::xticks_d])
212    ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
213    plt.xticks(rotation='60', **font_x)
214
215    for ja in range(nb_exp):
216        fplot = nmp.ma.masked_where( xF[:,ja]==0., xF[:,ja] )
217        plt.plot(vtime, fplot, '-', color=l_color[ja], \
218                 linestyle=l_style[ja], linewidth=l_width[ja], label=list_ext[ja], alpha=0.6 ) #zorder=10+ja)
219
220    fmin, fmax = round_bounds( nmp.min(xF[idx_okay]) , nmp.max(xF[idx_okay]), base=L_BASE[jv], prec=L_PREC[jv])
221    ax1.set_ylim(fmin, fmax) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
222    plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
223
224    ax1.grid(color='k', linestyle='-', linewidth=0.3)
225    plt.legend(loc='best', ncol=1, shadow=True, fancybox=True)
226    ax1.annotate(cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
227                 ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
228                 zorder=50, **font_inf)
229    plt.savefig(L_VARO[jv]+'_'+cforcing+'_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
230    plt.close(jv)
231    #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
232
233
234
235
236    def symetric_range( pmin, pmax ):
237        # Returns a symetric f-range that makes sense for the anomaly of "f" we're looking at...
238        from math import floor, copysign, log, ceil
239        zmax = max( abs(pmax) , abs(pmin) )
240        romagn = floor(log(zmax, 10)) ; # order of magnitude of the anomaly  we're dealing with
241        rmlt = 10.**(int(romagn)) / 2.
242        frng = copysign( ceil(abs(zmax)/rmlt)*rmlt , zmax)
243        return frng
244
245
246    if L_ANOM[jv]:
247
248        for ja in range(nb_exp): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1)
249
250        if nmp.sum(nmp.abs(xFa[:,:])) == 0.0:
251            print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!')
252            print('          ==> skipping anomaly plot...')
253
254        else:
255
256            yrng = symetric_range( nmp.min(xFa) , nmp.max(xFa) )
257
258            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259            fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k')
260            ax1 = plt.axes([0.09, 0.23, 0.9, 0.7])
261            ax1.set_xticks(vtime[::xticks_d])
262            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S'))
263            plt.xticks(rotation='60', **font_x)
264
265            for ja in range(nb_exp):
266                fplot = nmp.ma.masked_where( xFa[:,ja]==0., xFa[:,ja] )
267                plt.plot(vtime, fplot, '-', color=l_color[ja], \
268                         linewidth=l_width[ja], label=list_exp[ja], alpha=0.6) #, zorder=10+ja)
269
270            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[Nt-1])
271            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']')
272            ax1.grid(color='k', linestyle='-', linewidth=0.3)
273            plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True)
274            ax1.annotate('Anomaly of '+cvar_lnm+', station: '+cforcing, xy=(0.5, 1.04), xycoords='axes fraction', \
275                         ha='center', bbox={'facecolor':'w', 'alpha':1., 'pad':10}, \
276                         zorder=50, **font_inf)
277            plt.savefig(L_VARO[jv]+'_'+cforcing+'_anomaly_'+crealm+'.'+fig_ext, dpi=int(rDPI), transparent=False)
278            plt.close(10+jv)
279            #~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Note: See TracBrowser for help on using the repository browser.