#!/usr/bin/env python #=============================================================================== # le plot des sections West se fait : # # pour une section et une seule si # la lecture des champs se fait dans les fichiers du domaine BDY PERU (option pas MAJ pour les 3 autres regions) # # pour plusieurs sections si # la lecture des champs se fait dans les fichiers extraits du domaine TROP # # => des qu'il y a une exp bdy dans la liste des exp, alors le plot des sections from TROP est impossible sauf celles inclues dans le domaine regional #----------------------- # Nomenclature #----------------------- # upA : A upwelling area # upB : B upwelling area # upBas : B upwelling area - along shore # lon # lat # depth # i # j # k # thetao: TT, Tmin, Tmax # uo : UU, Umin, Umax #----------------------- # Decoupage des zones: #----------------------- # PERU : "cold tong" et upwelling peruvien # BENGUELA : on coupe la zone en 2 a 15S.... point critique Nord Sud # SENEGAL : Meme zone pour upA et upB # CALIFORNIA : de Santa Barbara jusqu'a la pointe Sud de la Baja California # de Santa Barbara jusqu'a Portland (meme delta de latitude) # 3 degre de largeur depuis la cote #=============================================================================== #=============================================================================== # Options #=============================================================================== debug=False # ATTENTION: si True pas dimpression images et taper pylab si on veut voir les plot de debug #pylab domain='peru'; domain_MAJ='PERU'; #domain='benguela'; domain_MAJ='BENGUELA'; #domain='senegal'; domain_MAJ='SENEGAL'; #domain='california'; domain_MAJ='CALIFORNIA'; cmip5_res='075'; cmip5_2D_interp_method='cdo' # '075' '025' 'cdo_cdo' 'scipy_cdo' #cmip5_res='025'; cmip5_2D_interp_method='scipy' # '075' '025' 'cdo_cdo' 'scipy_cdo' cmip5_interp_method=cmip5_2D_interp_method+'_cdo' kdim=75 lonlatdiffseuil=5e-06 mm=0 #----------------------------------------------------------------- # flags a nettoyer #----------------------------------------------------------------- #plot_WBDY=Obsolete #### "TT1_WBDY" correspond a la section de l'EUC extraite dans le domaine peru (from PERU) #### a priori 98.5W mais peut se regler avec iWBDY (iWBDY=2 pour 98.5W dans le 075 degre) ### Option d'avant les sections trop => ne fonctionne que pour le domaine peru ### Utile pour les exp BDY section_PERU='98.50W' # valeur unique a ajuster en fonction de iWBDY #if domain != 'peru' : # plot_WBDY=True #---------------------------------- # Obsolete depuis qu'il y a les 3 domaines. Les sections TROP correspondent a la section BDY de PERU mais pas aux autres domaines... # Conserver pour elargir a Benguela et senegal ou virer tous les plot_WBDY et iWBDY... #---------------------------------- #----------------------------------------------------------------- #---------- # SECTION TROP #---------- B_section_TROP=True # mask2 (fmask2), TT2 (fT2) & UU2 (fU2) : tmask, thetao & uo read in sections files - files extracted from trop mesh_mask & cmip5 file B_section_WBDY=True # mask1 (fmask1), TT1 (fT1) & UU1 (fU1) : tmask, thetao & uo read in domain files (peru, benguela...) - domain files extracted from trop files or generated by the regional model B_section_plot=True #---------- # ISOLEV # - 3D interpolation on trop grid (thetao) # - 2D interpolation on trop horizontal grid & cmip5 native levels # - cmip5 native grid #---------- B_T_iso=True # TT1 (fT1) : thetao read in domain files (peru, benguela...) - files extracted from trop files or generated by the regional model # TT10 (fT10) : thetao read in domain files (peru, benguela...) - files extracted from cmip5 NATIVE files B_T_iso_plot=True B_plot_cmip5_iso_native_lev=True # hfdsO & sst for trop2D & cmip5 native lev grid B_plot_cmip5_iso_native_grid=True # hfdsO & sst for cmip5 native grid # Heat Flux (Ocean) : index 200 B_hfdsO=True # hfdsO1 (fH201) : flux read in domain files (peru, benguela...) - files extracted from trop files or generated by the regional model # hfdsO10 (fH210) : flux read in domain files (peru, benguela...) - files extracted from cmip5 NATIVE files # Heat Flux (Atmosphere) : index 300 B_hfdsA=False # hfdsA (fH300) B_iso= (B_T_iso | B_section_WBDY | B_hfdsO | B_plot_cmip5_iso_native_lev | B_plot_cmip5_iso_native_grid) #----------------------------------------------------------------- #------------------------------------------------------------------------------- # Config & experience #------------------------------------------------------------------------------- exp1_list_peru12=['p12n_bdy_tr12n_quik','p12n_bdy_tr12n_quik_Tp1.5degre','p12n_bdy_soda3','p12n_bdy_soda3_T-1.5degre','p12n_bdy_tr12n_quik_wnd_erai','erai70e','erai72c','erai73c'] exp1_list_trop12_nemo=['tr12_quik','tr12_erainocrt','tr12n_n02','tr12n_n03'] exp1_list_trop025_nemo=['tr025n_gm00','tr025n_gm01','tr025_quik','tr025_quikubs','tr025_quiktau','tr025_erainocrt'] exp1_list_trop075_nemo=['tr075n_newflx','tr075n_n02nogm','tr075_quik','tr075n_quik','tr075_erainocrt','tr075n_erainocrt'] exp1_list_trop025_now=['tr025_cpl12L60_sol094'] exp1_list_trop075_now=['tr075_cpl12L60_sol093','tr075_cpl12L60_sol094','tr075now_11L60m_swrad04','tr075now_12L60_sol094obcokalb3','tr075now_44L60_sol099obcokalb3zhao3nn'] exp1_list_trop_all=exp1_list_trop12_nemo + exp1_list_trop025_nemo + exp1_list_trop075_nemo + exp1_list_trop025_now + exp1_list_trop075_now exp1_list_all_glorys=['GLORYS2V3_ORCA025_R20130808_t1','GLORYS2V3_ORCA025_R20130808_t2'] #------------------------------ # list 04/10/2019 at 12:23:44 cmip5_list_thetao=['cmip5_ACCESS1-0','cmip5_ACCESS1-3','cmip5_bcc-csm1-1','cmip5_bcc-csm1-1-m','cmip5_BNU-ESM','cmip5_CanCM4','cmip5_CanESM2','cmip5_CCSM4','cmip5_CESM1-BGC','cmip5_CESM1-CAM5-1-FV2','cmip5_CESM1-CAM5','cmip5_CESM1-FASTCHEM','cmip5_CESM1-WACCM','cmip5_CMCC-CESM','cmip5_CMCC-CM','cmip5_CMCC-CMS','cmip5_CNRM-CM5','cmip5_CNRM-CM5-2','cmip5_CSIRO-Mk3-6-0','cmip5_EC-EARTH','cmip5_FGOALS-g2','cmip5_GFDL-CM2p1','cmip5_GFDL-CM3','cmip5_GFDL-ESM2G','cmip5_GFDL-ESM2M','cmip5_GISS-E2-H','cmip5_GISS-E2-H-CC','cmip5_GISS-E2-R','cmip5_GISS-E2-R-CC','cmip5_HadCM3','cmip5_HadGEM2-AO','cmip5_HadGEM2-CC','cmip5_HadGEM2-ES','cmip5_IPSL-CM5A-LR','cmip5_IPSL-CM5A-MR','cmip5_IPSL-CM5B-LR','cmip5_MIROC-ESM','cmip5_MIROC-ESM-CHEM','cmip5_MPI-ESM-LR','cmip5_MPI-ESM-MR','cmip5_MPI-ESM-P','cmip5_MRI-CGCM3','cmip5_MRI-ESM1','cmip5_NorESM1-M','cmip5_NorESM1-ME'] cmip5_list_uo=['cmip5_ACCESS1-0','cmip5_ACCESS1-3','cmip5_bcc-csm1-1','cmip5_bcc-csm1-1-m','cmip5_BNU-ESM','cmip5_CanESM2','cmip5_CCSM4','cmip5_CESM1-BGC','cmip5_CESM1-CAM5-1-FV2','cmip5_CESM1-CAM5','cmip5_CESM1-FASTCHEM','cmip5_CESM1-WACCM','cmip5_CMCC-CESM','cmip5_CMCC-CM','cmip5_CMCC-CMS','cmip5_CNRM-CM5','cmip5_CNRM-CM5-2','cmip5_CSIRO-Mk3-6-0','cmip5_EC-EARTH','cmip5_FGOALS-g2','cmip5_GFDL-CM2p1','cmip5_GFDL-CM3','cmip5_GFDL-ESM2G','cmip5_GFDL-ESM2M','cmip5_GISS-E2-H','cmip5_GISS-E2-H-CC','cmip5_GISS-E2-R','cmip5_GISS-E2-R-CC','cmip5_HadCM3','cmip5_HadGEM2-CC','cmip5_HadGEM2-ES','cmip5_IPSL-CM5A-LR','cmip5_IPSL-CM5A-MR','cmip5_IPSL-CM5B-LR','cmip5_MIROC-ESM','cmip5_MIROC-ESM-CHEM','cmip5_MPI-ESM-LR','cmip5_MPI-ESM-MR','cmip5_MPI-ESM-P','cmip5_MRI-CGCM3','cmip5_MRI-ESM1','cmip5_NorESM1-M','cmip5_NorESM1-ME'] #cmip5_list_hfdsO=['cmip5_ACCESS1-0','cmip5_ACCESS1-3','cmip5_bcc-csm1-1','cmip5_bcc-csm1-1-m','cmip5_BNU-ESM','cmip5_CCSM4','cmip5_CMCC-CESM','cmip5_CMCC-CM','cmip5_CMCC-CMS','cmip5_CNRM-CM5','cmip5_CNRM-CM5-2','cmip5_CSIRO-Mk3-6-0','cmip5_EC-EARTH','cmip5_FGOALS-g2','cmip5_GFDL-ESM2G','cmip5_GISS-E2-R','cmip5_GISS-E2-R-CC','cmip5_MIROC-ESM','cmip5_MIROC-ESM-CHEM','cmip5_MPI-ESM-LR','cmip5_MPI-ESM-MR','cmip5_MPI-ESM-P','cmip5_MRI-CGCM3','cmip5_MRI-ESM1','cmip5_NorESM1-M','cmip5_NorESM1-ME'] # Sans MIROC cmip5_list_hfdsO=['cmip5_ACCESS1-0','cmip5_ACCESS1-3','cmip5_bcc-csm1-1','cmip5_bcc-csm1-1-m','cmip5_BNU-ESM','cmip5_CCSM4','cmip5_CMCC-CESM','cmip5_CMCC-CM','cmip5_CMCC-CMS','cmip5_CNRM-CM5','cmip5_CNRM-CM5-2','cmip5_CSIRO-Mk3-6-0','cmip5_EC-EARTH','cmip5_FGOALS-g2','cmip5_GFDL-ESM2G','cmip5_GISS-E2-R','cmip5_GISS-E2-R-CC', 'cmip5_MPI-ESM-LR','cmip5_MPI-ESM-MR','cmip5_MPI-ESM-P','cmip5_MRI-CGCM3','cmip5_MRI-ESM1','cmip5_NorESM1-M','cmip5_NorESM1-ME'] #------------------------------ # OLD LIST #cmip5_list_hfdsO=['cmip5_ACCESS1-0','cmip5_ACCESS1-3','cmip5_BNU-ESM','cmip5_CCSM4','cmip5_CMCC-CESM','cmip5_CMCC-CM','cmip5_CMCC-CMS','cmip5_CNRM-CM5-2','cmip5_CSIRO-Mk3-6-0' ,'cmip5_FGOALS-g2','cmip5_GFDL-ESM2G','cmip5_GISS-E2-R','cmip5_GISS-E2-R-CC','cmip5_MPI-ESM-LR','cmip5_MPI-ESM-MR','cmip5_MPI-ESM-P', 'cmip5_NorESM1-M','cmip5_NorESM1-ME'] #------------------------------------------------------------------------------- #exp1_list=\ #exp1_list_peru12 + \ #exp1_list_trop12_nemo + \ #exp1_list_trop025_nemo + \ #exp1_list_trop075_nemo + \ #exp1_list_trop025_now + \ #exp1_list_trop075_now + \ #exp1_list_all_cmip5_thetao + \ #exp1_list_all_cmip5_hfdsO + \ #exp1_list_all_glorys + \ exp1_list=['cmip5_CMCC-CM'] #exp1_list=['cmip5_ACCESS1-0'] #exp1_list=['cmip5_GISS-E2-H'] #exp1_list=['cmip5_CCSM4'] #exp1_list=['cmip5_EC-EARTH','cmip5_MRI-CGCM3','cmip5_MRI-ESM1'] #exp1_list=['cmip5_HadGEM2-AO'] #exp1_list=['cmip5_IPSL-CM5A-LR'] #exp1_list=['cmip5_CMCC-CM','cmip5_IPSL-CM5A-LR'] #exp1_list=['cmip5_IPSL-CM5A-LR','cmip5_CMCC-CM'] #exp1_list=['cmip5_CMCC-CM','cmip5_CMCC-CM'] #exp1_list=['cmip5_CMCC-CESM'] #exp1_list=['cmip5_BNU-ESM'] #exp1_list=['cmip5_EC-EARTH'] #exp1_list=['tr075n_newflx','tr025n_gm00','cmip5_IPSL-CM5A-LR','cmip5_CSIRO-Mk3-6-0'] #exp1_list=cmip5_list_uo #exp1_list=cmip5_list_hfdsO #exp1_list=cmip5_list_thetao + exp1_list_all_glorys #exp1_list=exp1_list_trop_all + exp1_list_all_glorys #exp1_list=exp1_list_trop_all + exp1_list_all_glorys + cmip5_list_uo #------------------------------------------------------------------------------- #------------------------------------------------------------------------------- # Sections #------------------------------------------------------------------------------- if (domain == 'peru') | (domain == 'california') : ocean='PACIFIC' section_list=['98.5W','85W','110.5W','120.25W','130W','140.5W','170.5W'] # section_list=['98.5W'] else : ocean='ATLANTIC' # section_list=['40W','30.25W','20.5W','10W','0.25W'] # island plante sur 0.25W? section_list=['30.25W','20.5W','10W'] # section_list=['40W','0.25W'] uo_max=False # print dans EUC_Upwelling_characterize_output_SENEGAL.txt # la temperature du max de l'EUC pour les sections au lieu de # la temperature moyennee en fonction du seuil defini par seuil_uo_EUC # Pour sections TROP pacifique et atlantique # profondeur pour carte lon/lat thetao kEUC1=10; kEUC2=40 zmax=34 seuil_uo_EUC=20 # temperature de l'EUC moyennee sur les n% du max de uo de l'EUC section_latlim1=-30 section_latlim2=10 depthlim1=-330 depthlim2=0 depth_text_mean=-100 depth_text_max_uo=-150 depth_text_coeur_EUC=-175 depth_text_seuil_EUC=-200 k_text_mean=20 k_text_max_uo=30 k_text_coeur_EUC=35 k_text_seuil_EUC=40 #------------------------------------------------------------------------------- # iso #------------------------------------------------------------------------------- k_iso_list=[0 ,18 ,22 ,25 ,1424 ,29 ] lev_list= ['lev1 (SST)','54m','87m','120m','29-114m','181m'] knum =[0 ,1 ,2 ,3 ,4 ,5 ] #knum =[0 ,1 ,2 ,3 ,5 ] #knum =[ 3 ,5 ] #knum =[ 4 ] #knum =[0 ] upB_as_size=4 # nombre de points de la largeur du 075 de la bande along shore de l'upwelling. # Attention si < au nombre de points de la up area alors bug visible sur le plot # des points de mer apparaissent alors comme des points de terre (blancs) Umin_max=0.5 Tmin_WBDY=12; Tmax_WBDY=28 Tmin_WBDY_seuil=14; Tmax_WBDY_seuil=22 if domain == 'peru' : lonlim1=-100 lonlim2=-70 latlim1=-30 latlim2=10 Tmin_iso_upBas_list=[20,13,13,12,13,10] Tmax_iso_upBas_list=[27,20,18,16,19,14] Tmin_iso_list=[15,12,12,11,11,10] Tmax_iso_list=[29,25,23,21,24,19] HFDS_min=0 HFDS_max=140 HFDSz_min=0 HFDSz_max=140 elif domain == 'benguela' : lonlim1=7.25 #lon T [ 464 ] 7.25 lonlim2=16.25 #lon T [ 476 ] 16.25 latlim1=-25 #lat T [ 34 ] -25.3771 latlim2=-4 #lat T [ 63 ] -4.49538 Tmin_iso_upBas_list=[21,15,13,11,14,9] Tmax_iso_upBas_list=[28,22,19,16,25,13] Tmin_iso_list=[15,11,10,9,10,8] Tmax_iso_list=[28,22,19,17,26,15] HFDS_min=0 HFDS_max=140 HFDSz_min=40 HFDSz_max=120 elif domain == 'senegal' : lonlim1=-22.75 #lon T [ 424 ] -22.75 lonlim2=-13.0 #lon T [ 437 ] -13.0 latlim1=8.22164 #lat T [ 80 ] 8.22164 latlim2=17.711 #lat T [ 93 ] 17.711 Tmin_iso_upBas_list=[21,16,14,12,15,10] Tmax_iso_upBas_list=[29,22,19,16,25,13] Tmin_iso_list=[21,15,13,12,15,11] Tmax_iso_list=[29,24,20,16,25,12] HFDS_min=0 HFDS_max=140 HFDSz_min=40 HFDSz_max=120 elif domain == 'california' : lonlim1=-132.25 #lon T [ 278 ] -132.25 lonlim2=-109.75 #lon T [ 308 ] -109.75 latlim1=22.637 #lat T [ 100 ] 22.637 latlim2=45.8778 #lat T [ 138 ] 45.8778 Tmin_iso_upBas_list=[12,10,9,8,10,7] Tmax_iso_upBas_list=[16,15,13,11,14,9] Tmin_iso_list=[12,10,9,8,10,7] Tmax_iso_list=[25,23,22,20,22,18] HFDS_min=0 HFDS_max=140 HFDSz_min=40 HFDSz_max=120 #------------------------------------------------------------------------------- # Divers #------------------------------------------------------------------------------- linewidth1=0.5 over1='y' under1='m' un=1 yb_nemo=2003; ye_nemo=2007 yb_now=1989; ye_now=2009 yb_cmip5=1989; ye_cmip5=2005 yb_glorys1='1993'; ye_glorys1='2009' yb_glorys2='2003'; ye_glorys2='2007' # Image type (none, png, eps) #img='none' img='png' # necessite mpl.use('Agg') #img='eps' # necessite mpl.use('ps') #img='pdf' # necessite mpl.use('pdf') if debug: img='none' # compute middle of lon & lat to make pcolor cells centered on T point #------------------------------------------------------------------------------- def flon_Tcentered(lon9,lat9) : #------------------------------------------------------------------------------- if lon9.ndim == 2 : jdim9,idim9=lon9.shape lon9T=np.copy(lon9) lat9T=np.copy(lat9) for jj in range(1,jdim9) : for ii in range(1,idim9) : lon9T[jj,ii]=(lon9[jj,ii-1]+lon9[jj,ii])/2 lat9T[jj,ii]=(lat9[jj-1,ii]+lat9[jj,ii])/2 elif lon9.ndim == 1 : idim9=lon9.shape[0] jdim9=lat9.shape[0] lon9T=np.copy(lon9) lat9T=np.copy(lat9) for ii in range(1,idim9) : lon9T[ii]=(lon9[ii-1]+lon9[ii])/2 for jj in range(1,jdim9) : lat9T[jj]=(lat9[jj-1]+lat9[jj])/2 # return {'lonT' : lon9T, 'latT' : lat9T} return (lon9T,lat9T) #------------------------------------------------------------------------------- def f_plot_iso(ifig99,cm99,lon99,lat99,var99,var99_str,unit99,min99,max99,min_upBas99,max_upBas99,pond99_upA,pond99_upB,pond99_upBas,exp99,img_name99,depth_iso_str99) : #------------------------------------------------------------------------------- #----------------- # Plot full domain #----------------- lon99T,lat99T=flon_Tcentered(lon99,lat99) plt.figure(ifig99) plt.pcolor(lon99T,lat99T,var99,cmap=cm99) plt.clim(min99,max99) diff99=max99-min99 if unit99 == 'C' : # temperature if diff99 < 7 : deltav99=np.linspace(min99,max99,(diff99)*2+1, endpoint=True) else : deltav99=np.linspace(min99,max99,(diff99)+1, endpoint=True) if var99_str == 'Flux Net Down (hfdsO)' : deltav99=np.linspace(min99,max99,(diff99)/20+1, endpoint=True) print "full domain : min=%.1f, max=%.1f, ticks:"%(min99,max99),deltav99 plt.colorbar(ticks=deltav99,extend='both') plt.xlim(lonlim1,lonlim2) plt.xlabel('Longitude') plt.ylim(latlim1,latlim2) plt.ylabel('Latitude') plt.grid() currentAxis = plt.gca() # trace upB currentAxis.add_patch(Rectangle((lonupB1, latupB1), lonupB2-lonupB1, latupB2-latupB1,alpha=1,edgecolor='blue',facecolor='none',linestyle='-',lw=1,hatch='/')) plt.text(lonupB1,latupB1,'Zone upB', verticalalignment='bottom',horizontalalignment='left',color='blue') # trace upA currentAxis.add_patch(Rectangle((lonupA1, latupA1), lonupA2-lonupA1, latupA2-latupA1,alpha=1,edgecolor='green',facecolor='none',linestyle='-',lw=1,hatch='\\')) plt.text(lonupA1,latupA2,'Zone upA', verticalalignment='top',horizontalalignment='left',color='green') # min,max sur toute la zone mm99="\\Large (min:%.1f max:%.1f)" %(var99.min(),var99.max()) plt.text(lonlim1,latlim1,mm99,verticalalignment='bottom',horizontalalignment='left') # latlimA=latlim1; latlimB=latlim2 vertA='bottom'; vertB='top' if (domain == 'peru') : latlimA=latlim2; latlimB=latlim1 vertA='top'; vertB='bottom' # upA pond mean plt.text(lonlim2,latlimA-2,'\\textbf{%.1f'%(pond99_upA)+' '+unit99,fontsize=48,color='green',verticalalignment=vertA,horizontalalignment='right') # upB pond mean plt.text(lonlim2,latlimB+1,'\\textbf{%.1f'%(pond99_upB)+' '+unit99,fontsize=48,color='blue',verticalalignment=vertB,horizontalalignment='right') # domain plt.text(lonlim1,latlim2-1,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # image created the... plt.text(lonlim2,latlim1,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title & save exp99_tex=exp99.replace('_','\\_') title99='\\large \\textbf{'+exp99_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{'+var99_str+'} ('+unit99+') \\textbf{'+depth_iso_str99+'} (T point)' plt.title(title99) if debug : plt.show() else : plt.savefig(img_name99) plt.close(ifig99) #--------------------- # Plot Zone upB along shore #--------------------- var99_upBas=var99[jupB1:jupB2+un:,iupB1:iupB2+un]*maskupBas img_name99_ZOOM=img_name99.replace('.'+img,'_ZOOM.'+img) plt.figure(ifig99+1) plt.pcolor(lon99T[iupB1:iupB2+un],lat99T[jupB1:jupB2+un],var99_upBas, cmap=cm99,edgecolors='blue', linewidth=0.3, linestyle='solid') plt.clim(min_upBas99,max_upBas99) diff_upBas99=max_upBas99-min_upBas99 if unit99 == 'C' : # temperature if diff_upBas99 < 7 : deltav_upBas99=np.linspace(min_upBas99,max_upBas99,(diff_upBas99)*2+1, endpoint=True) else : deltav_upBas99=np.linspace(min_upBas99,max_upBas99,(diff_upBas99)+1, endpoint=True) if var99_str == 'Flux Net Down (hfdsO)' : deltav_upBas99=np.linspace(min_upBas99,max_upBas99,(diff_upBas99)/20+1, endpoint=True) print "along shore : min=%.1f, max=%.1f, ticks:"%(min_upBas99,max_upBas99),deltav_upBas99 plt.colorbar(ticks=deltav_upBas99,extend='both') plt.xlim(lonupB1,lonupB2) plt.ylim(latupB1,latupB2) plt.grid() # domain plt.text(lonupB1,latupB2-0.2,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # min,max sur toute la zone mm99_upBas="\\Large (min:%.1f max:%.1f)" %(var99_upBas.min(),var99_upBas.max()) plt.text(lonupB1,latupB1,mm99_upBas,color='k',verticalalignment='bottom',horizontalalignment='left') # upBas pond mean plt.text(lonupB2,latupB2-0.5,'\\textbf{%.1f'%(pond99_upBas)+' '+unit99,fontsize=48,color='blue',verticalalignment='top',horizontalalignment='right') # image created the... plt.text(lonupB2,latupB1,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title & save plt.title(title99) if debug : plt.show() else : plt.savefig(img_name99_ZOOM) plt.close(ifig99+1) #=============================================================================== # Main #=============================================================================== #------------------------------------------------------------------------------- # python Modules import #------------------------------------------------------------------------------- import matplotlib as mpl print_image=True if img=='png' : mpl.use('Agg') elif img=='eps' : mpl.use('ps') elif img=='pdf' : mpl.use('pdf') else : print_image=False import matplotlib.pyplot as plt from netCDF4 import Dataset as netcdf import numpy as np import numpy.ma as ma import sys import os from matplotlib.font_manager import FontProperties import datetime import shutil import matplotlib.offsetbox as offsetbox from matplotlib.patches import Rectangle import island # fonts plt.rc('font',size=10) mpl.rc('text', usetex=True) # False pour desactiver latex. mpl.rc('font', family='serif') cm_bwr=plt.cm.bwr cm_bwr.set_over(over1) cm_bwr.set_under(under1) cm_jet=plt.cm.jet cm_jet.set_over(over1) cm_jet.set_under(under1) #------------------------------------------------------------------------------- scriptname=os.path.basename(sys.argv[0]) scriptname_root=scriptname.split('.')[0] local_time1=datetime.datetime.strftime(datetime.datetime.now(), '%Y-%m-%d') # '2018-09-26' local_time2=datetime.datetime.strftime(datetime.datetime.now(), '%d-%b-%Y') # '2018-09-26' local_time3=datetime.datetime.strftime(datetime.datetime.now(), '%Y_%m_%d-%H:%M:%S') # '2018_09_26-16:38:04' local_time4=datetime.datetime.strftime(datetime.datetime.now(), '%d-%b-%Y') # '19-Sep-2018' if uo_max : img_dir=domain_MAJ+'_EUC_uo_max' else : img_dir=domain_MAJ+'_EUC_uo_seuil_'+str(seuil_uo_EUC)+'pct' if (os.path.exists(img_dir)) : shutil.move(img_dir,img_dir+'_OLD_'+local_time3) os.makedirs(img_dir,0755) if debug : print '\nBEGIN Read fileds: ',datetime.datetime.strftime(datetime.datetime.now(), '[%H:%M:%S]\n') #------------------------------------------------------------------------------- # initialize the mean temperature & flux ascii file (out2) #------------------------------------------------------------------------------- out2=img_dir+'.txt.'+cmip5_2D_interp_method+'_'+cmip5_res if (os.path.exists(out2)) : shutil.move(out2,out2+'_OLD_'+local_time3) fout2 =open(out2, 'a') #------------ # First line #------------ if B_section_WBDY : text2a='%-7s'%('T EUC ') else : text2a='' if B_section_TROP : for ii in range(len(section_list)) : text2a=text2a+'%-7s'%('S EUC ') if B_hfdsO : text2a=text2a+'%-21s'%('HFDSO (W/m2)') if B_T_iso : for iik in range(len(knum)) : ik=knum[iik] text2a=text2a+'T %-19s'%(lev_list[ik]) text2a=text2a+'%-6s'%('id') text2a=text2a+'%-10s'%('exp name') text2a=text2a+'\n' #------------ # Second line #------------ if B_section_WBDY : text2a=text2a+'%-7s'%(section_PERU) if B_section_TROP : for ii in range(len(section_list)) : text2a=text2a+'%-7s'%(section_list[ii]) if B_hfdsO : text2a=text2a+'%-7s%-7s%-7s'%('upA','upB','upBas') # for hfdsO if B_T_iso : for ll in range(len(knum)) : text2a=text2a+'%-7s%-7s%-7s'%('upA','upB','upBas') text2a=text2a+'\n' fout2.write(text2a) #------------------------------------------------------------------------------- for exp1 in exp1_list : if ('p12' in exp1) | ('erai7' in exp1) : if B_section_TROP : print "\n\n\n exp BDY dans exp1_list ("+exp1+") => ATTENTION: TROP section < WBDY section indispensable \n\n\n" id_c5=0 id_gl=0 id_oi=0 #=============================================================================== # CHRIS EXP1 #=============================================================================== for exp1 in exp1_list : print '\n\n',80*'=','\n ',exp1,' \n',80*'=' if 'cmip5' in exp1 : exp2=exp1+'_tr'+cmip5_res+'_'+cmip5_interp_method else : exp2=exp1 exp2_tex=exp2.replace('_','\\_') #------------------------------------------------------------------------------- # Define project / config & Resolution for each exp1 #------------------------------------------------------------------------------- project='pulsation' yb1=yb_nemo; ye1=ye_nemo if ('p12' in exp1) | ('erai7' in exp1) : config='peru12_nemo';resolution='12'; elif 'tr12' in exp1 : config='trop12_nemo/'+domain;resolution='12'; elif '025_cpl' in exp1 : config='trop025_now/'+domain;resolution='025'; yb1=yb_now; ye1=ye_now elif '025' in exp1 : config='trop025_nemo/'+domain;resolution='025'; elif ('075now' in exp1) | ('075_cpl' in exp1) : config='trop075_now/'+domain; resolution='075'; yb1=yb_now; ye1=ye_now else : config='trop075_nemo/'+domain; resolution='075'; if 'cmip5' in exp1 : config='cmip5_interp/tr'+cmip5_res+'/'+cmip5_interp_method+'/'+domain; resolution=cmip5_res; yb1=yb_cmip5; ye1=ye_cmip5 project='cmip5' if 'GLORYS2V3_ORCA025_R20130808' in exp1 : project='glorys' if '_t1' in exp1 : yb1=yb_glorys1; ye1=ye_glorys1 exp2=exp1.replace('_t1','') if '_t2' in exp1 : yb1=yb_glorys2; ye1=ye_glorys2 exp2=exp1.replace('_t2','') #------------------------------------------------------------------------------- period=str(yb1)+'0101_'+str(ye1)+'1231' dir='/ccc/scratch/cont005/ra0542/hourdinc/' # irene #dir='/loceanfs/pulsation/cholod/' # ciclad dir=dir+config+'/' #------------------------------------------------------------------------------- # Area A & B Coordinates definition & Mask for domain... for 075, 025 & 12 resolution #------------------------------------------------------------------------------- # 075 coordinates : if domain == 'peru' : # iupB1=23; iupB2=33; jupB1=24; jupB2=33; # lonT[23]=-82.75; lonT[33]=-75.25; latT[24]=-13.3767 ; latT[33]=-6.7344 iupB1=21; iupB2=33; jupB1=24; jupB2=33; # lonT[23]=-84.25; lonT[33]=-75.25; latT[24]=-13.3767 ; latT[33]=-6.7344 # iupA1=13; iupA2=33; jupA1=31; jupA2=47; # lonT[13]=-90.25; lonT[33]=-75.25; latT[31]=-8.22164 ; latT[47]=3.74733 iupA1=13; iupA2=31; jupA1=31; jupA2=47; # lonT[13]=-90.25; lonT[33]=-76.75; latT[31]=-8.22164 ; latT[47]=3.74733 iWBDY=2; # iWBDY=2 correspond a 98.5W elif domain == 'benguela' : iupA1=0; iupA2=12; jupA1=0; jupA2=14; # iupB1=0; iupB2=12; jupB1=14; jupB2=29; # iWBDY=2; # iWBDY=2 correspond a 98.5W pour le PERU=> a adapter pour le BENGUELA si on veut elif domain == 'senegal' : iupB1=1; iupB2=13; jupB1=1; jupB2=12; # ATTENTION!! jupB2=13 car sinon pas de detection du premier point de mer # car dans le bas du domaine la mer touche le bord EST iupA1=1; iupA2=13; jupA1=1; jupA2=12; # iWBDY=2; # iWBDY=2 correspond a 98.5W pour le PERU=> a adapter pour le SENEGAL si on veut elif domain == 'california' : iupA1=13; iupA2=30; jupA1=0; jupA2=18; # iupB1=6; iupB2=18; jupB1=18; jupB2=37; # iWBDY=2; # iWBDY=2 correspond a 98.5W pour le PERU=> a adapter pour CALIFORNIA si on veut # => coordinates for 025 & 12 fmask1='/ccc/work/cont005/ra0542/hourdinc/trop'+resolution+'_nemo/'+domain+'/inputs/mesh_mask_TROP'+resolution+'_'+domain+'.nc' ilat1=27 # pour 075 indice correspondant a lat=30.0234 (lat Sud PERU) ilat2=83 # pour 075 indice correspondant a lat=10.4417 (lat Nord PERU) jEUC1=37; jEUC2=47 # +/-3.7473257 degres if resolution=='075' : coef=1 elif resolution=='025' : coef=3 ilat1=(ilat1 * 3)+1 ilat2=(ilat2 * 3)+1 # ilat1=82 # pour 025 indice correspondant a lat=30.0234 (lat Sud PERU) # ilat2=250 # pour 025 indice correspondant a lat=10.4417 (lat Nord PERU) elif resolution=='12' : # cas particulier des bdy peru qui necessite un mash mask avec des zeros au frontieres. if domain == 'peru' : # fmask1='/ccc/work/cont005/ra0542/hourdinc/peru12_nemo/inputs/grd/mesh_mask_p12n_bdy_soda3.nc' fmask1_OLD='/ccc/work/cont005/ra0542/hourdinc/peru12_nemo/inputs/grd/mesh_mask_peru12.L75_bathy_Stephane.nc' coef=9 ilat1=((ilat1 * 3)+1)*3 ilat2=((ilat2 * 3)+1)*3 # ilat1=246 # pour 12 indice correspondant a lat=30.0234 (lat Sud PERU) # ilat2=750 # pour 12 indice correspondant a lat=10.4417 (lat Nord PERU) iupB1=iupB1*coef iupB2=iupB2*coef jupB1=jupB1*coef jupB2=jupB2*coef iupA1=iupA1*coef iupA2=iupA2*coef jupA1=jupA1*coef jupA2=jupA2*coef iWBDY=iWBDY*coef jEUC1=jEUC1*coef jEUC2=jEUC2*coef #------------------------------------------------------------------------------- # exp1 indice definition #------------------------------------------------------------------------------- if exp1=='p12n_bdy_tr12n_quik' : id='12q1' if exp1=='p12n_bdy_tr12n_quik_Tp1.5degre' : id='12q2' if exp1=='p12n_bdy_soda3' : id='12q3' if exp1=='p12n_bdy_soda3_T-1.5degre' : id='12q4' if exp1=='p12n_bdy_tr12n_quik_wnd_erai' : id='12e1' if exp1=='erai70e' : id='12e2' fmask1=fmask1_OLD if exp1=='erai72c' : id='12q5' fmask1=fmask1_OLD if exp1=='erai73c' : id='12q6' if exp1=='tr12_quik' : id='13q1' if exp1=='tr12_erainocrt' : id='13e1' if exp1=='tr12n_n02' : id='13e2' if exp1=='tr12n_n03' : id='13e3' if exp1=='tr025n_gm00' : id='25e1' if exp1=='tr025n_gm01' : id='25e2' if exp1=='tr025_quik' : id='25q1' if exp1=='tr025_quikubs' : id='25q2' if exp1=='tr025_quiktau' : id='25q3' if exp1=='tr025_erainocrt' : id='25e3' if exp1=='tr075n_newflx' : id='75e1' if exp1=='tr075n_n02nogm' : id='75e2' if exp1=='tr075_quik' : id='75q1' if exp1=='tr075n_quik' : id='75q2' if exp1=='tr075_erainocrt' : id='75e3' if exp1=='tr075n_erainocrt' : id='75e4' if exp1=='tr025_cpl12L60_sol094' : id='25c1' if exp1=='tr075_cpl12L60_sol093' : id='75c1' if exp1=='tr075_cpl12L60_sol094' : id='75c2' if exp1=='tr075now_11L60m_swrad04' : id='75c3' if exp1=='tr075now_12L60_sol094obcokalb3' : id='75c4' if exp1=='tr075now_44L60_sol099obcokalb3zhao3nn' : id='75c5' if project == 'cmip5' : id_c5 = id_c5+1 id='c5%02d'%(id_c5) if project == 'glorys' : id_gl = id_gl+1 id='gl%02d'%(id_gl) #=============================================================================== # CHRIS ISO (mask1) #=============================================================================== if B_iso : #------------------------------------------------------------------------------- # mask1: read trop3D grid mesh_mask from domain file (peru, senegal...) #------------------------------------------------------------------------------- print 80*'-','\n on 3D trop grid : \n',80*'-' print 'domain '+domain_MAJ+' : reading mask1...' print ' fmask1 :', fmask1 ncmask1=netcdf(fmask1,'r') tmask1=ncmask1.variables['tmask'][0,:,:,:] umask1=ncmask1.variables['umask'][0,:,:,:] e1t1=ncmask1.variables['e1t'][0,:,:] e2t1=ncmask1.variables['e2t'][0,:,:] # e2u1=ncmask1.variables['e2u'][0,:,:] e3t1=ncmask1.variables['e3t'][0,:,:,:] # e3u1=ncmask1.variables['e3u'][0,:,:,:] glamt1=ncmask1.variables['glamt'][0,0,:] gphit1=ncmask1.variables['gphit'][0,:,0] glamu1=ncmask1.variables['glamu'][0,0,:] gphiu1=ncmask1.variables['gphiu'][0,:,0] deptht1=ncmask1.variables['gdept_0'][0,:] depthw1=ncmask1.variables['gdepw_0'][0,:] ncmask1.close() tmask1m=ma.masked_values(tmask1,0) umask1m=ma.masked_values(umask1,0) lonupB1=glamt1[iupB1] lonupB2=glamt1[iupB2] latupB1=gphit1[jupB1] latupB2=gphit1[jupB2] lonupA1=glamt1[iupA1] lonupA2=glamt1[iupA2] latupA1=gphit1[jupA1] latupA2=gphit1[jupA2] # print "zone upA: lon(%.2f => %.2f ) lat(%.2f => %.2f)"%(lonupA1,lonupA2,latupA1,latupA2) # print "zone upB: lon(%.2f => %.2f ) lat(%.2f => %.2f)"%(lonupB1,lonupB2,latupB1,latupB2) #------------------------------------------------------------------------------- # weight for mean computation in the domain #------------------------------------------------------------------------------- e1t1_3D=np.tile(e1t1,(kdim,1,1)) e2t1_3D=np.tile(e2t1,(kdim,1,1)) pond1a=e1t1_3D*e2t1_3D*e3t1 pond1a_sst=e1t1*e2t1 #=============================================================================== # CHRIS ISO on trop3D grid (TT1, UU1) #=============================================================================== if B_T_iso | B_section_WBDY : #------------------------------------------------------------------------------- # TT1: read thetao T3D file from domain file (peru, senegal...) #------------------------------------------------------------------------------- if project == 'pulsation' : dir1='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/outputs/'+exp1+'/' fT1=dir1+period+'/'+exp1+'_5d_'+period+'_grid_T_3D_'+domain+'_mean.nc' elif project == 'cmip5' : dir1='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/outputs/' fT1=dir1+'/'+exp1+'_1m_'+period+'_thetao_grid_T_3D_tr'+resolution+'-interp_'+cmip5_interp_method+'_'+domain+'_mean.nc' elif project == 'glorys' : dir1='/ccc/work/cont005/ra0542/hourdinc/data/'+exp2+'/'+domain+'/outputs/' fT1=dir1+'/'+exp2+'_1m_'+period+'_grid_T_3D_tr025_'+domain+'_mean.nc' ncT1=netcdf(fT1,'r') print 'domain '+domain_MAJ+' : reading thetao on 3D trop grid (TT1)...' print ' fT1 :', fT1 lon1T=ncT1.variables['nav_lon'][0,:] lat1T=ncT1.variables['nav_lat'][:,0] if abs(glamt1-lon1T).max()>lonlatdiffseuil: print '\n\n glamt1 != lon1T more than %.0e => EXIT'%(lonlatdiffseuil); print glamt1-lon1T ; sys.exit(0) if not((glamt1==lon1T).all()): print '\n CHRIS WARNING glamt1 !~= lon1T \n' if abs(gphit1-lat1T).max()>lonlatdiffseuil: print '\n\n gphit1 != lat1T more than %.0e => EXIT'%(lonlatdiffseuil); print gphit1-lat1T ; sys.exit(0) if not((gphit1==lat1T).all()): print '\n CHRIS WARNING gphit1 !~= lat1T \n' if (project == 'cmip5') & (cmip5_interp_method == 'scipy_cdo') : TT1=ncT1.variables['thetao'][:,:,:] else : if (project == 'cmip5') : TT1=ncT1.variables['thetao'][0,:,:,:] else : TT1=ncT1.variables['thetao'][0,:,:,:]*tmask1m # pond1 & mask_cmip5 75 & native_lev(iso) if project == 'cmip5' : mask_cmip5=ma.copy(TT1[:,:,:]) mask_cmip5[mask_cmip5.mask==False]=1 pond1=pond1a*mask_cmip5 else : pond1=pond1a*tmask1m if B_section_WBDY : pond1_WBDY=pond1[:,:,iWBDY] ncT1.close() if B_section_WBDY : TT1_WBDY=TT1[:,:,iWBDY] # section WBDY #------------------------------------------------------------------------------- # UU1 : read U3D file from domain file (peru, senegal...) #------------------------------------------------------------------------------- if project == 'pulsation' : dir1='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/outputs/'+exp1+'/' fU1=dir1+period+'/'+exp1+'_5d_'+period+'_grid_U_3D_'+domain+'_mean.nc' elif project == 'cmip5' : dir1='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/outputs/' fU1=dir1+'/'+exp1+'_1m_'+period+'_uo_grid_U_3D_tr'+resolution+'-interp_'+cmip5_interp_method+'_'+domain+'_mean.nc' print 'domain '+domain_MAJ+' : reading uo (UU1)...' print ' fU1 :', fU1 ncU1=netcdf(fU1,'r') lon1U=ncU1.variables['nav_lon'][0,:] if project == 'cmip5' : print '\n CHRIS WARNING :in cmip5 files, gridU=gridT => glamu1 (grid U of trop mesh_mask) != lon1U (gridT of cmip5 U file) \n',glamu1-lon1U,'\n\n' else : if abs(glamu1-lon1U).max()>lonlatdiffseuil: print '\n\n glamu1 != lon1U more than %.0e => EXIT'%(lonlatdiffseuil); print glamu1-lon1U ; sys.exit(0) lat1U=ncU1.variables['nav_lat'][:,0] if abs(gphiu1-lat1U).max()>lonlatdiffseuil: print '\n\n gphiu1 != lat1U more than %.0e => EXIT'%(lonlatdiffseuil); print gphiu1-lat1U ; sys.exit(0) if (project == 'cmip5') & (cmip5_interp_method == 'scipy_cdo') : UU1a=ncU1.variables['uo'][:,:,:] else : UU1a=ncU1.variables['uo'][0,:,:,:] UU1=UU1a*umask1m ncU1.close() UU1_WBDY=UU1[:,:,iWBDY] # section WBDY lonWBDY1_str="%.2fW"%(-glamt1[iWBDY]) # section WBDY str if lonWBDY1_str != section_PERU : print "\n\n section_PERU ("+section_PERU+") != lonWBDY1_str ("+lonWBDY1_str+") => Modifier iWBDY ou section_PERU => EXIT\n\n" ; sys.exit(0) #Calcul les poids pour chacune des profondeurs... #------------------------------------------------------------------------------- # ISOLEV lev1 (ou premier niveau sans Nan pour cmip5), 54m, 120m, 29-114m - thetao (C) #------------------------------------------------------------------------------- ii=0 pond1_mean_TT1_iso_upA_list=[] pond1_mean_TT1_iso_upB_list=[] pond1_mean_TT1_iso_upBas_list=[] # along shore cell1=pond1[:,:,:]*TT1[:,:,:] for iik in range(len(knum)) : ik=knum[iik] k_iso=k_iso_list[ik] # if (project == 'cmip5') & (k_iso == 0) : ik = 0 while not( ma.any(TT1[ik,:,:])): # while not( ma.any(TT1[0,ik,:,:])): print 'level ',ik ik += 1 k_iso=ik if k_iso == 1424 : # pond upA (cumul sur z=14:24) pond1_iso_upA=pond1[14:24,jupA1:jupA2+un:,iupA1:iupA2+un] cell1_iso_upA=cell1[14:24,jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB (cumul sur z=14:24) pond1_iso_upB=pond1[14:24,jupB1:jupB2+un:,iupB1:iupB2+un] cell1_iso_upB=cell1[14:24,jupB1:jupB2+un:,iupB1:iupB2+un] # TT1_iso (cumul sur z=14:24) pour calcul Tmoyen des niveaux 14:24 pond1_iso=pond1[14:24,:,:] cell1_iso=pond1_iso*TT1[14:24,:,:] TT1_iso=cell1_iso.sum(0)/pond1_iso.sum(0) depth_iso_str="-%.0f-%.0fm (k=14-24)"%(depthw1[14],depthw1[24]) depth_iso_str_fig="%.0fm-%.0fm_k_14-24"%(depthw1[14],depthw1[24]) else : # pond upA pond1_iso_upA=pond1[k_iso,jupA1:jupA2+un:,iupA1:iupA2+un] cell1_iso_upA=cell1[k_iso,jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB pond1_iso_upB=pond1[k_iso,jupB1:jupB2+un:,iupB1:iupB2+un] cell1_iso_upB=cell1[k_iso,jupB1:jupB2+un:,iupB1:iupB2+un] # TT1_iso TT1_iso=TT1[k_iso,:,:] depth_iso_str="-%.0fm (k=%d)"%(deptht1[k_iso],k_iso) depth_iso_str_fig="%.0fm_k_%d"%(deptht1[k_iso],k_iso) pond1_mean_TT1_iso_upA=cell1_iso_upA.sum()/pond1_iso_upA.sum() pond1_mean_TT1_iso_upB=cell1_iso_upB.sum()/pond1_iso_upB.sum() pond1_mean_TT1_iso_upA_list.append(pond1_mean_TT1_iso_upA) pond1_mean_TT1_iso_upB_list.append(pond1_mean_TT1_iso_upB) if k_iso == ik : depth_iso_str="sst (k=%d)"%(ik) # depth_iso_str_fig="sst_k_%d"%(ik) depth_iso_str_fig="sst" # ZOOM MASK + TT1 TT1_iso_upB=TT1_iso[jupB1:jupB2+un:,iupB1:iupB2+un] # Detection de la premiere maille d'ocean d'est en ouest: jdimupB,idimupB=TT1_iso_upB.shape maskupB1=np.zeros((jdimupB,idimupB)) maskupB2=np.zeros((jdimupB,idimupB)) maskupB1[TT1_iso_upB.mask]=1 B=maskupB1[:,0:idimupB-1]+maskupB1[:,1:idimupB] for jj in range(jdimupB) : iit1=np.where(B[jj,:]==1) iit=iit1[0][0] # print k_iso,jj,iit maskupB2[jj,iit-upB_as_size*coef+1:iit+1]=1 # largeur along shore en mailles pour le 075 # maskupB2[jj,iit-upB_as_size*coef+1:iit-1]=1 # largeur along shore en mailles pour le 075 # # et on evite les 2 points le long de la cote pour CMIP car depth1d mais pas 3d!!! maskupBas=ma.masked_values(maskupB2,0) # pond upBas pond1_iso_upBas=pond1_iso_upB*maskupBas cell1_iso_upBas=cell1_iso_upB*maskupBas pond1_mean_TT1_iso_upBas=cell1_iso_upBas.sum()/pond1_iso_upBas.sum() pond1_mean_TT1_iso_upBas_list.append(pond1_mean_TT1_iso_upBas) # TT1_iso TT1_iso_upBas=TT1_iso_upB*maskupBas # if k_iso == 1424 : # TT1_iso_upBas=TT1_iso_upBas*tmask1m[24,jupB1:jupB2+un:,iupB1:iupB2+un] # mean sur la Zone upB along shore # pond1_mean_TT1_iso_upBas_list.append(TT1_iso_upBas.mean()) if not(B_T_iso_plot) : print "COMPUTE pond_mean for thetao level iik=%1d K_iso=%-5d depth=%-18s upA = %.1f upB = %.1f upBas = %.1f"%(iik,k_iso,depth_iso_str,pond1_mean_TT1_iso_upA,pond1_mean_TT1_iso_upB,pond1_mean_TT1_iso_upBas) if B_T_iso_plot : print 'PLOT thetao level iik=%1d K_iso=%-5d depth=%-18s upA = %4.1f upB = %4.1f upBas = %4.1f'%(iik,k_iso,depth_iso_str,pond1_mean_TT1_iso_upA,pond1_mean_TT1_iso_upB,pond1_mean_TT1_iso_upBas) Tmin_iso=Tmin_iso_list[iik] Tmax_iso=Tmax_iso_list[iik] Tmin_iso_upB=Tmin_iso_upBas_list[iik] Tmax_iso_upB=Tmax_iso_upBas_list[iik] # print 'IK:",ik,"IIK",iik,Tmin_iso,Tmax_iso,Tmin_iso_upB,Tmax_iso_upB ifig99=17 cm99=cm_jet lon99=glamt1[:] lat99=gphit1[:] var99=TT1_iso[:,:] var99_str='thetao' unit99='C' min99=Tmin_iso_list[iik] max99=Tmax_iso_list[iik] min_upBas99=Tmin_iso_upBas_list[iik] max_upBas99=Tmax_iso_upBas_list[iik] pond99_upA=pond1_mean_TT1_iso_upA pond99_upB=pond1_mean_TT1_iso_upB pond99_upBas=pond1_mean_TT1_iso_upBas exp99=exp2 img_name99=img_dir+'/'+exp2+'_'+domain_MAJ+'_thetao_'+depth_iso_str_fig+'.'+img depth_iso_str99=depth_iso_str f_plot_iso(ifig99,cm99,lon99,lat99,var99,var99_str,unit99,min99,max99,min_upBas99,max_upBas99,pond99_upA,pond99_upB,pond99_upBas,exp99,img_name99,depth_iso_str99) #=============================================================================== # CHRIS ISO on trop2D grid (hfdsO1, sst1) # # hfdsO on the trop2D grid # sst (only for cmip5) on the trop2D grid & cmip5 native lev #=============================================================================== if B_hfdsO | B_plot_cmip5_iso_native_lev : #------------------------------------------------------------------------------- # if cmip5 : read cmip5 trop2D native lev sst => make cmip5_mask for hfdsO or plot sst # for all : compute pond1_sst to evaluate the mean for hfdsO & sst #------------------------------------------------------------------------------- if project == 'cmip5' : #------------------------------------------------------------------------------- # read cmip5 trop2D native lev sst => make cmip5_mask for hfdsO or plot sst #------------------------------------------------------------------------------- dir1='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/outputs/' fT1_sst=dir1+'/'+exp1+'_1m_'+period+'_thetao_grid_T_3D_tr'+resolution+'-interp_'+cmip5_2D_interp_method+'-lev_natif_'+domain+'_mean.nc' ncT1_sst=netcdf(fT1_sst,'r') if (cmip5_interp_method == 'scipy_cdo') : sst1=ncT1_sst.variables['thetao'][0,:,:] else : sst1=ncT1_sst.variables['thetao'][0,0,:,:] lev1_sst=ncT1_sst.variables['lev'][0] # make mask for cmip5 hfdsO mask_cmip5_sst=ma.copy(sst1[:,:]) mask_cmip5_sst[mask_cmip5_sst.mask==False]=1 pond1_sst=pond1a_sst*mask_cmip5_sst ncT1_sst.close() else : sst1=TT1[0,:,:] lev1_sst=deptht1[0] pond1_sst=pond1a_sst*tmask1m[0,:,:] # pond upA pond1_sst_upA=pond1_sst[jupA1:jupA2+un:,iupA1:iupA2+un] cell_sst1_upA=pond1_sst_upA*sst1[jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB pond1_sst_upB=pond1_sst[jupB1:jupB2+un:,iupB1:iupB2+un] cell_sst1_upB=pond1_sst_upB*sst1[jupB1:jupB2+un:,iupB1:iupB2+un] # pond1_mean_sst1_upA=cell_sst1_upA.sum()/pond1_sst_upA.sum() pond1_mean_sst1_upB=cell_sst1_upB.sum()/pond1_sst_upB.sum() # ZOOM MASK + TT1 sst1_upB=sst1[jupB1:jupB2+un:,iupB1:iupB2+un] # Detection de la premiere maille d'ocean d'est en ouest: jdimupB,idimupB=sst1_upB.shape maskupB1=np.zeros((jdimupB,idimupB)) maskupB2=np.zeros((jdimupB,idimupB)) maskupB1[sst1_upB.mask]=1 B=maskupB1[:,0:idimupB-1]+maskupB1[:,1:idimupB] for jj in range(jdimupB) : iit1=np.where(B[jj,:]==1) iit=iit1[0][0] # print k_iso,jj,iit maskupB2[jj,iit-upB_as_size*coef+1:iit+1]=1 # largeur along shore en mailles pour le 075 # maskupB2[jj,iit-upB_as_size*coef+1:iit-1]=1 # largeur along shore en mailles pour le 075 # # et on evite les 2 points le long de la cote pour CMIP car depth1d mais pas 3d!!! maskupBas=ma.masked_values(maskupB2,0) # pond upBas pond1_sst_upBas=pond1_sst_upB*maskupBas cell_sst1_upBas=cell_sst1_upB*maskupBas pond1_mean_sst1_upBas=cell_sst1_upBas.sum()/pond1_sst_upBas.sum() if B_plot_cmip5_iso_native_lev & ('cmip5' in exp1) & B_T_iso : print "PLOT sst (C) native LEV depth=%-10.1f upA = %.1f upB = %.1f upBas = %.1f"%(lev1_sst,pond1_mean_sst1_upA,pond1_mean_sst1_upB,pond1_mean_sst1_upBas) ifig99=39 cm99=cm_jet lon99=glamt1[:] lat99=gphit1[:] var99=sst1[:,:] var99_str='sst' unit99='C' min99=Tmin_iso_list[0] max99=Tmax_iso_list[0] min_upBas99=Tmin_iso_upBas_list[0] max_upBas99=Tmax_iso_upBas_list[0] pond99_upA=pond1_mean_sst1_upA pond99_upB=pond1_mean_sst1_upB pond99_upBas=pond1_mean_sst1_upBas exp99=exp1+'_tr'+cmip5_res+'_'+cmip5_2D_interp_method+'_native_lev' img_name99=img_dir+'/'+exp99+'_'+domain_MAJ+'_sst.'+img depth_iso_str99='' f_plot_iso(ifig99,cm99,lon99,lat99,var99,var99_str,unit99,min99,max99,min_upBas99,max_upBas99,pond99_upA,pond99_upB,pond99_upBas,exp99,img_name99,depth_iso_str99) else : print "COMPUTE pond_mean for sst (C) native depth=%-.1fm upA = %.1f upB = %.1f upBas = %.1f\n"%(lev1_sst,pond1_mean_sst1_upA,pond1_mean_sst1_upB,pond1_mean_sst1_upBas) if B_hfdsO : #------------------------------------------------------------------------------- # read hfdsO #------------------------------------------------------------------------------- if project == 'pulsation' : # CHRIS a implementer dir1='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/outputs/'+exp1+'/' # fT1=dir1+period+'/'+exp1+'_5d_'+period+'_grid_T_3D_'+domain+'_mean.nc' elif project == 'cmip5' : dir1='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/outputs/' fH201=dir1+'/'+exp1+'_1m_'+period+'_hfdsO_grid_T_2D_tr'+cmip5_res+'-interp_'+cmip5_2D_interp_method+'_'+domain+'_mean.nc' ncH201=netcdf(fH201,'r') if (project == 'cmip5') : if (cmip5_interp_method == 'scipy_cdo') : hfdsO1=ncH201.variables['hfds'][:,:] else : hfdsO1=ncH201.variables['hfds'][0,:,:] ncH201.close() elif project == 'glorys' : # CHRIS a implementer dir1='/ccc/work/cont005/ra0542/hourdinc/data/'+exp2+'/'+domain+'/outputs/' # fT1=dir1+'/'+exp2+'_1m_'+period+'_grid_T_3D_tr025_'+domain+'_mean.nc' print '\n',80*'-','\n on 2D trop grid & native LEV : \n',80*'-' print 'domain %s : reading hfdsO on 2D trop grid (hfdsO1)...'%(domain_MAJ) print ' fH201 :', fH201 print 'domain %s : reading thetao on 2D trop grid & native levels (sst1)...'%(domain_MAJ) print ' fT1_sst :', fT1_sst # pond upA pond1_sst_upA=pond1_sst[jupA1:jupA2+un:,iupA1:iupA2+un] cell_hfdsO1_upA=pond1_sst_upA*hfdsO1[jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB pond1_sst_upB=pond1_sst[jupB1:jupB2+un:,iupB1:iupB2+un] cell_hfdsO1_upB=pond1_sst_upB*hfdsO1[jupB1:jupB2+un:,iupB1:iupB2+un] # pond1_mean_hfdsO1_upA=cell_hfdsO1_upA.sum()/pond1_sst_upA.sum() pond1_mean_hfdsO1_upB=cell_hfdsO1_upB.sum()/pond1_sst_upB.sum() # pond upBas pond1_sst_upBas=pond1_sst_upB*maskupBas cell_hfdsO1_upBas=cell_hfdsO1_upB*maskupBas pond1_mean_hfdsO1_upBas=cell_hfdsO1_upBas.sum()/pond1_sst_upBas.sum() if B_plot_cmip5_iso_native_lev : print "\nPLOT hfdsO (W/m2) native LEV depth=%-10.1f upA = %4.1f upB = %4.1f upBas = %4.1f"%(lev1_sst,pond1_mean_hfdsO1_upA,pond1_mean_hfdsO1_upB,pond1_mean_hfdsO1_upBas) ifig99=37 cm99=cm_jet lon99=glamt1[:] lat99=gphit1[:] var99=hfdsO1[:,:] var99_str='Flux Net Down (hfdsO)' unit99='W/m2' min99=HFDS_min max99=HFDS_max min_upBas99=HFDSz_min max_upBas99=HFDSz_max pond99_upA=pond1_mean_hfdsO1_upA pond99_upB=pond1_mean_hfdsO1_upB pond99_upBas=pond1_mean_hfdsO1_upBas exp99=exp1+'_tr'+cmip5_res+'_'+cmip5_2D_interp_method img_name99=img_dir+'/'+exp99+'_'+domain_MAJ+'_hfdsO.'+img depth_iso_str99='' f_plot_iso(ifig99,cm99,lon99,lat99,var99,var99_str,unit99,min99,max99,min_upBas99,max_upBas99,pond99_upA,pond99_upB,pond99_upBas,exp99,img_name99,depth_iso_str99) else : print "\nCOMPUTE pond_mean for hfdsO (W/m2) native depth=%-.1fm upA = %.1f upB = %.1f upBas = %.1f"%(lev1_sst,pond1_mean_hfdsO1_upA,pond1_mean_hfdsO1_upB,pond1_mean_hfdsO1_upBas) #=============================================================================== # CHRIS ISO on cmip5 native grid (hfdsO1, sst1) # # only for cmip5 !!! # # hfdsO & sst #=============================================================================== if B_plot_cmip5_iso_native_grid & ('cmip5' in exp1) : if B_T_iso : #------------------------------------------------------------------------------- # plot SST NATIVE grid #------------------------------------------------------------------------------- # open & read NATIVE grid thetao file dir10='/ccc/scratch/cont005/ra0542/hourdinc/cmip5/outputs/plot_native_grid/thetao/' fT10=dir10+'/'+exp1+'_1m_'+period+'_thetao_grid_native_with_lon_2D_shifted_mean.nc' ncT10=netcdf(fT10,'r') lon10=ncT10.variables['lon'][:] lon10=lon10-360 lat10=ncT10.variables['lat'][:] depth10t=ncT10.variables['lev'][:] TT10=ncT10.variables['thetao'][:,:,:] ncT10.close() print '\n',80*'-','\n on NATIVE GRID : \n',80*'-' print 'domain '+domain_MAJ+' : reading thetao on NATIVE grid (TT10)...' print ' fT10 :', fT10 #----------------- # Plot SST full domain #----------------- print "PLOT sst (C) on native GRID depth=%-.1f\n"%(depth10t[0]) plt.figure(21) Tmin_iso=Tmin_iso_list[0] Tmax_iso=Tmax_iso_list[0] v=np.linspace(Tmin_iso,Tmax_iso,Tmax_iso-Tmin_iso+1, endpoint=True) lon10T,lat10T=flon_Tcentered(lon10,lat10) plt.pcolor(lon10T[:],lat10T[:],TT10[0,:,:],cmap=cm_jet) plt.clim(Tmin_iso,Tmax_iso) plt.colorbar(ticks=v,extend='both') plt.xlim(lonlim1,lonlim2) plt.xlabel('Longitude') plt.ylim(latlim1,latlim2) plt.ylabel('Latitude') plt.grid() currentAxis = plt.gca() # trace upB currentAxis.add_patch(Rectangle((lonupB1, latupB1), lonupB2-lonupB1, latupB2-latupB1,alpha=1,edgecolor='blue',facecolor='none',linestyle='-',lw=1,hatch='/')) plt.text(lonupB1,latupB1,'Zone upB', verticalalignment='bottom',horizontalalignment='left',color='blue') # trace upA currentAxis.add_patch(Rectangle((lonupA1, latupA1), lonupA2-lonupA1, latupA2-latupA1,alpha=1,edgecolor='green',facecolor='none',linestyle='-',lw=1,hatch='\\')) plt.text(lonupA1,latupA2,'Zone upA', verticalalignment='top',horizontalalignment='left',color='green') # domain plt.text(lonlim1+1.5,latlim2-1.5,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonlim1+1.5,latlim1+1.5,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='bottom',horizontalalignment='left') # image created the... plt.text(lonlim2,latlim1,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title ### carefull ... exp2 is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1 is used for title (cmip5_IPSL-CM5A-LR_thetao...) because we plot the NATIVE GRID exp1_tex=exp1.replace('_','\\_') # depth_verif_str10="lev1 -%.0fm"%(depth10t[0]) depth_verif_str10=" -%.0fm"%(depth10t[0]) # title21='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{thetao} (C) \\textbf{'+depth_verif_str10+'}' title21='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{sst} (C) \\textbf{'+depth_verif_str10+'} (T point)' plt.title(title21) # save img if debug : plt.show() else : plt.savefig(img_dir+'/'+exp1+'_native_grid_'+domain_MAJ+'_sst.'+img) #--------------------- # Plot SST Zone upB along shore #--------------------- plt.figure(22) Tmin_iso_upB=Tmin_iso_upBas_list[0] Tmax_iso_upB=Tmax_iso_upBas_list[0] v=np.linspace(Tmin_iso_upB,Tmax_iso_upB,Tmax_iso_upB-Tmin_iso_upB+1, endpoint=True) plt.pcolor(lon10T[:],lat10T[:],TT10[0,:,:],cmap=cm_jet,edgecolors='blue', linewidth=0.3, linestyle='solid') plt.clim(Tmin_iso_upB,Tmax_iso_upB) plt.colorbar(ticks=v,extend='both') plt.xlim(lonupB1,lonupB2) plt.xlabel('Longitude') plt.ylim(latupB1,latupB2) plt.ylabel('Latitude') plt.grid() currentAxis = plt.gca() # domain plt.text(lonupB1,latupB2,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonupB1,latupB1,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='bottom',horizontalalignment='left') # image created the... plt.text(lonupB2,latupB1,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title ### carefull ... exp2 is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1 is used for title (cmip5_IPSL-CM5A-LR_thetao...) because we plot the NATIVE GRID title22='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{sst} (C) \\textbf{'+depth_verif_str10+'} (T point)' plt.title(title22) # save img if debug : plt.show() else : plt.savefig(img_dir+'/'+exp1+'_native_grid_'+domain_MAJ+'_sst_ZOOM.'+img) if B_hfdsO : #------------------------------------------------------------------------------- # plot HFDS O NATIVE grid #------------------------------------------------------------------------------- # NATIVE grid HFDS O file dir210='/ccc/scratch/cont005/ra0542/hourdinc/cmip5/outputs/plot_native_grid/hfdsO/' fH210=dir210+'/'+exp1+'_1m_'+period+'_hfdsO_grid_native_with_lon_2D_shifted_mean.nc' ncH210=netcdf(fH210,'r') lon210=ncH210.variables['lon'][:] lon210=lon210-360 lat210=ncH210.variables['lat'][:] hfdsO10=ncH210.variables['hfds'][:,:] ncH210.close() print 'domain '+domain_MAJ+' : reading hfdsO on NATIVE grid (hfdsO10)...' print ' fH210 :', fH210 #----------------- # Plot HFDS O full domain #----------------- print "PLOT hfdsO (W/m2) on native GRID" plt.figure(23) lon210T,lat210T=flon_Tcentered(lon210,lat210) v=np.linspace(HFDS_max,HFDS_min,(HFDS_max-HFDS_min)/20+1, endpoint=True) plt.pcolor(lon210T[:],lat210T[:],hfdsO10[:,:],cmap=cm_jet) plt.clim(HFDS_min,HFDS_max) plt.colorbar(ticks=v,extend='both') plt.xlim(lonlim1,lonlim2) plt.xlabel('Longitude') plt.ylim(latlim1,latlim2) plt.ylabel('Latitude') plt.grid() currentAxis = plt.gca() # trace upB currentAxis.add_patch(Rectangle((lonupB1, latupB1), lonupB2-lonupB1, latupB2-latupB1,alpha=1,edgecolor='blue',facecolor='none',linestyle='-',lw=1,hatch='/')) plt.text(lonupB1,latupB1,'Zone upB', verticalalignment='bottom',horizontalalignment='left',color='blue') # trace upA currentAxis.add_patch(Rectangle((lonupA1, latupA1), lonupA2-lonupA1, latupA2-latupA1,alpha=1,edgecolor='green',facecolor='none',linestyle='-',lw=1,hatch='\\')) plt.text(lonupA1,latupA2,'Zone upA', verticalalignment='top',horizontalalignment='left',color='green') # domain plt.text(lonlim1+1.5,latlim2-1.5,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonlim1+1.5,latlim1+1.5,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='bottom',horizontalalignment='left') # image created the... plt.text(lonlim2,latlim1,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title ### carefull ... exp2 is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1 is used for title (cmip5_IPSL-CM5A-LR_thetao...) because we plot the NATIVE GRID exp1_tex=exp1.replace('_','\\_') title23='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{Flux Net Down (hfdsO)} (W/m2) (T point)' plt.title(title23) # save img if debug : plt.show() else : plt.savefig(img_dir+'/'+exp1+'_native_grid_'+domain_MAJ+'_hfdsO.'+img) #--------------------- # Plot HFDS O Zone upB along shore #--------------------- plt.figure(24) v=np.linspace(HFDS_max,HFDS_min,(HFDS_max-HFDS_min)/20+1, endpoint=True) plt.pcolor(lon210T[:],lat210T[:],hfdsO10[:,:],cmap=cm_jet,edgecolors='blue', linewidth=0.3, linestyle='solid') plt.clim(HFDSz_min,HFDSz_max) plt.colorbar(ticks=v,extend='both') plt.xlim(lonupB1,lonupB2) plt.xlabel('Longitude') plt.ylim(latupB1,latupB2) plt.ylabel('Latitude') plt.grid() currentAxis = plt.gca() # domain plt.text(lonupB1,latupB2,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonupB1,latupB1,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='bottom',horizontalalignment='left') # image created the... plt.text(lonupB2,latupB1,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title ### carefull ... exp2 is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1 is used for title (cmip5_IPSL-CM5A-LR_thetao...) because we plot the NATIVE GRID title24='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{Flux Net Down (hfdsO)} (W/m2) (T point)' plt.title(title24) # save img if debug : plt.show() else : plt.savefig(img_dir+'/'+exp1+'_native_grid_'+domain_MAJ+'_hfdsO_ZOOM.'+img) #=============================================================================== # CHRIS SECTION # get mask1, TT1 & UU1 from domain part of the script below # & read mask2, TT2 & UU2 from trop sections #=============================================================================== if B_section_TROP | B_section_WBDY : print 80*'-','\n on SECTIONS trop grid : \n',80*'-' #------------------------------------------------------------------------------- # Initialize for sections #------------------------------------------------------------------------------- section_from_domain1=domain_MAJ gphit4_list=[] deptht4_list=[] TTUUmax_list=[] UUU_WBDY_list=[] TTT_WBDY_list=[] pond_WBDY_list=[] lonWBDY_str_list=[] section_from_domain_list=[] pond_mean_TTT_WBDY=[] #------------------------------------------------------------------------------- # Fill the T & U sections array for BDY & TROP sections #------------------------------------------------------------------------------- if B_section_WBDY : #------------------------------------------------------------------------------- # UU1 TT1 already read in domain part #------------------------------------------------------------------------------- # Pour eviter les sections > WBDY: if config != 'peru12_nemo' : ##### if ('p12' in exp1) | ('erai7' in exp1) : gphit4_list.append(gphit1) deptht4_list.append(deptht1) UUU_WBDY_list.append(UU1_WBDY) TTT_WBDY_list.append(TT1_WBDY) pond_WBDY_list.append(pond1_WBDY) lonWBDY_str_list.append(lonWBDY1_str) section_from_domain_list.append(section_from_domain1) if B_section_TROP : for section in section_list: #------------------------------------------------------------------------------- # mask2 : read TROP section mask #------------------------------------------------------------------------------- print 'TROP Section '+section+' : reading mask (mask2) & thetao (TT2) & uo (UU2)...' fmask2='/ccc/work/cont005/ra0542/hourdinc/trop'+resolution+'_nemo/inputs/mesh_mask_TROP'+resolution+'_section_'+section+'.nc' print ' fmask2 :', fmask2 ncmask2=netcdf(fmask2,'r') tmask2=ncmask2.variables['tmask'][0,:,ilat1:ilat2+un,0] umask2=ncmask2.variables['umask'][0,:,ilat1:ilat2+un,0] e1t2=ncmask2.variables['e1t'][0,ilat1:ilat2+un,0] e2t2=ncmask2.variables['e2t'][0,ilat1:ilat2+un,0] # e2u2=ncmask2.variables['e2u'][0,ilat1:ilat2+un,0] e3t2=ncmask2.variables['e3t'][0,:,ilat1:ilat2+un,0] # e3u2=ncmask2.variables['e3u'][0,:,ilat1:ilat2+un,0] glamt2=ncmask2.variables['glamt'][0,0,:] gphit2=ncmask2.variables['gphit'][0,ilat1:ilat2+un,0] glamu2=ncmask2.variables['glamu'][0,0,:] gphiu2=ncmask2.variables['gphiu'][0,ilat1:ilat2+un,0] deptht2=ncmask2.variables['gdept_0'][0,:] depthw2=ncmask2.variables['gdepw_0'][0,:] ncmask2.close() tmask2m=ma.masked_values(tmask2,0) umask2m=ma.masked_values(umask2,0) #------------------------------------------------------------------------------- # weight for TROP section mean computation #------------------------------------------------------------------------------- e1t2_2D=np.tile(e1t2,(kdim,1)) e2t2_2D=np.tile(e2t2,(kdim,1)) pond2_WBDY=e1t2_2D*e2t2_2D*e3t2*tmask2m #------------------------------------------------------------------------------- # read T3D file from TROP (TT2) #------------------------------------------------------------------------------- if project == 'pulsation' : dir2='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/../outputs/'+exp1+'/' fT2=dir2+'section_'+section+'/'+exp1+'_5d_'+period+'_grid_T_3D_section_'+section+'_mean.nc' elif project == 'cmip5' : dir2='/ccc/scratch/cont005/ra0542/hourdinc/'+config+'/../outputs/../' # fT2=dir2+'section_'+section+'/'+exp1+'_1m_'+period+'_thetao_grid_T_3D_tr075_section_'+section+'_mean.nc' fT2=dir2+'section_'+section+'/'+exp1+'_1m_'+period+'_thetao_grid_T_3D_tr'+resolution+'-interp_'+cmip5_interp_method+'_section_'+section+'_mean.nc' elif project == 'glorys' : dir2='/ccc/work/cont005/ra0542/hourdinc/data/'+exp2+'/' fT2=dir2+'section_'+section+'/'+exp2+'_1m_'+period+'_grid_T_3D_tr025_section_'+section+'_mean.nc' print ' fT2 :', fT2 ncT2=netcdf(fT2,'r') lon2T=ncT2.variables['nav_lon'][0,:] lat2T=ncT2.variables['nav_lat'][ilat1:ilat2+un,0] if abs(glamt2-lon2T).max()>lonlatdiffseuil: print '\n\n glamt2 != lon2T more than %.0e => EXIT'%(lonlatdiffseuil); print glamt2-lon2T ; sys.exit(0) if not((glamt2==lon2T).all()): print '\n CHRIS WARNING glamt2 !~= lon2T \n' if abs(gphit2-lat2T).max()>lonlatdiffseuil: print '\n\n gphit2 != lat2T more than %.0e => EXIT'%(lonlatdiffseuil); print gphit2-lat2T ; sys.exit(0) if not((gphit2==lat2T).all()): print '\n CHRIS WARNING gphit2 !~= lat2T \n' if (project == 'cmip5') & (cmip5_interp_method == 'scipy_cdo') : TT2a=ncT2.variables['thetao'][:,ilat1:ilat2+un,0] else : TT2a=ncT2.variables['thetao'][0,:,ilat1:ilat2+un,0] # if project == 'cmip5' : # if not((exp1 == 'cmip5_EC-EARTH') | (exp1 == 'cmip5_GFDL-CM2p1')) : # TT2a=TT2a-273.15 TT2_WBDY=TT2a*tmask2m ncT2.close() #------------------------------------------------------------------------------- # read U3D file from TROP (UU2) #------------------------------------------------------------------------------- if project == 'pulsation' : fU2=dir2+'section_'+section+'/'+exp1+'_5d_'+period+'_grid_U_3D_section_'+section+'_mean.nc' elif project == 'cmip5' : fU2=dir2+'section_'+section+'/'+exp1+'_1m_'+period+'_uo_grid_U_3D_tr'+resolution+'-interp_'+cmip5_interp_method+'_section_'+section+'_mean.nc' elif project == 'glorys' : fU2=dir2+'section_'+section+'/'+exp2+'_1m_'+period+'_grid_U_3D_tr025_section_'+section+'_mean.nc' print ' fU2 :', fU2 ncU2=netcdf(fU2,'r') lon2U=ncU2.variables['nav_lon'][0,:] if project == 'cmip5' : print "\n CHRIS WARNING :in cmip5 files, gridU=gridT => glamu2=%.2f (grid U of trop mesh_mask) != lon2U=%.2f (gridT of cmip5 U file) => diff=%.3f \n\n"%(glamu2[0],lon2U[0],glamu2[0]-lon2U[0]) else : if not( 'GLORYS2V3_ORCA025' in exp1) : if abs(glamu2-lon2U).max()>lonlatdiffseuil: print "\n\n glamu2 != lon2U more than ",lonlatdiffseuil,"\n\n",glamu2-lon2U,"\n\n => EXIT\n\n" ; sys.exit(0) else : print "\n CHRIS WARNING :in GLORYS files, gridU=gridT => glamu2=%.2f (grid U of trop mesh_mask) != lon2U=%.2f (gridT of GLORYS U file) => diff=%.3f \n\n"%(glamu2[0],lon2U[0],glamu2[0]-lon2U[0]) lat2U=ncU2.variables['nav_lat'][ilat1:ilat2+un,0] if abs(gphiu2-lat2U).max()>lonlatdiffseuil: print "\n\n gphiu2 != lat2U more than ",lonlatdiffseuil,"\n\n",gphiu2-lat2U,"\n\n => EXIT\n\n" ; sys.exit(0) if (project == 'cmip5') & (cmip5_interp_method == 'scipy_cdo') : UU2a=ncU2.variables['uo'][:,ilat1:ilat2+un,0] else : UU2a=ncU2.variables['uo'][0,:,ilat1:ilat2+un,0] UU2_WBDY=UU2a*umask2m ncU2.close() lonWBDY2_str="%.2fW"%(-glamt2[0]) # section WBDY str section_from_domain2='TROP' gphit4_list.append(gphit2) deptht4_list.append(deptht2) UUU_WBDY_list.append(UU2_WBDY) TTT_WBDY_list.append(TT2_WBDY) pond_WBDY_list.append(pond2_WBDY) lonWBDY_str_list.append(lonWBDY2_str) section_from_domain_list.append(section_from_domain2) #------------------------------------------------------------------------------- # For each section (BDY & TROP sections) # find EUC center #------------------------------------------------------------------------------- for ii in range(len(UUU_WBDY_list)) : gphit4=gphit4_list[ii] deptht4=deptht4_list[ii] UUU_WBDY=UUU_WBDY_list[ii] TTT_WBDY=TTT_WBDY_list[ii] pond_WBDY=pond_WBDY_list[ii] lonWBDY_str=lonWBDY_str_list[ii] section_from_domain=section_from_domain_list[ii] if B_section_plot : print 'PLOT section '+lonWBDY_str+' from '+section_from_domain+'...' else : print 'COMPUTE section '+lonWBDY_str+' from '+section_from_domain+'...' #---------------------------------------------- # MAX de uo sur tout le domaine: #---------------------------------------------- kUUmax,jUUmax=np.argwhere(UUU_WBDY == UUU_WBDY.max())[0] TTUUmax=TTT_WBDY[kUUmax,jUUmax] TTUUmax_str='\\Large Max uo:[%.2f m/s %.1f C]'%(UUU_WBDY.max(),TTUUmax) #---------------------------------------------- # Detection du max et de la zone forte / seuil #---------------------------------------------- seuil=UUU_WBDY.max()*seuil_uo_EUC/100 seuil_str='\\Large Selection: \\textbf{%.0f'%(seuil_uo_EUC)+'\%} de uo max' vv1=np.copy(UUU_WBDY) vv1[vv1 >= seuil]=1 vv1[np.isnan(vv1)]=1 # pour mettre les NaN a zero quand interpolation faite avec scipy (pas de fill_value a 1e+20) vv1[vv1 < seuil]=0 if debug : plt.figure(3000+mm) plt1a=plt.pcolormesh(vv1[:,:],cmap=cm_jet) # plt.ylim(yf,y0) # plt.xlim(35,45) # plt.ylim(54,10) ax = plt.gca() ax.invert_yaxis() plt.grid() plt.title(exp2_tex+' '+lonWBDY_str) plt.colorbar() plt.show() # input('PRESS ENTER TO CONTINUE') mm=mm+1 #------------------------------- # import island # La matrice complete # data1='' # for kk in range(vv1.shape[0]) : # for jj in range(vv1.shape[1]) : # print kk,jj, "%1d"%(vv1[kk,jj]) # data1=data1+"%1d"%(vv1[kk,jj]) # data1=data1+"\n" # regs1 = island.locate_regions(data1) # print len(regs) # print_data1=island.print_regions(regs1,data1) # La matrice mise a zero hors zone de l'EUC vv2=np.copy(vv1) vv2[0:kEUC1,:]=0 vv2[kEUC2:-1,:]=0 vv2[:,0:jEUC1]=0 vv2[:,jEUC2:-1]=0 if debug : plt.figure(4000+mm) plt1a=plt.pcolormesh(vv2[:,:],cmap=cm_jet) # plt.ylim(yf,y0) # plt.xlim(35,45) # plt.ylim(54,10) ax = plt.gca() ax.invert_yaxis() plt.grid() plt.title(exp2_tex+' '+lonWBDY_str) plt.colorbar() plt.show() # input('PRESS ENTER TO CONTINUE') mm=mm+1 data2='' for kk in range(vv2.shape[0]) : for jj in range(vv2.shape[1]) : # print kk,jj, "%1d"%(vv2[kk,jj]) data2=data2+"%1d"%(vv2[kk,jj]) data2=data2+"\n" regs2 = island.locate_regions(data2) # print 'len(regs2)',len(regs2) print_data2=island.print_regions(regs2,data2,debug) # print "\n\n\n" if len(regs2) > 8 : print "\n\n Le nombre d'iles (",str(len(regs)),") ne doit pas etre > a 9 sinon 2 digits => EXIT\n\n" ; sys.exit(0) incrn=0 vv3=np.zeros(vv2.shape) for kk in range(vv2.shape[0]) : for jj in range(vv2.shape[1]) : incr=incrn+kk*vv1.shape[1]+jj if incr < len(print_data2)-1 : # print kk,jj, incr, print_data2[incr] if print_data2[incr]=='\n' : incrn=incrn+1 # print 'test',incr elif print_data2[incr]=='.' : vv3[kk,jj]=0 else : vv3[kk,jj]=float(print_data2[incr]) # raw_input('PRESS ENTER TO CONTINUE') #---------------------------------------------- # MAX dans la zone de l'EUC #---------------------------------------------- UUU_WBDY_EUC=ma.copy(UUU_WBDY*vv2) kUUmax_EUC,jUUmax_EUC=np.argwhere(UUU_WBDY_EUC == UUU_WBDY_EUC.max())[0] TTUUmax_EUC=TTT_WBDY[kUUmax_EUC,jUUmax_EUC] TTUUmax_EUC_str='\\Large Coeur EUC:[%.2f m/s %.1f C]'%(UUU_WBDY_EUC.max(),TTUUmax_EUC) TTUUmax_list.append(TTUUmax_EUC) iEUC=vv3[kUUmax_EUC,jUUmax_EUC] print " EUC max : j=",jUUmax_EUC," k=",kUUmax_EUC, "... indice:",iEUC,"\n" vv3[vv3!=iEUC]=0 vv3[vv3==iEUC]=1 vv33=ma.masked_values(vv3,0) #---------------------------------------------- # ponderate mean computation T & U #---------------------------------------------- TTT_WBDY3=vv33*TTT_WBDY pond_TTT_WBDY3=pond_WBDY[:,:]*vv33 cell_TTT_WBDY3=pond_TTT_WBDY3*TTT_WBDY3 pond_mean_TTT_WBDY3=cell_TTT_WBDY3.sum()/pond_TTT_WBDY3.sum() pond_mean_TTT_WBDY.append(pond_mean_TTT_WBDY3) UUU_WBDY3=vv33[0:zmax+un,:]*UUU_WBDY[0:zmax+un,:] pond_UUU_WBDY3=pond_WBDY[0:zmax+un,:]*vv33[0:zmax+un,:] cell_UUU_WBDY3=pond_UUU_WBDY3*UUU_WBDY3 pond_mean_UUU_WBDY3=cell_UUU_WBDY3.sum()/pond_UUU_WBDY3.sum() if B_section_plot : #----------------------------------------------- # fig 11 : uo Plot 0-300m #----------------------------------------------- plt.figure(11) UUU_WBDY2=UUU_WBDY[0:zmax+un,:] plt.pcolor(gphit4[:],-deptht4[0:zmax+un],UUU_WBDY2, cmap=cm_bwr) plt.clim(-Umin_max,Umin_max) plt.colorbar(extend='both') plt.xlim(section_latlim1,section_latlim2) plt.xlabel('Latitude') plt.ylim(depthlim1,depthlim2) plt.ylabel('depth (m)') plt.grid() cs11=plt.contour(gphit4[:],-deptht4[0:zmax+un],UUU_WBDY2, [-0.4,0,0.4],linewidths=linewidth1,colors='k') # plt.clabel(cs11) plt.plot([0,0], [0,-deptht4[zmax+un]], 'k', lw=1,linestyle='dashed') # Min Max mmUUU_WBDY2="\\Large (min:%.2f max:%.2f)" %(UUU_WBDY2.min(),UUU_WBDY2.max()) plt.text(section_latlim1,depthlim1,mmUUU_WBDY2,verticalalignment='bottom',horizontalalignment='left') # Max coeur EUC plt.plot(gphit4[jUUmax],-deptht4[kUUmax],'go') plt.text(section_latlim1,depth_text_max_uo,'\\textbf{'+TTUUmax_str+'}',verticalalignment='bottom',horizontalalignment='left',color='g') # Max uo sur tout le domaine plt.plot(gphit4[jUUmax_EUC],-deptht4[kUUmax_EUC],'mo') plt.text(section_latlim1,depth_text_coeur_EUC,'\\textbf{'+TTUUmax_EUC_str+'}',verticalalignment='bottom',horizontalalignment='left',color='m') # which ocean plt.text(section_latlim2,depthlim2-5,'\\Large \\textbf{'+ocean+'}',verticalalignment='top',horizontalalignment='right',color='k') # section from PERU/TROP plt.text(section_latlim2,depthlim1,"\\Large Section from "+section_from_domain,verticalalignment='bottom',horizontalalignment='right') # image created the... plt.text(section_latlim2,depthlim1+30,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title & save title11='\\large \\textbf{'+exp2_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{uo} (m/s) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title11) # if debug : plt.show() else : plt.savefig(img_dir+'/'+exp2+'_'+lonWBDY_str+'_from_'+section_from_domain+'_uo_0-300m.'+img) plt.close(11) #------------------------------------------------------------------------------- # fig 111 : uo Plot de la selection de la zone forte : 0-300m #------------------------------------------------------------------------------- plt.figure(111) plt.pcolor(gphit4[:],-deptht4[0:zmax+un],UUU_WBDY3, cmap=cm_bwr) plt.clim(-Umin_max,Umin_max) plt.colorbar(extend='both') plt.xlim(section_latlim1,section_latlim2) plt.xlabel('Latitude') plt.ylim(depthlim1,depthlim2) plt.ylabel('depth (m)') plt.grid() cs111=plt.contour(gphit4[:],-deptht4[0:zmax+un],UUU_WBDY2, [-0.4,0,0.4],linewidths=linewidth1,colors='k') # plt.clabel(cs111) plt.plot([0,0], [0,-deptht4[zmax+un]], 'k', lw=1,linestyle='dashed') mmUUU_WBDY3="\\Large (min:%.2f max:%.2f)" %(UUU_WBDY3.min(),UUU_WBDY3.max()) # Min Max plt.text(section_latlim1,depthlim1,mmUUU_WBDY3,verticalalignment='bottom',horizontalalignment='left') # Max coeur EUC plt.plot(gphit4[jUUmax],-deptht4[kUUmax],'go') plt.text(section_latlim1,depth_text_max_uo,'\\textbf{'+TTUUmax_str+'}',verticalalignment='bottom',horizontalalignment='left',color='g') # Max uo sur tout le domaine plt.plot(gphit4[jUUmax_EUC],-deptht4[kUUmax_EUC],'mo') plt.text(section_latlim1,depth_text_coeur_EUC,'\\textbf{'+TTUUmax_EUC_str+'}',verticalalignment='bottom',horizontalalignment='left',color='m') # Selection % of uo plt.text(section_latlim1,depth_text_seuil_EUC,seuil_str,verticalalignment='bottom',horizontalalignment='left') # uo ponderate mean plt.text(section_latlim1,0,'\\textbf{%.2f'%(pond_mean_UUU_WBDY3)+'m/s}',fontsize=48,color='magenta', verticalalignment='top',horizontalalignment='left') # section from PERU/TROP plt.text(section_latlim2,depthlim1,"\\Large Section from "+section_from_domain,verticalalignment='bottom',horizontalalignment='right') # which ocean plt.text(section_latlim2,depthlim2-5,'\\Large \\textbf{'+ocean+'}',verticalalignment='top',horizontalalignment='right',color='k') # image created the... plt.text(section_latlim2,depthlim1+30,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title & save title111='\\large \\textbf{'+exp2_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{uo} (m/s) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title111) if debug : plt.show() else : plt.savefig(img_dir+'/'+exp2+'_'+lonWBDY_str+'_from_'+section_from_domain+'_uo_select_EUC_%.0f_pourc_0-300m.'%(seuil_uo_EUC)+img) plt.close(111) #------------------------------------------------------------------------------- # fig 1111 : uo Plot de la selection de la zone forte : 75 levels #------------------------------------------------------------------------------- plt.figure(1111) UUU_WBDY4=vv33*UUU_WBDY plt.pcolor(gphit4[:],np.arange(75),UUU_WBDY4[:,:], cmap=cm_bwr) ax = plt.gca() ax.invert_yaxis() plt.clim(-Umin_max,Umin_max) plt.colorbar(extend='both') plt.xlim(section_latlim1,section_latlim2) plt.xlabel('Latitude') plt.ylim(75,0) plt.ylabel('Level') plt.grid() plt.plot([0,0], [0,75], 'k', lw=1,linestyle='dashed') cs1111=plt.contour(gphit4[:],np.arange(75),UUU_WBDY[:,:], [-0.4,0,0.4],linewidths=linewidth1,colors='k') # plt.clabel(cs111) # Min Max plt.text(section_latlim1,75,mmUUU_WBDY3,verticalalignment='bottom',horizontalalignment='left') # Max coeur EUC plt.plot(gphit4[jUUmax],kUUmax,'go') plt.text(section_latlim1,k_text_max_uo,'\\textbf{'+TTUUmax_str+'}',verticalalignment='bottom',horizontalalignment='left',color='g') # Max uo sur tout le domaine plt.plot(gphit4[jUUmax_EUC],kUUmax_EUC,'mo') plt.text(section_latlim1,k_text_coeur_EUC,'\\textbf{'+TTUUmax_EUC_str+'}',verticalalignment='bottom',horizontalalignment='left',color='m') # Selection % of uo plt.text(section_latlim1,k_text_seuil_EUC,seuil_str,verticalalignment='bottom',horizontalalignment='left') # uo ponderate mean plt.text(section_latlim1,0,'\\textbf{%.2f'%(pond_mean_UUU_WBDY3)+'m/s}',fontsize=48,color='magenta', verticalalignment='top',horizontalalignment='left') # section from PERU/TROP plt.text(section_latlim2,75,"\\Large Section from "+section_from_domain,verticalalignment='bottom',horizontalalignment='right') # which ocean plt.text(section_latlim2,1,'\\Large \\textbf{'+ocean+'}',verticalalignment='top',horizontalalignment='right',color='k') # image created the... plt.text(section_latlim2,70,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title & save title1111='\\large \\textbf{'+exp2_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{uo} (m/s) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title1111) if debug : plt.show() else : plt.savefig(img_dir+'/'+exp2+'_'+lonWBDY_str+'_from_'+section_from_domain+'_uo_select_EUC_%.0f_pourc_75_levels.'%(seuil_uo_EUC)+img) plt.close(1111) #------------------------------------------------------------------------------- # Section WBDY - EUC thetao (C) #------------------------------------------------------------------------------- #------------------------------------------------------------------------------- # fig 12 : thetao Plot 0-300m #------------------------------------------------------------------------------- plt.figure(12) v=np.linspace(Tmin_WBDY,Tmax_WBDY,Tmax_WBDY-Tmin_WBDY+1, endpoint=True) TTT_WBDY2=TTT_WBDY[0:zmax+un,:] plt.pcolor(gphit4[:],-deptht4[0:zmax+un],TTT_WBDY2, cmap=cm_jet) plt.clim(Tmin_WBDY,Tmax_WBDY) plt.colorbar(ticks=v,extend='both') plt.xlim(section_latlim1,section_latlim2) plt.xlabel('Latitude') plt.ylim(depthlim1,depthlim2) plt.ylabel('depth (m)') plt.grid() cs12=plt.contour(gphit4[:],-deptht4[0:zmax+un],UUU_WBDY2, [-0.4,0,0.4],linewidths=linewidth1,colors='k') # plt.clabel(cs12) plt.plot([0,0], [0,-deptht4[zmax+un]], 'k', lw=1,linestyle='dashed') # Min Max mmTTT_WBDY2="\\Large (min:%.1f max:%.1f)" %(TTT_WBDY2.min(),TTT_WBDY2.max()) plt.text(section_latlim1,depthlim1,mmTTT_WBDY2,verticalalignment='bottom',horizontalalignment='left') # Max coeur EUC plt.plot(gphit4[jUUmax],-deptht4[kUUmax],'go') plt.text(section_latlim1,depth_text_max_uo,'\\textbf{'+TTUUmax_str+'}',verticalalignment='bottom',horizontalalignment='left',color='g') # Max uo sur tout le domaine plt.plot(gphit4[jUUmax_EUC],-deptht4[kUUmax_EUC],'mo') plt.text(section_latlim1,depth_text_coeur_EUC,'\\textbf{'+TTUUmax_EUC_str+'}',verticalalignment='bottom',horizontalalignment='left',color='m') # section from PERU/TROP plt.text(section_latlim2,depthlim1,"\\Large Section from "+section_from_domain,verticalalignment='bottom',horizontalalignment='right') # which ocean plt.text(section_latlim2,depthlim2-5,'\\Large \\textbf{'+ocean+'}',verticalalignment='top',horizontalalignment='right',color='k') # image created the... plt.text(section_latlim2,depthlim1+30,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title & save title12='\\large \\textbf{'+exp2_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{thetao} (C) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title12) if debug : plt.show() else : plt.savefig(img_dir+'/'+exp2+'_'+lonWBDY_str+'_from_'+section_from_domain+'_thetao_0-300m.'+img) plt.close(12) #------------------------------------------------------------------------------- # fig 122 : thetao Plot de la selection de la zone forte : 0-300m #------------------------------------------------------------------------------- plt.figure(122) v=np.linspace(Tmin_WBDY_seuil,Tmax_WBDY_seuil,Tmax_WBDY_seuil-Tmin_WBDY_seuil+1, endpoint=True) plt.pcolor(gphit4[:],-deptht4[0:zmax+un],TTT_WBDY3[0:zmax+un,:], cmap=cm_jet) plt.clim(Tmin_WBDY_seuil,Tmax_WBDY_seuil) plt.colorbar(ticks=v,extend='both') plt.xlim(section_latlim1,section_latlim2) plt.xlabel('Latitude') plt.ylim(depthlim1,depthlim2) plt.ylabel('depth (m)') plt.grid() cs122=plt.contour(gphit4[:],-deptht4[0:zmax+un],UUU_WBDY2, [-0.4,0,0.4],linewidths=linewidth1,colors='k') # plt.clabel(cs122) plt.plot([0,0], [0,-deptht4[zmax+un]], 'k', lw=1,linestyle='dashed') # Min Max mmTTT_WBDY3="\\Large (min:%.1f max:%.1f)" %(TTT_WBDY3.min(),TTT_WBDY3.max()) plt.text(section_latlim1,depthlim1,mmTTT_WBDY3,verticalalignment='bottom',horizontalalignment='left') # Max coeur EUC plt.plot(gphit4[jUUmax],-deptht4[kUUmax],'go') plt.text(section_latlim1,depth_text_max_uo,'\\textbf{'+TTUUmax_str+'}',verticalalignment='bottom',horizontalalignment='left',color='g') # Max uo sur tout le domaine plt.plot(gphit4[jUUmax_EUC],-deptht4[kUUmax_EUC],'mo') plt.text(section_latlim1,depth_text_coeur_EUC,'\\textbf{'+TTUUmax_EUC_str+'}',verticalalignment='bottom',horizontalalignment='left',color='m') # Selection % of uo plt.text(section_latlim1,depth_text_seuil_EUC,seuil_str,verticalalignment='bottom',horizontalalignment='left') # thetao & uo ponderate mean plt.text(section_latlim1+0.5,depthlim2-5,'\\textbf{%.1f'%(pond_mean_TTT_WBDY3)+'C',fontsize=48,color='blue', verticalalignment='top',horizontalalignment='left') plt.text(section_latlim1+0.5,depthlim1+10,'\\textbf{%.2f m/s}'%(pond_mean_UUU_WBDY3),fontsize=48,color='magenta', verticalalignment='bottom',horizontalalignment='left') # section from PERU/TROP plt.text(section_latlim2,depthlim1,"\\Large Section from "+section_from_domain,verticalalignment='bottom',horizontalalignment='right') # which ocean plt.text(section_latlim2,depthlim2-5,'\\Large \\textbf{'+ocean+'}',verticalalignment='top',horizontalalignment='right',color='k') # image created the... plt.text(section_latlim2,depthlim1+30,'\\large \\textbf{\\it{'+local_time4+'}}',verticalalignment='bottom',horizontalalignment='right',color='black') # Title & save title122='\\large \\textbf{'+exp2_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{thetao} (C) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title122) if debug : plt.show() else : plt.savefig(img_dir+'/'+exp2+'_'+lonWBDY_str+'_from_'+section_from_domain+'_thetao_select_EUC_%.0f_pourc_0-300m.'%(seuil_uo_EUC)+img) plt.close(122) #------------------------------------------------------------------------------- # Fill the mean temperature & flux ascii file (out2) #------------------------------------------------------------------------------- text2c='' if B_section_WBDY | B_section_TROP : for ii in range(len(UUU_WBDY_list)) : if uo_max : # print la Temperature au point du max uo de l'EUC dans le fichier de sortie EUC_Upwelling_characterize_output_SENEGAL.txt text2c=text2c+'%-7.1f'%(TTUUmax_list[ii]) else : # print la Temperature ponderee sur le coeur de l'EUC ... text2c=text2c+'%-7.1f'%(pond_mean_TTT_WBDY[ii]) if B_hfdsO : text2c=text2c+\ '%-7.1f'%(pond1_mean_hfdsO1_upA)+\ '%-7.1f'%(pond1_mean_hfdsO1_upB)+\ '%-7.1f'%(pond1_mean_hfdsO1_upBas) if B_T_iso : for ll in range(len(pond1_mean_TT1_iso_upA_list)) : text2c=text2c+\ '%-7.1f'%(pond1_mean_TT1_iso_upA_list[ll])+\ '%-7.1f'%(pond1_mean_TT1_iso_upB_list[ll])+\ '%-7.1f'%(pond1_mean_TT1_iso_upBas_list[ll]) text2c=text2c+'%-6s'%(id)+'%-37s'%(exp1)+'\n' fout2.write(text2c) #------------------------------- plt.close(11) # plt.close(110) plt.close(111) plt.close(1111) plt.close(12) plt.close(120) plt.close(122) plt.close(13) plt.close(130) plt.close(133) plt.close(14) plt.close(140) plt.close(143) plt.close(15) plt.close(150) plt.close(16) plt.close(160) plt.close(21) plt.close(22) plt.close(23) plt.close(24) #------------------------------------------------------------------------------- # close the mean temperature & flux ascii file (out2) #------------------------------------------------------------------------------- if (B_section_TROP | B_section_WBDY) & (B_T_iso | B_hfdsO) : text2b=\ "cmip5_res='%3s'\n"%(cmip5_res)+\ "cmip5_2D_interp_method='%s'\n"%(cmip5_2D_interp_method)+\ "Zone upA: longitude: %.4f / %.4f\n"%(lonupA1,lonupA2)+\ "Zone upA: latitude: %.4f / %.4f\n"%(latupA1,latupA2)+\ "Zone upB: longitude: %.4f / %.4f\n"%(lonupB1,lonupB2)+\ "Zone upB: latitude: %.4f / %.4f\n"%(latupB1,latupB2)+\ "ZOOM pour detection max EUC: latitude: %.4f / %.4f\n"%(gphit1[jEUC1],gphit1[jEUC2])+\ "ZOOM pour detection max EUC: depth: %.4f / %.4f\n"%(deptht1[kEUC1],deptht1[kEUC2])+\ "Section from "+domain_MAJ+": "+str(yb1)+'-'+str(ye1)+'\n'+\ "Section from TROP: "+str(yb1)+'-'+str(ye1)+'\n' if uo_max : text2b=text2b+\ "Temperature des sections: valeur au point uo max \n" else : text2b=text2b+\ "Temperature des sections: valeur moyenne du coeur de l'EUC: seuil %.0f"%(seuil_uo_EUC)+"% de uo max\n" fout2.write(text2b) fout2.close() # shutil.move(out2,scriptname_root+'_output.txt')