Changeset 11996


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

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

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk
Files:
2 added
3 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90

    r11963 r11996  
    505505      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt 
    506506 
     507       
     508      !! Time to call the user-selected bulk parameterization for 
     509      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more... 
     510      SELECT CASE( nblk )         
     511          
     512      CASE( np_NCAR      ) 
     513         CALL turb_ncar     (     rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                         & 
     514            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     515          
     516      CASE( np_COARE_3p0 ) 
     517         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     518            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
     519            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
     520          
     521      CASE( np_COARE_3p6 ) 
     522         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     523            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
     524            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
     525          
     526      CASE( np_ECMWF     ) 
     527         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     528            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
     529            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
     530 
     531      CASE DEFAULT 
     532         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     533      END SELECT 
    507534 
    508535      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    509  
    510          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    511  
    512          CASE( np_COARE_3p0 ) 
    513             CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    514                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
    515                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    516  
    517          CASE( np_COARE_3p6 ) 
    518             CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    519                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
    520                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    521  
    522          CASE( np_ECMWF     ) 
    523             CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
    524                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     & 
    525                &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) ) 
    526  
    527          CASE DEFAULT 
    528             CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' ) 
    529          END SELECT 
    530  
    531536         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq: 
    532537         WHERE ( fr_i < 0.001_wp ) 
     
    536541         ELSEWHERE 
    537542            ! we forget about the update... 
    538             zst(:,:) = zqsb(:,:) !LB: using what we backed up before skin-algo 
    539             zsq(:,:) = zqla(:,:) !LB:  "   "   " 
     543            zst(:,:) = zqsb(:,:) !#LB: using what we backed up before skin-algo 
     544            zsq(:,:) = zqla(:,:) !#LB:  "   "   " 
    540545         END WHERE 
    541  
    542          !LB: Update of tsk, the "official" array for skin temperature 
    543          tsk(:,:) = zst(:,:) 
    544  
    545       ELSE !IF( ln_skin_cs .OR. ln_skin_wl ) 
    546  
    547          SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
    548             ! 
    549          CASE( np_NCAR      ) 
    550             CALL turb_ncar    (      rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                         & 
    551                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    552  
    553          CASE( np_COARE_3p0 ) 
    554             CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    555                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    556  
    557          CASE( np_COARE_3p6 ) 
    558             CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
    559                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    560  
    561          CASE( np_ECMWF     ) 
    562             CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
    563                &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    564  
    565          CASE DEFAULT 
    566             CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    567          END SELECT 
    568  
    569       ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl ) 
     546         tsk(:,:) = zst(:,:) !#LB: Update of tsk, the "official" array for skin temperature 
     547      END IF 
     548       
     549      !ELSE !IF( ln_skin_cs .OR. ln_skin_wl ) 
     550      ! 
     551      !   SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     552      ! 
     553      !CASE( np_COARE_3p0 ) 
     554      !      CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     555      !         &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     556      ! 
     557      !   CASE( np_COARE_3p6 ) 
     558      !CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, & 
     559      !         &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     560      ! 
     561      !   CASE( np_ECMWF     ) 
     562      !      CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  & 
     563      !         &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
     564      !      !  
     565      !   CASE DEFAULT 
     566      !      CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
     567      !   END SELECT! 
     568      ! 
     569      !ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl ) 
    570570 
    571571      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/tests/STATION_ASF/EXPREF/launch_sasf.sh

    r11930 r11996  
    4040rsync -avP ${FORC_DIR}/Station_PAPA_50N-145W*.nc ${WORK_DIR}/ 
    4141 
    42 for CASE in "ECMWF" "COARE3p6" "NCAR"; do 
     42for CASE in "ECMWF-noskin" "COARE3p6-noskin" "ECMWF" "COARE3p6" "NCAR"; do 
    4343 
    4444    echo ; echo 
  • 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.