Ignore:
Timestamp:
2019-11-28T12:28:31+01:00 (11 months ago)
Author:
laurent
Message:

Only 1 call to each algo in "sbcblk.F90"; + added "no-skin" cases in STATION_ASF

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/plot_station_asf.py

    r11962 r11996  
    55 
    66import sys 
    7 #from os import path as path 
     7from os import path as path 
    88#from string import replace 
    99import math 
     
    2626sys.setdefaultencoding('utf8') 
    2727 
    28 cy1     = '2016' ; # First year  
     28cy1     = '2016' ; # First year 
    2929cy2     = '2018' ; # Last year 
    3030 
     
    4040 
    4141L_ALGOS = [ 'COARE3p6' , 'ECMWF'   , 'NCAR' ] 
     42l_xtrns = [ '-noskin'  , '-noskin' ,  ''    ] ; # string to add to algo name (L_ALGOS) to get version without skin params turned on 
    4243l_color = [  '#ffed00' , '#008ab8' , '0.4'  ] ; # colors to differentiate algos on the plot 
    4344l_width = [     3      ,    2      ,  1     ] ; # line-width to differentiate algos on the plot 
    4445l_style = [    '-'     ,   '-'     , '--'   ] ; # line-style 
    4546 
    46 L_VNEM  = [ 'qla'  , 'qsb'   , 'qt'   , 'qlw' , 'taum'  , 'dt_skin'  ] 
    47 L_VARO  = [  'Qlat'    ,   'Qsen'    ,  'Qnet'    ,  'Qlw'    , 'Tau'   , 'dT_skin'  ] ; # name of variable on figure 
    48 L_VARL  = [r'$Q_{lat}$',r'$Q_{sens}$',r'$Q_{net}$',r'$Q_{lw}$',r'$|\tau|$',r'$\Delta T_{skin}$' ] ; # name of variable in latex mode 
    49 L_VUNT  = [  r'$W/m^2$'   ,  r'$W/m^2$'    ,  r'$W/m^2$'   ,  r'$W/m^2$'  , r'$N/m^2$' ,   'K'     ] 
    50 L_VMAX  = [    75.     ,    75.      ,   800.     ,    25.    ,   1.2   ,   -0.7     ] 
    51 L_VMIN  = [  -250.     ,  -125.      ,  -400.     ,  -150.    ,   0.    ,    0.7     ] 
    52 L_ANOM  = [  True      ,   True      ,   True     ,   True    ,  True   ,  False     ] 
     47L_VNEM  = [   'qla'     ,     'qsb'     ,     'qt'     ,   'qlw'     ,  'taum'     ,    'dt_skin'         ] 
     48L_VARO  = [   'Qlat'    ,    'Qsen'     ,     'Qnet'   ,   'Qlw'     ,  'Tau'      ,    'dT_skin'         ] ; # name of variable on figure 
     49L_VARL  = [ r'$Q_{lat}$', r'$Q_{sens}$' , r'$Q_{net}$' , r'$Q_{lw}$' , r'$|\tau|$' , r'$\Delta T_{skin}$' ] ; # name of variable in latex mode 
     50L_VUNT  = [ r'$W/m^2$'  , r'$W/m^2$'    , r'$W/m^2$'   , r'$W/m^2$'  , r'$N/m^2$'  ,      'K'             ] 
     51L_VMAX  = [     75.     ,     75.       ,    800.      ,     25.     ,    1.2      ,      -0.7            ] 
     52L_VMIN  = [   -250.     ,   -125.       ,   -400.      ,   -150.     ,    0.       ,       0.7            ] 
     53L_ANOM  = [   True      ,    True       ,    True      ,    True     ,   True      ,      False           ] 
     54 
     55#L_VNEM  = [   'qlw'      ] 
     56#L_VARO  = [   'Qlw'      ] ; # name of variable on figure 
     57#L_VARL  = [ r'$Q_{lw}$'  ] ; # name of variable in latex mode 
     58#L_VUNT  = [ r'$W/m^2$'   ] 
     59#L_VMAX  = [     25.      ] 
     60#L_VMIN  = [   -150.      ] 
     61#L_ANOM  = [    True      ] 
     62 
     63 
    5364 
    5465nb_algos = len(L_ALGOS) ; print(nb_algos) 
     
    5667# Getting arguments: 
    5768narg = len(sys.argv) 
    58 if narg != 3: 
    59     print 'Usage: '+sys.argv[0]+' <DIR_OUT_SASF> <year>'; sys.exit(0) 
     69if narg != 2: 
     70    print 'Usage: '+sys.argv[0]+' <DIR_OUT_SASF>'; sys.exit(0) 
    6071cdir_data = sys.argv[1] 
    61 cyyyy     = sys.argv[2] 
    62  
    63  
    64  
    65 # Files to be read: 
    66 cf_in = [] 
    67 for ja in range(nb_algos): cf_in.append(cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc') 
     72 
     73 
     74 
     75# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     76# Populating and checking existence of files to be read 
     77# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     78def chck4f(cf): 
     79    cmesg = 'ERROR: File '+cf+' does not exist !!!' 
     80    if not path.exists(cf): print cmesg ; sys.exit(0) 
     81 
     82###cf_in = nmp.empty((), dtype="S10") 
     83cf_in = [] ; cf_in_ns = [] 
     84for ja in range(nb_algos): 
     85    cfi = cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc' 
     86    chck4f(cfi) 
     87    cf_in.append(cfi) 
     88# Same but without skin params: 
     89for ja in range(nb_algos): 
     90    cfi = cdir_data+'/output/'+'STATION_ASF-'+L_ALGOS[ja]+l_xtrns[ja]+'_1h_'+cy1+'0101_'+cy2+'1231_gridT.nc' 
     91    chck4f(cfi) 
     92    cf_in_ns.append(cfi) 
    6893print('Files we are goin to use:') 
    6994for ja in range(nb_algos): print(cf_in[ja]) 
     95print('   --- same without cool-skin/warm-layer:') 
     96for ja in range(nb_algos): print(cf_in_ns[ja]) 
     97#----------------------------------------------------------------- 
    7098 
    7199 
     
    97125 
    98126 
    99  
    100  
    101  
    102 for jv in range(nb_var): 
    103     print('\n *** Treating variable: '+L_VARO[jv]+' !') 
    104  
    105     for ja in range(nb_algos): 
    106         id_in = Dataset(cf_in[ja]) 
    107         xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain! 
    108         if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name 
    109         id_in.close()     
    110         #print('    => '+L_ALGOS[ja]+'  => ok!') 
    111  
    112  
    113     fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k') 
    114  
    115     ax1 = plt.axes([0.07, 0.22, 0.9, 0.75]) 
    116  
    117     ax1.set_xticks(vtime[::xticks_d]) 
    118     ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S')) 
    119     plt.xticks(rotation='60') 
    120  
    121     for ja in range(nb_algos): 
    122         plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja) 
    123  
    124     ax1.set_ylim(L_VMIN[jv], L_VMAX[jv]) ; ax1.set_xlim(vtime[0],vtime[nbr-1]) 
    125     plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']') 
    126      
    127     ax1.grid(color='k', linestyle='-', linewidth=0.3) 
    128     plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True) 
    129     ax1.annotate(cvar_lnm, xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf) 
    130     plt.savefig(L_VARO[jv]+'.'+fig_ext, dpi=int(rDPI), transparent=False) 
    131     plt.close(jv) 
    132  
    133  
    134  
    135     if L_ANOM[jv]: 
    136  
    137         for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1) 
    138  
    139         # Want a symetric y-range that makes sense for the anomaly we're looking at: 
    140         rmax = nmp.max(xFa) ; rmin = nmp.min(xFa) 
    141         rmax = max( abs(rmax) , abs(rmin) ) 
    142         romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with 
    143         rmlt = 10.**(int(romagn)) / 2. 
    144         yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax) 
    145         #print 'yrng = ', yrng ;  #sys.exit(0) 
    146          
    147         fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k') 
     127for ctest in ['skin','noskin']: 
     128 
     129    for jv in range(nb_var): 
     130        print('\n *** Treating variable: '+L_VARO[jv]+' ('+ctest+') !') 
     131 
     132        for ja in range(nb_algos): 
     133            # 
     134            if ctest == 'skin':   id_in = Dataset(cf_in[ja]) 
     135            if ctest == 'noskin': id_in = Dataset(cf_in_ns[ja]) 
     136            xF[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain! 
     137            if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name 
     138            id_in.close() 
     139 
     140        fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k') 
     141 
    148142        ax1 = plt.axes([0.07, 0.22, 0.9, 0.75]) 
    149          
     143 
    150144        ax1.set_xticks(vtime[::xticks_d]) 
    151145        ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S')) 
    152146        plt.xticks(rotation='60') 
    153          
     147 
    154148        for ja in range(nb_algos): 
    155             plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja) 
    156  
    157         ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1]) 
     149            plt.plot(vtime, xF[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja) 
     150 
     151        ax1.set_ylim(L_VMIN[jv], L_VMAX[jv]) ; ax1.set_xlim(vtime[0],vtime[nbr-1]) 
    158152        plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']') 
     153 
    159154        ax1.grid(color='k', linestyle='-', linewidth=0.3) 
    160155        plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True) 
    161         ax1.annotate('Anomaly of '+cvar_lnm, xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf) 
    162         plt.savefig(L_VARO[jv]+'_anomaly.'+fig_ext, dpi=int(rDPI), transparent=False) 
    163         plt.close(10+jv) 
    164  
     156        ax1.annotate(cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf) 
     157        plt.savefig(L_VARO[jv]+'_'+ctest+'.'+fig_ext, dpi=int(rDPI), transparent=False) 
     158        plt.close(jv) 
     159 
     160 
     161 
     162        if L_ANOM[jv]: 
     163 
     164            for ja in range(nb_algos): xFa[:,ja] = xF[:,ja] - nmp.mean(xF,axis=1) 
     165 
     166            if nmp.sum(xFa[:,:]) == 0.0: 
     167                print('     Well! Seems that for variable '+L_VARO[jv]+', choice of algo has no impact a all!') 
     168                print('          ==> skipping anomaly plot...') 
     169 
     170            else: 
     171 
     172                # Want a symetric y-range that makes sense for the anomaly we're looking at: 
     173                rmax = nmp.max(xFa) ; rmin = nmp.min(xFa) 
     174                rmax = max( abs(rmax) , abs(rmin) ) 
     175                romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with 
     176                rmlt = 10.**(int(romagn)) / 2. 
     177                yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax) 
     178                #print 'yrng = ', yrng ;  #sys.exit(0) 
     179 
     180                fig = plt.figure(num = 10+jv, figsize=size_fig, facecolor='w', edgecolor='k') 
     181                ax1 = plt.axes([0.07, 0.22, 0.9, 0.75]) 
     182 
     183                ax1.set_xticks(vtime[::xticks_d]) 
     184                ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S')) 
     185                plt.xticks(rotation='60') 
     186 
     187                for ja in range(nb_algos): 
     188                    plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linewidth=l_width[ja], label=L_ALGOS[ja], zorder=10+ja) 
     189 
     190                ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1]) 
     191                plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']') 
     192                ax1.grid(color='k', linestyle='-', linewidth=0.3) 
     193                plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True) 
     194                ax1.annotate('Anomaly of '+cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf) 
     195                plt.savefig(L_VARO[jv]+'_'+ctest+'_anomaly.'+fig_ext, dpi=int(rDPI), transparent=False) 
     196                plt.close(10+jv) 
     197 
     198 
     199 
     200 
     201# Difference skin vs noskin: 
     202xFns = nmp.zeros((nbr,nb_algos)) 
     203 
     204for jv in range(nb_var-1): 
     205    print('\n *** Treating variable: '+L_VARO[jv]+' ('+ctest+') !') 
     206 
     207    for ja in range(nb_algos-1): 
     208        id_in = Dataset(cf_in[ja]) 
     209        xF[:,ja]   = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain! 
     210        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name 
     211        id_in.close() 
     212        # 
     213        id_in = Dataset(cf_in_ns[ja]) 
     214        xFns[:,ja] = id_in.variables[L_VNEM[jv]][:,1,1] # only the center point of the 3x3 spatial domain! 
     215        if ja == 0: cvar_lnm = id_in.variables[L_VNEM[jv]].long_name 
     216        id_in.close() 
     217 
     218        xFa[:,ja] = xF[:,ja] - xFns[:,ja]  ; # difference! 
     219 
     220 
     221    # Want a symetric y-range that makes sense for the anomaly we're looking at: 
     222    rmax = nmp.max(xFa) ; rmin = nmp.min(xFa) 
     223    rmax = max( abs(rmax) , abs(rmin) ) 
     224    romagn = math.floor(math.log(rmax, 10)) ; # order of magnitude of the anomaly  we're dealing with 
     225    rmlt = 10.**(int(romagn)) / 2. 
     226    yrng = math.copysign( math.ceil(abs(rmax)/rmlt)*rmlt , rmax) 
     227    print 'yrng = ', yrng ;  #sys.exit(0) 
     228 
     229 
     230 
     231 
     232    for ja in range(nb_algos-1): 
     233 
     234        calgo = L_ALGOS[ja] 
     235 
     236        if nmp.sum(xFa[:,ja]) == 0.0: 
     237            print('     Well! Seems that for variable '+L_VARO[jv]+', and algo '+calgo+', skin param has no impact') 
     238            print('          ==> skipping difference plot...') 
     239 
     240        else: 
     241 
     242            fig = plt.figure(num = jv, figsize=size_fig, facecolor='w', edgecolor='k') 
     243            ax1 = plt.axes([0.07, 0.22, 0.9, 0.75]) 
     244 
     245            ax1.set_xticks(vtime[::xticks_d]) 
     246            ax1.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M:%S')) 
     247            plt.xticks(rotation='60') 
     248 
     249            plt.plot(vtime, xFa[:,ja], '-', color=l_color[ja], linestyle=l_style[ja], linewidth=l_width[ja], label=None, zorder=10+ja) 
     250 
     251            ax1.set_ylim(-yrng,yrng) ; ax1.set_xlim(vtime[0],vtime[nbr-1]) 
     252            plt.ylabel(L_VARL[jv]+' ['+L_VUNT[jv]+']') 
     253 
     254            ax1.grid(color='k', linestyle='-', linewidth=0.3) 
     255            #plt.legend(bbox_to_anchor=(0.45, 0.2), ncol=1, shadow=True, fancybox=True) 
     256            ax1.annotate(cvar_lnm+' ('+ctest+')', xy=(0.3, 0.97), xycoords='axes fraction',  bbox={'facecolor':'w', 'alpha':1., 'pad':10}, zorder=50, **font_inf) 
     257            plt.savefig('diff_skin-noskin_'+L_VARO[jv]+'_'+calgo+'_'+ctest+'.'+fig_ext, dpi=int(rDPI), transparent=False) 
     258            plt.close(jv) 
Note: See TracChangeset for help on using the changeset viewer.