- Timestamp:
- Jul 7, 2025, 1:39:13 PM (31 hours ago)
- Location:
- trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4
- Files:
-
- 1 added
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4/vmr_ch4.py
r3832 r3833 10 10 import datetime 11 11 from mpl_toolkits.basemap import Basemap, shiftgrid 12 from FV3_utils import * 12 from FV3_utils import * # import name 13 from input import * # import 'name' 13 14 14 15 ############################ 15 # basename="../../Xhistins2072" 16 filename1="../"+name+"_A.nc" 16 filename1="../"+name+"_A.nc" # define name in input.py 17 17 filename2="../"+name+".nc" 18 filename3="../../phisinit.nc" 19 var="tsurf" #variable 20 phisinit="phisinit" #variable 18 var="vmr_ch4" 19 var2="phisinit" #variable 21 20 vari="u" 22 21 varj="v" 23 tint=[30,32] #Time must be as written in the input file 22 font=26 23 tint=[30,31] #Time must be as written in the input file 24 24 #tintstr=["03:00","09:00","15:00","21:00"] #Time must be as written in the input file 25 25 26 font=26 26 # altitude in km 27 altitude=1 27 28 28 29 nc1=Dataset(filename1) 29 30 nc2=Dataset(filename2) 30 nc3=Dataset(filename3)31 31 32 32 lat=getvar(nc1,"latitude") … … 35 35 tim=getvar(nc1,"Time") 36 36 ############################ 37 38 def swinglon(myvar): 39 # changer les longitudes pour mettre TR au centre 40 vec=shape(myvar) 41 myvarbis=np.zeros(vec,dtype='f') 42 # i lat : pas de changement 43 # j lon : 44 for i in range(vec[0]): 45 for j in range(vec[1]): 46 if j < int(vec[1]/2.) : 47 myvarbis[i,j]=myvar[i,j+int(vec[1]/2.)] 48 # myvar2bis[i,j]=myvar2[i,j+int(vec[1]/2)] 49 else: 50 myvarbis[i,j]=myvar[i,j-int(vec[1]/2.)] 51 # myvar2bis[i,j]=myvar2[i,j-int(vec[1]/2)] 52 return myvarbis 37 53 38 54 def getwinds(lon,lat,vecx,vecy): … … 44 60 color='black' # arrow color 45 61 pivot='mid' # arrow around middle of box. Alternative : tip 46 scale= 100 # scale arrow : plus grand = fleche plus petite62 scale=200 # scale arrow : plus grand = fleche plus petite 47 63 width=0.002 # width arrow 48 linewidths=0. 5# epaisseur contour arrow64 linewidths=0.1 # epaisseur contour arrow 49 65 edgecolors='k' # couleur contour arrow 50 66 … … 86 102 def getfigvar(i): 87 103 pal=get_cmap(name="jet") 88 lev=np.linspace(37,51,15) 89 # newlon=lon+180 90 newlon=lon 104 lev=np.linspace(0.2,0.8,50) 105 newlon=lon+180 91 106 CF=mpl.contourf(newlon, lat, myvarbis,lev,cmap=pal,extend='both') 92 107 yticks=[-90,-60,-30,0,30,60,90] 93 108 xticks=[0,60,120,180,240,300,360] 94 cbar=mpl.colorbar(CF,shrink=1, format="%1. 0f")95 cbar.ax.set_title(" Tsurf [K]",y=1.04,fontsize=font)109 cbar=mpl.colorbar(CF,shrink=1, format="%1.2f") 110 cbar.ax.set_title("VMR CH4 [%]",y=1.04,fontsize=font) 96 111 97 112 for t in cbar.ax.get_yticklabels(): 98 113 t.set_fontsize(font) 99 114 100 c=mpl.contour(newlon, lat, phisinit2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5)115 c=mpl.contour(newlon, lat, myvar2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5) 101 116 mpl.clabel(c, fmt='%2.1f',inline=1, colors='k', fontsize=23,inline_spacing=1) 102 117 #mpl.title('Local Time at Sputnik Planum='+str(i*3)+'H00',fontsize=font) … … 105 120 mpl.xticks(xticks,fontsize=font) 106 121 mpl.yticks(yticks,fontsize=font) 107 getwinds(newlon,lat,u,v)122 #getwinds(newlon,lat,u,v) 108 123 109 #def getnumalt(choicealt,alt):110 #numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))111 #numalt=numalt[0][0]112 #return numalt124 def getnumalt(choicealt,alt): 125 numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt))) 126 numalt=numalt[0][0] 127 return numalt 113 128 114 129 ####################### 115 # numalt=getnumalt(30,alt) 116 numalt=np.searchsorted(alt[...],30) 130 numalt=getnumalt(altitude,alt) 117 131 print(('numalt =',numalt,'altitude=',alt[numalt])) 118 132 uini=getvar(nc1,vari,tint,tim)[:,numalt] 119 133 vini=getvar(nc1,varj,tint,tim)[:,numalt] 120 134 myvar=getvar(nc2,var,tint,tim) 121 phisinit2=getvar(nc3,phisinit) # phisinitmyvar2=phisinit2/0.6169/1000. # altitude km 135 myvar2=getvar(nc2,var2) # phisinit 136 myvar2=myvar2/0.6169/1000. # altitude km 122 137 nbfig=uini.shape[0] 123 138 print(("nbfig=",nbfig)) 124 phisinit2bis=np.copy(phisinit2)139 myvar2bis=swinglon(myvar2) 125 140 126 141 for i in range(nbfig): … … 129 144 v2=vini[i,:,:] 130 145 myv=myvar[i,:,:] 131 u= np.copy(u2)132 v= np.copy(v2)133 myvarbis= np.copy(myv)134 print(i )146 u=swinglon(u2) 147 v=swinglon(v2) 148 myvarbis=swinglon(myv) 149 print(i,"/",nbfig) 135 150 getfigvar(i) 136 151 mpl.savefig('mapwinds'+str('{0:03}'.format(i))+'.eps',dpi=200) … … 146 161 #mpl.subplots_adjust(hspace = .1) 147 162 163 148 164 #mpl.show() 149 165
Note: See TracChangeset
for help on using the changeset viewer.