#!/usr/bin/env python # 2019 (C. Hourdin) # 19 Juin 2023 (C. Hourdin) # # plot sections / Upwelling areas # mean of # #------------------------------------------------------------------------------- # Infos #------------------------------------------------------------------------------- # upA : A upwelling area # upB : B upwelling area # upBas1 : B upwelling area - along shore # upBas2 : B upwelling area - along shore # lon, lat, depth # i, j, k # thetao: TT, Tmin, Tmax # uo : UU, Umin, Umax ### areas: # 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 #=============================================================================== # Choices #=============================================================================== debug=False # ATTENTION: si True pas dimpression images et taper pylab si on veut voir les plot de debug #pylab resolution='025' # 'tr075' or 'tr025' : Pulsation Tropical grid resolution # 'reg1x1' : regular grid 1x1 degrees #project_list=['cmip5','cmip6', 'trop', 'GLORYS', 'soda3'] # cmip5, cmip6, trop, GLORYS, soda3 project_list=['cmip6'] #exp1_list="['GLORYS2V3_ORCA025_1993-2009_R20130808_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip5_ACCESS1-0_historical_19890101-20051231_Omon_1m_r1i1p1_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip5_BNU-ESM_historical_19890101-20051231_Omon_1m_r1i1p1_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip5_CCSM4_historical_19890101-20051231_Omon_1m_r1i1p1_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip5_NorESM1-M_historical_19890101-20051231_Omon_1m_r1i1p1_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip5_GISS-E2-R-CC_historical_19890101-20051231_Omon_1m_r1i1p1_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip5_GISS-E2-R_historical_19890101-20051231_Omon_1m_r1i1p1_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip6_CESM2_historical_19890101-20051231_Omon_1m_r1i1p1f1_gn_thetao_mean.grid_tr025-lev1_orig.nc']" exp1_list="['cmip6_CIESM_historical_19890101-20051231_Omon_1m_r1i1p1f1_gn_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip6_MRI-ESM2-0_historical_19890101-20051231_Omon_1m_r1i1p1f1_gn_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['cmip6_CESM2-WACCM-FV2_historical_19890101-20051231_Omon_1m_r1i1p1f1_gn_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['trop075_now-tr075now_long_sol094_c1m_1980_2012_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['trop025_nemo-tr025_blp_c1m_2001_2008_thetao_mean.grid_tr025-lev1_orig.nc','trop075_nemo-tr075_erainocrt_c1m_2001_2008_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['trop025_nemo-tr025_blp_c1m_2001_2008_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['trop12_nemo-tr12_quik_c1m_2001_2008_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="['trop025_nemo-tr025_blp_c1m_2001_2008_thetao_mean.grid_tr025-lev1_orig.nc']" #exp1_list="sorted(glob.glob('*grid_tr'+resolution+'-lev1_orig.nc'))" #exp1_list="sorted(glob.glob('cmip5_HadGEM2*_thetao_mean.grid_tr025-lev1_orig.nc'))" #exp1_list="sorted(glob.glob('*tr12n_n0*'))" #exp1_list="sorted(glob.glob('trop12*tr'+resolution+'*'))" #exp1_list="\ #glob.glob('cmip6_EC-Earth3*')+\ #glob.glob('cmip6_EC-Earth3-Veg*')+\ #glob.glob('cmip6_TaiESM1*')" #exp1_list="sorted(glob.glob('*tr'+resolution+'*'))" domain='peru'; domain_MAJ='PERU'; #domain='benguela'; domain_MAJ='BENGUELA'; #domain='senegal'; domain_MAJ='SENEGAL'; #domain='california'; domain_MAJ='CALIFORNIA'; kdim=75 lonlatdiffseuil=5e-06 print('debug :', debug) print('project_list :', project_list) print('resolution :', resolution) print('domain :', domain_MAJ) #------------------------------------------------------------------------------- # Plot #------------------------------------------------------------------------------- linewidth1=0.5 over1='y' under1='m' un=1 ii5c=0 ii6c=0 ii25n=0 ii25c=0 ii75n=0 ii75c=0 ii12n=0 iiglo=0 iisod=0 # Image type (none, png, eps) img='png' #img='png' # necessite mpl.use('Agg') #img='eps' # necessite mpl.use('ps') #img='pdf' # necessite mpl.use('pdf') if debug: img='none' mm=0 # beginning of plot increment number #----------------------------------------------------------------- # flags a nettoyer #----------------------------------------------------------------- F_dxdy=True # Need SST file on native grid dx=0.; dy=0. # Pour éviter un value undefined si grid native pas dispo # coordonnées ou on regarde la résolution : lon_dxdy=10; lon_dxdy_str='10°S' # longitude 0° lat_dxdy=-110; lat_dxdy_str='110W' # latitude 110W #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 # ATTENTION : 2 digits! sinon "Modifier iWBDY ou section_PERU" (cf test plus loin) #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... #---------------------------------- ########################### ### DEFAULT VALUES # don't touch !!!! B_hfdsO_plot=False B_isodepth_plot=False B_hfdsO=False B_section_plot=False hfdsO_exists=True ### DEFAULT VALUES ########################### #--------------------------------------- # sections (uo, thetao) - from 3D files #--------------------------------------- B_section_TROP=False # mask2 (fmask2), TT2 (fT2) & UU2 (fU2) : tmask, thetao & uo read in TROP sections files B_section_WBDY=False # 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 if B_section_TROP | B_section_WBDY : B_section_plot=False #------------------------------------------------------ # isodepth (thetao from 3D files / hfdsO from 2D file) #------------------------------------------------------ B_isodepth=True # TT1 (fT1) : thetao read in domain files (peru, benguela...) - files extracted from trop files or generated by the regional model if B_isodepth : B_isodepth_plot=True # thetao=f(depth) B_hfdsO=True if B_hfdsO : B_hfdsO_plot=True #-------------------------------------------------------- # cmip5 NATIVE GRID (_ng) separate plot (sst, Flux): #-------------------------------------------------------- #ng_Flux_list=['hfdsO','hfdsA'] ng_Flux_list=['hfdsO'] D_native_grid_sst_plot=False # read the 3D thetao cmip file # TT_ng (fT_ng) : thetao read in domain files (peru, benguela...) - files extracted from cmip5 NATIVE files D_native_grid_Flux_plot=False # read the 2D Heat Flux cmip5 file # hfds_ng (fH_ng) : flux read in domain files (peru, benguela...) - files extracted from cmip5 NATIVE files ### CHRIS... ne faut il pas retirer B_iso et executer tout le temps ce qui est sous la clé??? B_iso= (B_isodepth | B_section_WBDY | D_native_grid_sst_plot | D_native_grid_Flux_plot) #------------------------------------------------------------------------------- # Sections #------------------------------------------------------------------------------- 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 seuil_uo_EUC=20 # % du Max de uo de l'EUC sur lequel on calcul la moyennne de temperature if (domain == 'peru') | (domain == 'california') : ocean='PACIFIC' section_list=['85W','98.5W','110.5W','120.25W','130W','140.5W','170.5W'] # section_list=['85W','98.5W','110.5W'] # section_list=['110.5W','120.25W','130W'] 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'] # Pour sections TROP pacifique et atlantique # profondeur pour carte lon/lat thetao kEUC1=10; kEUC2=40 zmax=34 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 ,2026 ,29 ] lev_list= ['sst','54m','87m','120m','31-108m','69-133m','181m'] #knum =[0 ,1 ,2 ,3 ,4 ,5 ,6 ] knum =[0,1,2,3,4,5,6] #k_iso_list= [0 ,14 ,22 ,24 ,1424 ,2026 ,29 ] #lev_list= ['sst','31m','87m','108m','31-108m','69-133m','181m'] #k_iso_list= [0 ,20 ,22 ,26 ,1424 ,2026 ,29 ] #lev_list= ['sst','69m','87m','133m','31-108m','69-133m','181m'] #knum =[0 ,1 ,2 ,3 ,4 ,5 ,6 ] #knum =[0,2,4,5,6] ### Pour test de sensibilité en profondeur (sur PERU.txt et section.txt) ### Ne fonctionne pas pour les plots sans définir les échelles de couleur plus bas ### Tmin_iso_upBas_list & Tmin_iso_list ### Tmax_iso_upBas_list & Tmax_iso_list ### avant lev 12 (23m) risque de NaN ###k_iso_list=[0 ,8 ,9 ,10 ,11 ,12 ,13 ,14 ,15 ,16 ,17 ,18 ,19 ,20 ,21 ,22 ,23 ,24 ,25 ,26 ,27 ,28 ,29 ,1424 ,2026 ] ###lev_list= ['sst','12m','14m','17m','19m','23m','27m','31m','36m','41m','47m','54m','61m','69m','78m','87m','97m','108m','120m','133m','147m','163m','181m','31-108m','69-133m'] ###knum =[0 ,1 ,2 ,3 ,4 ,5 ,6 ,7 ,8 ,9 ,10 ,11 ,12 ,13 ,14 ,15 ,16 ,17 ,18 ,19 ,20 ,21 ,22 ,23 ,24 ] #k_iso_list=[0 ,12 ,13 ,14 ,15 ,16 ,17 ,18 ,19 ,20 ,21 ,22 ,23 ,24 ,25 ,26 ,27 ,28 ,29 ,1424 ,2026 ] #lev_list= ['sst','23m','27m','31m','36m','41m','47m','54m','61m','69m','78m','87m','97m','108m','120m','133m','147m','163m','181m','31-108m','69-133m'] #knum =[0 ,1 ,2 ,3 ,4 ,5 ,6 ,7 ,8 ,9 ,10 ,11 ,12 ,13 ,14 ,15 ,16 ,17 ,18 ,19 ,20 ] #k_iso_list=[0 ,1424 , 2026 ] #lev_list= ['sst','31-108m','69-133m'] #knum =[0 ,1 ,2 ] upB_as_size1=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) upB_as_size2=8 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,12,10] Tmax_iso_upBas_list=[27,20,18,16,19,18,14] Tmin_iso_list=[15,12,12,11,11,11,10] Tmax_iso_list=[29,25,23,21,24,23,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,13,9] Tmax_iso_upBas_list=[28,22,19,16,25,19,13] Tmin_iso_list=[15,11,10,9,10,10,8] Tmax_iso_list=[28,22,19,17,26,19,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,14,10] Tmax_iso_upBas_list=[29,22,19,16,25,19,13] Tmin_iso_list=[21,15,13,12,15,13,11] Tmax_iso_list=[29,24,20,16,25,20,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,9,7] T=max_iso_upBas_list=[16,15,13,11,14,13,9] Tmin_iso_list=[12,10,9,8,10,9,7] Tmax_iso_list=[25,23,22,20,22,21,18] HFDS_min=0 HFDS_max=140 HFDSz_min=40 HFDSz_max=120 #=============================================================================== # 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 glob 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 sys.path.insert(0,'./EUC_Upwelling_characterize/modules') import island import xarray as xr #------------------------------------------------------------------------------- # Fonts #------------------------------------------------------------------------------- plt.rc('font',size=10) mpl.rc('text', usetex=True) # False pour desactiver latex. mpl.rc('font', family='serif') cm_bwr=plt.cm.get_cmap("bwr").copy() cm_bwr.set_over(over1) cm_bwr.set_under(under1) cm_jet=plt.cm.get_cmap("jet").copy() cm_jet.set_over(over1) cm_jet.set_under(under1) #------------------------------------------------------------------------------- # Divers #------------------------------------------------------------------------------- scriptname=os.path.basename(sys.argv[0]) scriptname_root=scriptname.split('.')[0] script_dir=os.getcwd() computer=str.replace(os.popen('hostname').readlines()[0], '\n', '') 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' img_dir_section_name='section_EUC' img_dir_area=script_dir+'/EUC_Upwelling_characterize/images/'+domain_MAJ img_dir_native=script_dir+'/EUC_Upwelling_characterize/images/'+domain_MAJ+'_native_grid' img_dir_section=script_dir+'/EUC_Upwelling_characterize/images/section' if uo_max : txt_out_dir=script_dir+'/EUC_Upwelling_characterize/outputs_EUC_uo_max' else : txt_out_dir=script_dir+'/EUC_Upwelling_characterize/outputs_EUC_uo_seuil_'+str(seuil_uo_EUC)+'pc' if B_hfdsO : txt_out_dir=txt_out_dir+'_hfdsO' dir1_list='[txt_out_dir' if B_isodepth_plot | B_hfdsO_plot : dir1_list=dir1_list+', img_dir_area' if B_section_plot : dir1_list=dir1_list+', img_dir_section' if D_native_grid_sst_plot | D_native_grid_Flux_plot : dir1_list=dir1_list+', img_dir_native' dir1_list=dir1_list+']' if not(debug) : for dir1 in eval(dir1_list) : if (os.path.exists(dir1)) : shutil.move(dir1,dir1+'_OLD_'+local_time3) os.makedirs(dir1,0o755) if debug : print('\nBEGIN Read fields at',datetime.datetime.strftime(datetime.datetime.now(), '%H:%M:%S\n')) # 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_upBas1,pond99_upBas2,exp99_tex,img_name99,depth_iso_str99,dx99,dy99) : #------------------------------------------------------------------------------- #----------------- # Plot full domain #----------------- lon99T,lat99T=flon_Tcentered(lon99,lat99) plt.figure(ifig99) plt.pcolor(lon99T,lat99T,var99,cmap=cm99,shading='auto') plt.clim(min99,max99) diff99=max99-min99 if unit99 == 'C' : # temperature if diff99 < 7 : deltav99=np.linspace(min99,max99,int(diff99*2+1), endpoint=True) else : deltav99=np.linspace(min99,max99,int(diff99+1), endpoint=True) if 'Flux Net Down' in var99_str : deltav99=np.linspace(min99,max99,int(diff99/20+1), endpoint=True) print("full domain : min=%.1f, max=%.1f, ticks:"%(min99,max99),deltav99) # print("CHRISTOPHE ", pond99_upBas1, pond99_upBas2) 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='black',facecolor='none',linestyle='-',lw=1,hatch='\\')) plt.text(lonupA1,latupA2,'Zone upA', verticalalignment='top',horizontalalignment='left',color='black') # latlim 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,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='black',verticalalignment=vertA,horizontalalignment='right') # upB pond mean plt.text(lonlim2,latlimB+1,'\\textbf{%.1f'%(pond99_upB)+' '+unit99,fontsize=48,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='blue',verticalalignment=vertB,horizontalalignment='right') # native_dxdy plt.text(-100,-10,'\\textbf{dx %.2f}'%(dx99), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='bottom',horizontalalignment='left') plt.text(-100,-10,'\\textbf{dy %.2f}'%(dy99), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='top',horizontalalignment='left') plt.text(-100,-15,'(at '+lon_dxdy_str+' '+lat_dxdy_str+')', fontsize=12,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='black',verticalalignment='bottom',horizontalalignment='left') # 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',bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='black') # domain plt.text(lonlim1,latlim2-1.5,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonlim1,latlim2-3.5,'\\Large \\textbf{grid tr'+resolution+'}',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 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_dir_area+'/'+img_name99) if img == 'png' : img_dir=img_dir_area img_name=img_name99 with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(ifig99) #--------------------- # Plot Zone upB along shore #--------------------- var99_upBas1=var99[jupB1:jupB2+un:,iupB1:iupB2+un]*maskupBas1_surf var99_upBas2=var99[jupB1:jupB2+un:,iupB1:iupB2+un]*maskupBas2_surf img_name99_ZOOM=img_name99.replace('.'+img,'_ZOOM.'+img) plt.figure(ifig99+1) plt.pcolor(lon99T[iupB1:iupB2+un],lat99T[jupB1:jupB2+un],var99_upBas1, cmap=cm99,edgecolors='blue', linewidth=0.3, linestyle='solid',shading='auto') plt.clim(min_upBas99,max_upBas99) plt.pcolor(lon99T[iupB1:iupB2+un],lat99T[jupB1:jupB2+un],var99_upBas2, cmap=cm99,edgecolors='magenta', linewidth=0.5, linestyle='solid',alpha=0.7,shading='auto') 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,int(diff_upBas99*2+1), endpoint=True) else : deltav_upBas99=np.linspace(min_upBas99,max_upBas99,int(diff_upBas99+1), endpoint=True) if 'Flux Net Down' in var99_str : deltav_upBas99=np.linspace(min_upBas99,max_upBas99,int(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') plt.text(lonupB1,latupB2-0.5,'\\Large \\textbf{grid tr'+resolution+'}',color='black',verticalalignment='top',horizontalalignment='left') # upBas pond mean # plt.text(lonupB2,latupB2-0.5,'\\textbf{%.1f'%(pond99_upBas1)+' '+unit99,fontsize=48,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='blue',verticalalignment='top',horizontalalignment='left') color991='blue'; color992='magenta' if pond99_upBas1 > pond99_upBas2 : color991='magenta'; color992='blue' plt.text(lonupB2,latupB2-0.5,'\\textbf{%.1f'%(pond99_upBas1), fontsize=48,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color=color991,verticalalignment='top',horizontalalignment='right') plt.text(lonupB2,latupB2-1.5,'\\textbf{%.1f'%(pond99_upBas1-pond99_upBas2),fontsize=32,bbox=dict(facecolor='grey', edgecolor='grey',alpha=0.5),color=color991,verticalalignment='top',horizontalalignment='right') plt.text(lonupB1,latupB1+0.5,'\\textbf{%.1f'%(pond99_upBas2), fontsize=48,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color=color992,verticalalignment='bottom',horizontalalignment='left') # native_dxdy plt.text(-88.7,-10,'\\textbf{dx %.2f'%(dx99), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='bottom',horizontalalignment='left') plt.text(-88.7,-10,'\\textbf{dy %.2f'%(dy99), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='top',horizontalalignment='left') plt.text(-88.7,-11,'(at '+lon_dxdy_str+' '+lat_dxdy_str+')', fontsize=12,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='black',verticalalignment='bottom',horizontalalignment='left') # min,max sur upBas1 & upBas2 mm99_upBas1="\\Large (min:%.1f max:%.1f)" %(var99_upBas1.min(),var99_upBas1.max()) plt.text(lonupB2,latupB2,mm99_upBas1,verticalalignment='top',horizontalalignment='right',bbox=dict(facecolor='white', edgecolor='white', alpha=0.5),color=color991) mm99_upBas2="\\Large (min:%.1f max:%.1f)" %(var99_upBas2.min(),var99_upBas2.max()) plt.text(lonupB1,latupB1,mm99_upBas2,verticalalignment='bottom',horizontalalignment='left',bbox=dict(facecolor='white', edgecolor='white', alpha=0.5),color=color992) # 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_dir_area+'/'+img_name99_ZOOM) if img == 'png' : img_dir=img_dir_area img_name=img_name99_ZOOM with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(ifig99+1) #------------------------------------------------------------------------------- # Initialize the mean temperature & flux ascii files (out_area & out_section) #------------------------------------------------------------------------------- out_area=txt_out_dir+'/'+domain_MAJ+'.txt' out_section=txt_out_dir+'/'+img_dir_section_name+'.txt' #if (os.path.exists(out_area)) : shutil.move(out_area,out_area+'_OLD_'+local_time3) #if (os.path.exists(out_section)) : shutil.move(out_section,out_section+'_OLD_'+local_time3) fout_area =open(out_area, 'a') fout_section =open(out_section, 'a') #------------ # First line #------------ # section if B_section_WBDY : text2section='# %-12s'%('T_EUC') else : text2section='# ' if B_section_TROP : for ii in range(len(section_list)) : text2section=text2section+'%-12s'%('S_EUC ') # area text2area='# ' if B_hfdsO: text2area=text2area+'%-28s'%('hfdsO') if B_isodepth : for iik in range(len(knum)) : ik=knum[iik] text2area=text2area+'T_%-26s'%(lev_list[ik]) # section & area text2commonB='%-90s'%('exp') text2commonB=text2commonB+'%-6s'%('dx') text2commonB=text2commonB+'%-6s'%('dy') text2commonB=text2commonB+'%-6s'%('yb') text2commonB=text2commonB+'%-6s'%('ye') text2commonB=text2commonB+'%-6s'%('id') text2commonB=text2commonB+'\n# ' text2section=text2section+ text2commonB text2area=text2area+ text2commonB #------------ # Second line #------------ # section if B_section_WBDY : text2section=text2section+'%-12s'%(section_PERU) if B_section_TROP : for ii in range(len(section_list)) : text2section=text2section+'%-12s'%(section_list[ii]) text2section=text2section+text2commonB # area text2common='%-7s%-7s%-7s%-7s'%('upA','upB','upBas1','upBas2') if B_hfdsO: text2area=text2area+text2common if B_isodepth : for ll in range(len(knum)) : text2area=text2area+text2common text2area=text2area+text2commonB # section & area text2area=text2area+'\n\n' text2section=text2section+'\n\n' fout_area.write(text2area) fout_section.write(text2section) #------------------------------------------------------------------------------- # Area A & B Coordinates definition 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[21]=-84.25; lonT[33]=-75.25; latT[24]=-13.3767 ; latT[33]=-6.7344 iupB1=15; iupB2=33; jupB1=24; jupB2=33; # 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 # UpA centre sur l'upwelling cotier a l'ouest du domaine # iupA1=0; iupA2=13; jupA1=33; jupA2=46; # 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 #-------------------------------------- # 075 => coordinates for 025 & 12 #-------------------------------------- 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' : 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 #=============================================================================== # CHRIS ISO (mask1) #=============================================================================== if B_iso : #------------------------------------------------------------------------------- # mask1: read trop3D grid mesh_mask from domain file (peru, senegal...) #------------------------------------------------------------------------------- fmask1='/data/cholod/EUC_Upwelling/trop/trop'+resolution+'_nemo/inputs/mesh_mask_TROP'+resolution+'_'+domain+'.nc' print('\nflag B_iso :') 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)) w1a=e1t1_3D*e2t1_3D*e3t1 w1a_surf=e1t1*e2t1 #------------------------------------------------------------------------------- # For each project... #------------------------------------------------------------------------------- for project in project_list : #------------------------------------------------------------------------------- # For each exp1_file... #------------------------------------------------------------------------------- dir1='/data/cholod/EUC_Upwelling/'+project+'/05_lev1_orig/thetao' print('\ndir1 :', dir1) os.chdir(dir1) if debug : print ('directory : ',os.getcwd()) if debug : os.system('ls -1') for exp1_file in eval(exp1_list) : if ('p12' in exp1_file) | ('erai7' in exp1_file) : if B_section_TROP : print('\n\n\n exp BDY dans '+exp1_file+' => ATTENTION: TROP section < WBDY section indispensable \n\n\n') print('\n',80*'=','\n EXP1= ',exp1_file,' \n',80*'=') if 'cmip' in project : yb1= exp1_file.split('_Omon')[0][-17:-13] ye1= exp1_file.split('_Omon')[0][-8:-4] exp1_title = exp1_file.split('_historical_')[0] elif 'trop' in project : yb1= exp1_file.split('_thetao_mean')[0][-9:-5] ye1= exp1_file.split('_thetao_mean')[0][-4:] exp1_title= exp1_file.split('_thetao_mean')[0][0:-10] elif 'GLORYS' in project : yb1= exp1_file.split('_R20130808')[0][-9:-5] ye1= exp1_file.split('_R20130808')[0][-4:] exp1_title = exp1_file.split('_R20130808')[0][:] elif 'soda3' in project : yb1= exp1_file.split('_thetao_mean')[0][-9:-5] ye1= exp1_file.split('_thetao_mean')[0][-4:] exp1_title = exp1_file.split('_thetao_mean')[0][:] exp1_name= exp1_file.split('_thetao_mean')[0][:] exp1_tex=exp1_title.replace('_','\\_') print('exp1_name : ',exp1_name) print('exp1_title : ',exp1_title) print('exp1_tex : ',exp1_tex) print('yb1 : ',yb1) print('ye1 : ',ye1) #=============================================================================== # CHRIS ISO on cmip5 native grid (hfds_ng, TT_ng) # # only for cmip5 !!! # # hfds & sst #=============================================================================== if D_native_grid_sst_plot | (F_dxdy & (not('trop' in exp1_file))) : #------------------------------------------------------------------------------- # read SST NATIVE grid #------------------------------------------------------------------------------- # open & read NATIVE grid thetao file dir1_ng=dir1.replace('05_lev1_orig','02_pre-processing').replace('thetao','thetao_lon_linear') fT_ng=exp1_file.replace('grid_tr'+resolution+'-lev1_orig.nc','PRE-PROCESS.nc') ncT_ng=netcdf(dir1_ng+'/'+fT_ng,'r') lon_ng=ncT_ng.variables['lon'][:,:] print('lon_ng: de %7.1f a %7.1f'%(lon_ng.min(),lon_ng.max())) if ('CESM1' in exp1_file) | \ ('CCSM4' in exp1_file) |\ ('NorESM1' in exp1_file) |\ ('CESM2' in exp1_file) |\ ('CIESM' in exp1_file) |\ ('FIO' in exp1_file) |\ ('MCM' in exp1_file) |\ ('SAM0' in exp1_file) |\ ('TaiESM' in exp1_file) : lon_ng=lon_ng-720 else : lon_ng=lon_ng-360 lat_ng=ncT_ng.variables['lat'][:,:] deptht_ng=ncT_ng.variables['lev'][:] TT_ng=ncT_ng.variables['thetao'][:,:,:] ncT_ng.close() print('lon_ng: de %7.1f a %7.1f'%(lon_ng.min(),lon_ng.max())) print('\n',80*'-','\n on sst NATIVE GRID : \n',80*'-') print('domain '+domain_MAJ+' : reading thetao on sst NATIVE grid (TT_ng)...') print(' dir1_ng :', dir1_ng) print(' fT_ng :', fT_ng) if F_dxdy : #------------------------------------------------------------------------------- # read & write dxdy #------------------------------------------------------------------------------- if 'trop' in exp1_file : if 'trop12' in exp1_file : dx = 0.08 if 'trop025' in exp1_file : dx = 0.25 if 'trop075' in exp1_file : dx = 0.75 dy = dx else : # dx = lon_ng[0,-1]-lon_ng[0,-2] # dy = lat_ng[-1,0]-lat_ng[-2,0] diff_lat=np.absolute(lat_ng[:,0]-lon_dxdy) Yindex=diff_lat.argmin() diff_lon=np.absolute(lon_ng[Yindex,:]-lat_dxdy) Xindex=diff_lon.argmin() dx=abs(lon_ng[Yindex,Xindex]-lon_ng[Yindex,Xindex+1]) dy=abs(lat_ng[Yindex,Xindex]-lat_ng[Yindex+1,Xindex]) print(' dx : %.2f'%dx) print(' dy : %.2f'%dy) if D_native_grid_sst_plot : #------------------------------------------------------------------------------- # plot SST NATIVE grid #------------------------------------------------------------------------------- #----------------- # Plot SST full domain #----------------- print("PLOT sst (C) on native GRID depth=%-.1f\n"%(deptht_ng[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) lon_ngT,lat_ngT=flon_Tcentered(lon_ng,lat_ng) plt.pcolor(lon_ngT[:],lat_ngT[:],TT_ng[0,:,:],cmap=cm_jet,shading='auto') 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='black',facecolor='none',linestyle='-',lw=1,hatch='\\')) plt.text(lonupA1,latupA2,'Zone upA', verticalalignment='top',horizontalalignment='left',color='black') # domain plt.text(lonlim1,latlim2-1.5,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonlim1,latlim2-3.5,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='top',horizontalalignment='left') # native_dxdy plt.text(-100,-10,'\\textbf{dx %.2f'%(dx), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='bottom',horizontalalignment='left') plt.text(-100,-10,'\\textbf{dy %.2f'%(dy), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='top',horizontalalignment='left') plt.text(-100,-15,'(at '+lon_dxdy_str+' '+lat_dxdy_str+')', fontsize=12,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),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 ... exp1_file is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1_file is used for title (cmip5_IPSL-CM5A-LR_thetao...) because we plot the NATIVE GRID # depth_verif_str_ng="lev1 -%.0fm"%(deptht_ng[0]) depth_verif_str_ng=" -%.0fm"%(deptht_ng[0]) # title21='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{thetao} (C) \\textbf{'+depth_verif_str_ng+'}' title21='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{sst} (C) \\textbf{'+depth_verif_str_ng+'} (T point)' plt.title(title21) # save img if debug : plt.show() else : img_name=exp1_name+'_'+domain_MAJ+'_sst_native_grid.'+img plt.savefig(img_dir_native+'/'+img_name) if img == 'png' : img_dir=img_dir_native with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) #--------------------- # 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(lon_ngT[:],lat_ngT[:],TT_ng[0,:,:],cmap=cm_jet,edgecolors='blue', linewidth=0.3, linestyle='solid',shading='auto') 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() plt.text(lonupB1,latupB2-0.3,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonupB1,latupB2-0.6,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='top',horizontalalignment='left') # native_dxdy plt.text(-88.7,-10,'\\textbf{dx %.2f'%(dx), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='bottom',horizontalalignment='left') plt.text(-88.7,-10,'\\textbf{dy %.2f'%(dy), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='top',horizontalalignment='left') plt.text(-88.7,-11,'(at '+lon_dxdy_str+' '+lat_dxdy_str+')', fontsize=12,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),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 ... exp1_file is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1_file 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_str_ng+'} (T point)' plt.title(title22) # save img if debug : plt.show() else : img_name=exp1_name+'_'+domain_MAJ+'_sst_native_grid_ZOOM.'+img plt.savefig(img_dir_native+'/'+img_name) if img == 'png' : img_dir=img_dir_native with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) if D_native_grid_Flux_plot : for ng_Flux in ng_Flux_list : #------------------------------------------------------------------------------- # plot HFDS O NATIVE grid #------------------------------------------------------------------------------- # open & read NATIVE grid hfdsO file dir1_ng_hfds=dir1.replace('05_lev1_orig','02_pre-processing').replace('thetao','hfdsO_lon_linear') fH_ng=exp1_file.replace('grid_tr'+resolution+'-lev1_orig.nc','PRE-PROCESS.nc').replace('thetao','hfdsO') if os.path.exists(dir1_ng_hfds+'/'+fH_ng) : ncH_ng=netcdf(dir1_ng_hfds+'/'+fH_ng,'r') lon_ng=ncH_ng.variables['lon'][:,:] lat_ng=ncH_ng.variables['lat'][:] hfds_ng=ncH_ng.variables['hfds'][:,:] ncH_ng.close() print('\n',80*'-','\n on hfdsO NATIVE GRID : \n',80*'-') print(' domain '+domain_MAJ+' : reading hfds on hfdsO NATIVE grid (TT_ng)...') print('dir1_ng_hfds :', dir1_ng_hfds) print(' fH_ng :', fH_ng) ng_Flux_str=ng_Flux hfds_ng_B=np.copy(hfds_ng) #----------------- # Plot FLUX full domain #----------------- print('\nPLOT '+ng_Flux+' (W/m2) on native GRID' ) plt.figure(23) if ('CESM1' in exp1_file) | \ ('CCSM4' in exp1_file) |\ ('NorESM1' in exp1_file) |\ ('CESM2' in exp1_file) |\ ('CIESM' in exp1_file) |\ ('FIO' in exp1_file) |\ ('MCM' in exp1_file) |\ ('SAM0' in exp1_file) |\ ('TaiESM' in exp1_file) : lon_ng=lon_ng-720 else : lon_ng=lon_ng-360 lon_ngT,lat_ngT=flon_Tcentered(lon_ng,lat_ng) v=np.linspace(HFDS_max,HFDS_min,int((HFDS_max-HFDS_min)/20+1), endpoint=True) plt.pcolor(lon_ngT[:],lat_ngT[:],hfds_ng_B[:,:],cmap=cm_jet,shading='auto') 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='black',facecolor='none',linestyle='-',lw=1,hatch='\\')) plt.text(lonupA1,latupA2,'Zone upA', verticalalignment='top',horizontalalignment='left',color='black') # domain plt.text(lonlim1,latlim2-1.5,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonlim1,latlim2-3.5,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='top',horizontalalignment='left') # native_dxdy plt.text(-100,-10,'\\textbf{dx %.2f'%(dx), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='bottom',horizontalalignment='left') plt.text(-100,-10,'\\textbf{dy %.2f'%(dy), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='top',horizontalalignment='left') plt.text(-100,-15,'(at '+lon_dxdy_str+' '+lat_dxdy_str+')', fontsize=12,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),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 ... exp1_file is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1_file is used for title (cmip5_IPSL-CM5A-LR_thetao...) because we plot the NATIVE GRID depth_verif_str_ng=" -%.0fm"%(deptht_ng[0]) title23='\\large \\textbf{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{Flux Net Down ('+ng_Flux_str+')} (W/m2) \\textbf{'+depth_verif_str_ng+'} (T point)' plt.title(title23) # save img if debug : plt.show() else : img_name=exp1_name+'_'+domain_MAJ+'_'+ng_Flux_str+'_native_grid.'+img plt.savefig(img_dir_native+'/'+img_name) if img == 'png' : img_dir=img_dir_native with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(23) #--------------------- # Plot FLUX Zone upB along shore #--------------------- plt.figure(24) v=np.linspace(HFDS_max,HFDS_min,int((HFDS_max-HFDS_min)/20+1), endpoint=True) plt.pcolor(lon_ngT[:],lat_ngT[:],hfds_ng_B[:,:],cmap=cm_jet,edgecolors='blue', linewidth=0.3, linestyle='solid',shading='auto') 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-0.3,'\\Large \\textbf{'+domain_MAJ+'}',color='black',verticalalignment='top',horizontalalignment='left') # grid plt.text(lonupB1,latupB2-0.6,'\\Large \\textbf{NATIVE grid}',color='black',verticalalignment='top',horizontalalignment='left') # native_dxdy plt.text(-88.7,-10,'\\textbf{dx %.2f'%(dx), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='bottom',horizontalalignment='left') plt.text(-88.7,-10,'\\textbf{dy %.2f'%(dy), fontsize=24,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),color='green',verticalalignment='top',horizontalalignment='left') plt.text(-88.7,-11,'(at '+lon_dxdy_str+' '+lat_dxdy_str+')', fontsize=12,bbox=dict(facecolor='white', edgecolor='white',alpha=0.5),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 ... exp1_file is used for img savname (cmip5_IPSL-CM5A-LR_tr025_cdo_cdo..). Because it is convenient for the HTML table # exp1_file 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 ('+ng_Flux_str+')} (W/m2) \\textbf{'+depth_verif_str_ng+'} (T point)' plt.title(title24) # save img if debug : plt.show() else : img_name=exp1_name+'_'+domain_MAJ+'_'+ng_Flux_str+'_native_grid_ZOOM.'+img plt.savefig(img_dir_native+'/'+img_name) if img == 'png' : img_dir=img_dir_native with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(24) else : print('\n\nSorry the file we\'re looking for doesn\'t exist : \n', dir1_ng_hfds+'/'+fH_ng,'\n') # sys.exit(0) #=============================================================================== # CHRIS ISO on trop3D grid (TT1, UU1) #=============================================================================== if B_isodepth | B_section_WBDY : #------------------------------------------------------------------------------- # TT1: read thetao T3D file from domain file (peru, senegal...) #------------------------------------------------------------------------------- fT1='../'+domain_MAJ+'/thetao/'+exp1_file.replace('.nc', '.'+domain_MAJ+'.nc') print('\nflag B_isodepth | B_section_WBDY :') print(' fT1 :', fT1) print(' TT1 : '+domain_MAJ+' area : reading thetao...') ncT1=netcdf(fT1,'r') 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 1 glamt1 !~= lon1T' ) 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 1 gphit1 !~= lat1T' ) TT1=ncT1.variables['thetao'][:,:,:] sst=ncT1.variables['lev'][0] ncT1.close() #------------------------------------------------------------------------------- # w1 : weight from TT1 # (from masked TT1 => for all projects even if we don't have the mask) #------------------------------------------------------------------------------- mask_from_TT1=ma.copy(TT1[:,:,:]) mask_from_TT1[mask_from_TT1.mask==False]=1 w1=w1a*mask_from_TT1 w1_surf=w1a_surf*mask_from_TT1[0,:,:] if B_section_WBDY : w1_WBDY=w1[:,:,iWBDY] TT1_WBDY=TT1[:,:,iWBDY] # section WBDY #------------------------------------------------------------------------------- # UU1 : read U3D file from domain file (peru, senegal...) #------------------------------------------------------------------------------- print('\nflag B_section_WBDY :') fU1=fT1.replace('thetao', 'uo',2) print(' fU1 :', fU1) ncU1=netcdf(fU1,'r') print(' UU1 : '+domain_MAJ+' area : reading uo...\n') lon1U=ncU1.variables['nav_lon'][0,:] if 'cmip' in project : print('\n CHRIS WARNING 2 :in cmip5 files, gridU=gridT => glamu1 (grid U of trop mesh_mask) != lon1U (gridT of cmip5 U file) \n',glamu1-lon1U,'\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) UU1a=ncU1.variables['uo'][:,:,:] 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) #------------------------------------------------------------------------------- # maskupBas1_surf & 2 : mask for the ZOOM (at the surface) #------------------------------------------------------------------------------- ### Detect the first cell of land from East to West for upwelling B area TT1_sst_upB=TT1[0,jupB1:jupB2+un:,iupB1:iupB2+un] jdimupB,idimupB=TT1_sst_upB.shape maskupBall=np.zeros((jdimupB,idimupB)) maskupB1=np.zeros((jdimupB,idimupB)) maskupB2=np.zeros((jdimupB,idimupB)) maskupBall[TT1_sst_upB.mask]=1 upB_land_first_element=maskupBall[:,0:idimupB-1]+maskupBall[:,1:idimupB] for jj in range(jdimupB) : # find the first land element for latitude jj iit1=np.where(upB_land_first_element[jj,:]==1) iit=iit1[0][0] # print(k_iso,jj,iit) ii1=iit-upB_as_size1*coef+1 # find last land element for upBas1 (depend of resolution) ii2=iit-upB_as_size2*coef+1 # find last land element for upBas2 (depend of resolution) maskupB1[jj,ii1:iit+1]=1 # mask upBas1 (0 / 1) maskupB2[jj,ii2:ii1]=1 # mask upBas2 (0 / 1) maskupBas1_surf=ma.masked_values(maskupB1,0) # mask upBas1 (Nan / 1) maskupBas2_surf=ma.masked_values(maskupB2,0) # mask upBas1 (Nan / 1) #------------------------------------------------------------------------------- # TT1_iso : weighted temperature TT1 (for ISOLEV sst, 54m, 120m, 31-108m ...) #------------------------------------------------------------------------------- ii=0 wTT1_iso_mean_upA_list=[] wTT1_iso_mean_upB_list=[] wTT1_iso_mean_upBas1_list=[] # along shore wTT1_iso_mean_upBas2_list=[] # along shore wTT1a=w1[:,:,:]*TT1[:,:,:] for iik in range(len(knum)) : ik=knum[iik] k_iso=k_iso_list[ik] # if ('cmip' in project) & (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 : k1=14; k2=24 # pond upA (cumul sur z=k1:k2) w1_iso_upA=w1[k1:k2,jupA1:jupA2+un:,iupA1:iupA2+un] wTT1a_iso_upA=wTT1a[k1:k2,jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB (cumul sur z=k1:k2) w1_iso_upB=w1[k1:k2,jupB1:jupB2+un:,iupB1:iupB2+un] wTT1a_iso_upB=wTT1a[k1:k2,jupB1:jupB2+un:,iupB1:iupB2+un] # TT1_iso (cumul sur z=k1:k2) pour calcul TT1 mean des niveaux k1:k2 w1_iso=w1[k1:k2,:,:] wTT1a_iso=wTT1a[k1:k2,:,:] # compute TT1 for k1:k2 TT1_iso=wTT1a_iso.sum(0)/w1_iso.sum(0) depth_iso_str="-%.0f-%.0fm (k=%d-%d)"%(deptht1[k1],deptht1[k2],k1,k2) depth_iso_str_fig="%.0fm-%.0fm_k_%d-%d"%(deptht1[k1],deptht1[k2],k1,k2) elif k_iso == 2026 : k1=20; k2=26 # pond upA (cumul sur z=k1:k2) w1_iso_upA=w1[k1:k2,jupA1:jupA2+un:,iupA1:iupA2+un] wTT1a_iso_upA=wTT1a[k1:k2,jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB (cumul sur z=k1:k2) w1_iso_upB=w1[k1:k2,jupB1:jupB2+un:,iupB1:iupB2+un] wTT1a_iso_upB=wTT1a[k1:k2,jupB1:jupB2+un:,iupB1:iupB2+un] # TT1_iso (cumul sur z=k1:k2) pour calcul TT1 mean des niveaux k1:k2 w1_iso=w1[k1:k2,:,:] wTT1a_iso=wTT1a[k1:k2,:,:] # compute TT1 for k1:k2 TT1_iso=wTT1a_iso.sum(0)/w1_iso.sum(0) depth_iso_str="-%.0f-%.0fm (k=%d-%d)"%(deptht1[k1],deptht1[k2],k1,k2) depth_iso_str_fig="%.0fm-%.0fm_k_%d-%d"%(deptht1[k1],deptht1[k2],k1,k2) else : # pond upA w1_iso_upA=w1[k_iso,jupA1:jupA2+un:,iupA1:iupA2+un] wTT1a_iso_upA=wTT1a[k_iso,jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB w1_iso_upB=w1[k_iso,jupB1:jupB2+un:,iupB1:iupB2+un] wTT1a_iso_upB=wTT1a[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) wTT1_iso_mean_upA=wTT1a_iso_upA.sum()/w1_iso_upA.sum() wTT1_iso_mean_upB=wTT1a_iso_upB.sum()/w1_iso_upB.sum() wTT1_iso_mean_upA_list.append(wTT1_iso_mean_upA) wTT1_iso_mean_upB_list.append(wTT1_iso_mean_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] # pond upBas w1_iso_upBas1=w1_iso_upB*maskupBas1_surf wTT1a_iso_upBas1=wTT1a_iso_upB*maskupBas1_surf wTT1_iso_mean_upBas1=wTT1a_iso_upBas1.sum()/w1_iso_upBas1.sum() wTT1_iso_mean_upBas1_list.append(wTT1_iso_mean_upBas1) w1_iso_upBas2=w1_iso_upB*maskupBas2_surf wTT1a_iso_upBas2=wTT1a_iso_upB*maskupBas2_surf wTT1_iso_mean_upBas2=wTT1a_iso_upBas2.sum()/w1_iso_upBas2.sum() wTT1_iso_mean_upBas2_list.append(wTT1_iso_mean_upBas2) # TT1_iso TT1_iso_upBas1=TT1_iso_upB*maskupBas1_surf TT1_iso_upBas2=TT1_iso_upB*maskupBas2_surf # 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 # wTT1_iso_mean_upBas_list.append(TT1_iso_upBas.mean()) print("\n COMPUTE pond_mean for thetao level iik=%1d K_iso=%-5d depth=%-18s upA = %.1f upB = %.1f upBas = %.1f upBas2 = %.1f"%(iik,k_iso,depth_iso_str,wTT1_iso_mean_upA,wTT1_iso_mean_upB,wTT1_iso_mean_upBas1,wTT1_iso_mean_upBas2)) if B_isodepth_plot : ifig99=17 cm99=cm_jet lon99=glamt1[:] lat99=gphit1[:] var99=TT1_iso[:,:] var99_str='thetao' unit99='C' min99=Tmin_iso_list[ik] max99=Tmax_iso_list[ik] min_upBas99=Tmin_iso_upBas_list[ik] max_upBas99=Tmax_iso_upBas_list[ik] pond99_upA=wTT1_iso_mean_upA pond99_upB=wTT1_iso_mean_upB pond99_upBas1=wTT1_iso_mean_upBas1 pond99_upBas2=wTT1_iso_mean_upBas2 exp99_tex=exp1_tex img_name99=exp1_name+'_'+domain_MAJ+'_thetao_'+depth_iso_str_fig+'.'+img depth_iso_str99=depth_iso_str print(' flag B_isodepth_plot : plotting '+img_name99+'...') f_plot_iso(ifig99,cm99,lon99,lat99,var99,var99_str,unit99,min99,max99,min_upBas99,max_upBas99,pond99_upA,pond99_upB,pond99_upBas1,pond99_upBas2,exp99_tex,img_name99,depth_iso_str99,dx,dy) if B_hfdsO : #=============================================================================== # 2D SURF (hfdsO) #=============================================================================== # pond upA w1_surf_upA=w1_surf[jupA1:jupA2+un:,iupA1:iupA2+un] # pond upB w1_surf_upB=w1_surf[jupB1:jupB2+un:,iupB1:iupB2+un] # pond upBas w1_surf_upBas1=w1_surf_upB*maskupBas1_surf w1_surf_upBas2=w1_surf_upB*maskupBas2_surf dir1_hfdsO=dir1.replace('05_lev1_orig','03_interp-horiz').replace('thetao',domain_MAJ+'/hfdsO') fH1=exp1_file.replace('-lev1_orig.nc','.'+domain_MAJ+'.nc').replace('thetao','hfdsO') fH1_abs=dir1_hfdsO+'/'+fH1 print('\n dir1_hfdsO :', dir1_hfdsO) print( ' fH1_abs : ', fH1_abs) if os.path.exists(fH1_abs) : hfdsO_exists=True dsH1=xr.open_dataset(fH1_abs) hfdsO=dsH1.hfds # pond upA whfdsOa_upA=w1_surf_upA*hfdsO[jupA1:jupA2+un:,iupA1:iupA2+un] whfdsO_mean_upA=whfdsOa_upA.sum()/w1_surf_upA.sum() # pond upB whfdsOa_upB=w1_surf_upB*hfdsO[jupB1:jupB2+un:,iupB1:iupB2+un] whfdsO_mean_upB=whfdsOa_upB.sum()/w1_surf_upB.sum() # pond upBas whfdsOa_upBas1=whfdsOa_upB*maskupBas1_surf whfdsOa_upBas2=whfdsOa_upB*maskupBas2_surf whfdsO_mean_upBas1=whfdsOa_upBas1.sum()/w1_surf_upBas1.sum() whfdsO_mean_upBas2=whfdsOa_upBas2.sum()/w1_surf_upBas2.sum() print("\n COMPUTE pond_mean for hfdsO (W/m2) upA = %.1f upB = %.1f upBas = %.1f upBas2 = %.1f"%(whfdsO_mean_upA,whfdsO_mean_upB,whfdsO_mean_upBas1,whfdsO_mean_upBas2)) if B_hfdsO_plot : ifig99=37 cm99=cm_jet lon99=glamt1[:] lat99=gphit1[:] var99=hfdsO[:,:] 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=whfdsO_mean_upA pond99_upB=whfdsO_mean_upB pond99_upBas1=whfdsO_mean_upBas1 pond99_upBas2=whfdsO_mean_upBas2 exp99_tex=exp1_tex img_name99=exp1_name+'_'+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_upBas1,pond99_upBas2,exp99_tex,img_name99,depth_iso_str99,dx,dy) else : hfdsO_exists=False print('\n\nSorry the file we\'re looking for doesn\'t exist : \n', fH1_abs,'\n') # sys.exit(0) #=============================================================================== # 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 : #------------------------------------------------------------------------------- # 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 : print('\n',80*'-','\n flag B_section_WBDY : Reading WBDY files... \n',80*'-') #------------------------------------------------------------------------------- # UU1 TT1 already read in domain part #------------------------------------------------------------------------------- # Pour eviter les sections > WBDY: if config != 'peru12_nemo' : ##### if ('p12' in exp1_file) | ('erai7' in exp1_file) : gphit4_list.append(gphit1) deptht4_list.append(deptht1) UUU_WBDY_list.append(UU1_WBDY) TTT_WBDY_list.append(TT1_WBDY) pond_WBDY_list.append(w1_WBDY) lonWBDY_str_list.append(lonWBDY1_str) section_from_domain_list.append(section_from_domain1) print('section_'+lonWBDY1_str+' from '+section_from_domain1+' :') print(' fmask1, fT1, TT1, fU1, UU1 : read before (see above)') if B_section_TROP : print('\n',80*'-','\n flag B_section_TROP : Reading TROP files for each section... \n',80*'-') for section in section_list: #------------------------------------------------------------------------------- # mask2 : read TROP section mask #------------------------------------------------------------------------------- print('section_'+section+' from TROP:') mesh_mask2='mesh_mask_TROP'+resolution+'_section_'+section+'.nc' fmask2='/data/cholod/EUC_Upwelling/trop/trop'+resolution+'_nemo/inputs/'+mesh_mask2 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 #------------------------------------------------------------------------------- # TT2: read thetao T3D from TROP section #------------------------------------------------------------------------------- fT2='../section_'+section+'/thetao/'+exp1_file.replace('.nc','.section_'+section+'.nc') print(' fT2 :', fT2 ) print(' TT2 : reading thetao...') 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 3 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 3 gphit2 !~= lat2T \n') TT2a=ncT2.variables['thetao'][:,ilat1:ilat2+un,0] TT2_WBDY=TT2a*tmask2m ncT2.close() #------------------------------------------------------------------------------- # UU2: read uo U3D from TROP section #------------------------------------------------------------------------------- fU2='../section_'+section+'/uo/'+exp1_file.replace('.nc','.section_'+section+'.nc').replace('thetao', 'uo',2) print(' fU2 :', fU2 ) print(' UU2 : reading uo...') ncU2=netcdf(fU2,'r') lon2U=ncU2.variables['nav_lon'][0,:] if 'cmip' in project : print("\n CHRIS WARNING 4 :in cmip5 files, gridU=gridT => glamu2=%.2f (grid U of trop mesh_mask) != lon2U=%.2f (gridT of cmip5 U file) => diff=%.3f \n"%(glamu2[0],lon2U[0],glamu2[0]-lon2U[0])) else : if not( 'GLORYS2V3_ORCA025' in exp1_file) : 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 4 :in GLORYS files, gridU=gridT => glamu2=%.2f (grid U of trop mesh_mask) != lon2U=%.2f (gridT of GLORYS U file) => diff=%.3f \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) UU2a=ncU2.variables['uo'][:,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 #------------------------------------------------------------------------------- print('\n',80*'-','\n Find EUC center for TROP & BDY... \n',80*'-') 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] print('\nsection '+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,shading='auto') # plt.ylim(yf,y0) # plt.xlim(35,45) # plt.ylim(54,10) ax = plt.gca() ax.invert_yaxis() plt.grid() plt.title(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,shading='auto') # plt.ylim(yf,y0) # plt.xlim(35,45) # plt.ylim(54,10) ax = plt.gca() ax.invert_yaxis() plt.grid() # plt.title(exp1_tex+' '+lonWBDY_str) plt.title(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) 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,shading='auto') 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{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{uo} (m/s) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title11) # img_name=exp1_name+'_section_'+lonWBDY_str+'_from_'+section_from_domain+'_uo_0-300m.'+img if debug : plt.show() else : plt.savefig(img_dir_section+'/'+img_name) if img == 'png' : img_dir=img_dir_section with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(11) print('flag B_section_plot : plotting '+img_name+'...') #------------------------------------------------------------------------------- # 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,shading='auto') 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{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{uo} (m/s) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title111) img_name=exp1_name+'_section_'+lonWBDY_str+'_from_'+section_from_domain+'_uo_select_EUC_%.0f_percent_0-300m.'%(seuil_uo_EUC)+img if debug : plt.show() else : plt.savefig(img_dir_section+'/'+img_name) if img == 'png' : img_dir=img_dir_section with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(111) print('flag B_section_plot : plotting '+img_name+'...') #------------------------------------------------------------------------------- # 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,shading='auto') 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{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{uo} (m/s) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title1111) img_name=exp1_name+'_section_'+lonWBDY_str+'_from_'+section_from_domain+'_uo_select_EUC_%.0f_percent_75_levels.'%(seuil_uo_EUC)+img if debug : plt.show() else : plt.savefig(img_dir_section+'/'+img_name) if img == 'png' : img_dir=img_dir_section with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(1111) print('flag B_section_plot : plotting '+img_name+'...') #------------------------------------------------------------------------------- # 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,shading='auto') 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{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{thetao} (C) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title12) img_name=exp1_name+'_section_'+lonWBDY_str+'_from_'+section_from_domain+'_thetao_0-300m.'+img if debug : plt.show() else : plt.savefig(img_dir_section+'/'+img_name) if img == 'png' : img_dir=img_dir_section with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(12) print('flag B_section_plot : plotting '+img_name+'...') #------------------------------------------------------------------------------- # 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,shading='auto') 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{'+exp1_tex+ '} \\newline ('+str(yb1)+'-'+str(ye1)+') \\hspace{0.5cm} \\textbf{thetao} (C) \\hspace{0.5cm}section: \\textbf{'+lonWBDY_str+ '}' plt.title(title122) img_name=exp1_name+'_section_'+lonWBDY_str+'_from_'+section_from_domain+'_thetao_select_EUC_%.0f_percent_0-300m.'%(seuil_uo_EUC)+img if debug : plt.show() else : plt.savefig(img_dir_section+'/'+img_name) if img == 'png' : img_dir=img_dir_section with open(script_dir+'/EUC_Upwelling_characterize/modules/zip_script_png.py') as f3: exec(f3.read()) plt.close(122) print('flag B_section_plot : plotting '+img_name+'...') else : print('COMPUTE section '+lonWBDY_str+' from '+section_from_domain+'...') #------------------------------------------------------------------------------- # Fill the mean temperature & flux ascii file (out_area) #------------------------------------------------------------------------------- text2area_2=' ' text2section_2=' ' 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) text2section_2=text2section_2+'%-12.1f'%(TTUUmax_list[ii]) else : # print(la Temperature ponderee sur le coeur de l'EUC ...) text2section_2=text2section_2+'%-12.1f'%(pond_mean_TTT_WBDY[ii]) if B_hfdsO : if hfdsO_exists : text2area_2=text2area_2+\ '%-7.1f'%(whfdsO_mean_upA)+\ '%-7.1f'%(whfdsO_mean_upB)+\ '%-7.1f'%(whfdsO_mean_upBas1)+\ '%-7.1f'%(whfdsO_mean_upBas2) else : text2area_2=text2area_2+"none none none none " if B_isodepth : for ll in range(len(wTT1_iso_mean_upA_list)) : text2area_2=text2area_2+\ '%-7.1f'%(wTT1_iso_mean_upA_list[ll])+\ '%-7.1f'%(wTT1_iso_mean_upB_list[ll])+\ '%-7.1f'%(wTT1_iso_mean_upBas1_list[ll])+\ '%-7.1f'%(wTT1_iso_mean_upBas2_list[ll]) if 'cmip5' in exp1_name : id='id5c'+str(ii5c); ii5c=ii5c+1 elif 'cmip6' in exp1_name : id='id6c'+str(ii6c); ii6c=ii6c+1 elif 'trop025_nemo' in exp1_name : id='id25n'+str(ii25n); ii25n=ii25n+1 elif 'trop025_now' in exp1_name : id='id25c'+str(ii25c); ii25c=ii25c+1 elif 'trop075_nemo' in exp1_name : id='id75n'+str(ii75n); ii75n=ii75n+1 elif 'trop075_now' in exp1_name : id='id75c'+str(ii75c); ii75c=ii75c+1 elif 'trop12_nemo' in exp1_name : id='id12n'+str(ii12n); ii12n=ii12n+1 elif 'GLORYS' in exp1_name : id='idglo'+str(iiglo); iiglo=iiglo+1 elif 'soda3' in exp1_name : id='idsod'+str(iisod); iisod=iisod+1 text2common='%-90s'%(exp1_name)+'%-6.2f'%(dx)+'%-6.2f'%(dy)+'%-6s'%(yb1)+'%-6s'%(ye1)+'%-6s'%(id)+'\n' text2area_2=text2area_2+text2common text2section_2=text2section_2+text2common if (not(B_hfdsO)) | (B_hfdsO & hfdsO_exists) : fout_area.write(text2area_2); fout_area.flush() fout_section.write(text2section_2); fout_section.flush() #------------------------------- 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 (out_area) #------------------------------------------------------------------------------- if (B_section_TROP | B_section_WBDY) | B_isodepth : text2common=\ "\n\n# resolution='%3s'\n"%(resolution) text2area_3=text2common text2section_3=text2common text2area_3=text2area_3+\ "# 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) if B_hfdsO: text2area_3=text2area_3+\ "# Heat Flux: hfdsO\n" text2section_3=text2section_3+\ "# 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]) if uo_max : text2section_3=text2section_3+\ "# Temperature des sections: valeur au point uo max \n" else : text2section_3=text2section_3+\ "# Temperature des sections: valeur moyenne du coeur de l'EUC: seuil %.0f"%(seuil_uo_EUC)+"% de uo max\n" fout_area.write(text2area_3) fout_section.write(text2section_3) fout_area.close() fout_section.close() #shutil.move(out_area,scriptname_root+'_output.txt') #shutil.copy2(script_dir+'/'+scriptname,img_dir_area) #shutil.copy2(script_dir+'/'+scriptname,img_dir_section)