Changeset 3833 for trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4/vmr_ch4.py

Timestamp:

Jul 7, 2025, 1:39:13 PM (31 hours ago)

Author:

afalco

Message:

Pluto: updated plots scripts.
Fixed some issues with reading XIOS, etc.
Included display_netcdf.py tool from Mars PCM.
AF

Location:

trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4

Files:

• . (added)
• vmr_ch4.py (copied) (copied from trunk/LMDZ.PLUTO/util/script_figures/movie_winds/tsurf_winds.py) (7 diffs)

trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4/vmr_ch4.py

@@ -10 +10 @@
 import datetime
 from mpl_toolkits.basemap import Basemap, shiftgrid
-from FV3_utils import *
+from FV3_utils import * # import name
+from input import * # import 'name'
 
 ############################
-# basename="../../Xhistins2072"
-filename1="../"+name+"_A.nc"
+filename1="../"+name+"_A.nc" # define name in input.py
 filename2="../"+name+".nc"
-filename3="../../phisinit.nc"
-var="tsurf" #variable
-phisinit="phisinit" #variable
+var="vmr_ch4"
+var2="phisinit" #variable
 vari="u"
 varj="v"
-tint=[30,32] #Time must be as written in the input file
+font=26
+tint=[30,31] #Time must be as written in the input file
 #tintstr=["03:00","09:00","15:00","21:00"] #Time must be as written in the input file
 
-font=26
+# altitude in km
+altitude=1 
 
 nc1=Dataset(filename1)
 nc2=Dataset(filename2)
-nc3=Dataset(filename3)
 
 lat=getvar(nc1,"latitude")
@@ -35 +35 @@
 tim=getvar(nc1,"Time")
 ############################
+
+def swinglon(myvar):
+   # changer les longitudes pour mettre TR au centre
+   vec=shape(myvar)
+   myvarbis=np.zeros(vec,dtype='f')
+# i lat : pas de changement
+# j lon :
+   for i in range(vec[0]):
+       for j in range(vec[1]):
+           if j < int(vec[1]/2.) :
+             myvarbis[i,j]=myvar[i,j+int(vec[1]/2.)]
+             # myvar2bis[i,j]=myvar2[i,j+int(vec[1]/2)]
+           else:
+             myvarbis[i,j]=myvar[i,j-int(vec[1]/2.)]
+             # myvar2bis[i,j]=myvar2[i,j-int(vec[1]/2)]
+   return myvarbis
 
 def getwinds(lon,lat,vecx,vecy):
@@ -44 +60 @@
           color='black'    # arrow color
           pivot='mid'      # arrow around middle of box. Alternative : tip
-          scale=100     # scale arrow : plus grand = fleche plus petite
+          scale=200     # scale arrow : plus grand = fleche plus petite
           width=0.002      # width arrow
-          linewidths=0.5   # epaisseur contour arrow
+          linewidths=0.1   # epaisseur contour arrow
           edgecolors='k'   # couleur contour arrow
 
@@ -86 +102 @@
 def getfigvar(i):
     pal=get_cmap(name="jet")
-    lev=np.linspace(37,51,15)
-    # newlon=lon+180
-    newlon=lon 
+    lev=np.linspace(0.2,0.8,50)
+    newlon=lon+180
     CF=mpl.contourf(newlon, lat, myvarbis,lev,cmap=pal,extend='both')
     yticks=[-90,-60,-30,0,30,60,90]
     xticks=[0,60,120,180,240,300,360]
-    cbar=mpl.colorbar(CF,shrink=1, format="%1.0f")
-    cbar.ax.set_title("Tsurf [K]",y=1.04,fontsize=font)
+    cbar=mpl.colorbar(CF,shrink=1, format="%1.2f")
+    cbar.ax.set_title("VMR CH4 [%]",y=1.04,fontsize=font)
 
     for t in cbar.ax.get_yticklabels():
         t.set_fontsize(font)
 
-    c=mpl.contour(newlon, lat, phisinit2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5)
+    c=mpl.contour(newlon, lat, myvar2bis, 10,levels=np.linspace(-4,4,8), colors = 'k', linewidths = 3.5)
     mpl.clabel(c, fmt='%2.1f',inline=1, colors='k', fontsize=23,inline_spacing=1)
     #mpl.title('Local Time at Sputnik Planum='+str(i*3)+'H00',fontsize=font)
@@ -105 +120 @@
     mpl.xticks(xticks,fontsize=font)
     mpl.yticks(yticks,fontsize=font)
-    getwinds(newlon,lat,u,v)
+    #getwinds(newlon,lat,u,v)
 
-# def getnumalt(choicealt,alt):
-#     numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
-#     numalt=numalt[0][0]
-#     return numalt
+def getnumalt(choicealt,alt):
+    numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
+    numalt=numalt[0][0]
+    return numalt
 
 #######################
-# numalt=getnumalt(30,alt)
-numalt=np.searchsorted(alt[...],30)
+numalt=getnumalt(altitude,alt)
 print(('numalt =',numalt,'altitude=',alt[numalt]))
 uini=getvar(nc1,vari,tint,tim)[:,numalt]
 vini=getvar(nc1,varj,tint,tim)[:,numalt]
 myvar=getvar(nc2,var,tint,tim)
-phisinit2=getvar(nc3,phisinit) # phisinitmyvar2=phisinit2/0.6169/1000.  # altitude km
+myvar2=getvar(nc2,var2) # phisinit
+myvar2=myvar2/0.6169/1000.  # altitude km
 nbfig=uini.shape[0]
 print(("nbfig=",nbfig))
-phisinit2bis=np.copy(phisinit2)
+myvar2bis=swinglon(myvar2)
 
 for i in range(nbfig):
@@ -129 +144 @@
    v2=vini[i,:,:]
    myv=myvar[i,:,:]
-   u=np.copy(u2)
-   v=np.copy(v2)
-   myvarbis=np.copy(myv)
-   print(i)
+   u=swinglon(u2)
+   v=swinglon(v2)
+   myvarbis=swinglon(myv)
+   print(i,"/",nbfig)
    getfigvar(i)
    mpl.savefig('mapwinds'+str('{0:03}'.format(i))+'.eps',dpi=200)
@@ -146 +161 @@
 #mpl.subplots_adjust(hspace = .1)
 
+
 #mpl.show()