Ignore:
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:
1 added
1 copied

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.PLUTO/util/script_figures/movie_vmr_ch4/vmr_ch4.py

    r3832 r3833  
    1010import datetime
    1111from mpl_toolkits.basemap import Basemap, shiftgrid
    12 from FV3_utils import *
     12from FV3_utils import * # import name
     13from input import * # import 'name'
    1314
    1415############################
    15 # basename="../../Xhistins2072"
    16 filename1="../"+name+"_A.nc"
     16filename1="../"+name+"_A.nc" # define name in input.py
    1717filename2="../"+name+".nc"
    18 filename3="../../phisinit.nc"
    19 var="tsurf" #variable
    20 phisinit="phisinit" #variable
     18var="vmr_ch4"
     19var2="phisinit" #variable
    2120vari="u"
    2221varj="v"
    23 tint=[30,32] #Time must be as written in the input file
     22font=26
     23tint=[30,31] #Time must be as written in the input file
    2424#tintstr=["03:00","09:00","15:00","21:00"] #Time must be as written in the input file
    2525
    26 font=26
     26# altitude in km
     27altitude=1
    2728
    2829nc1=Dataset(filename1)
    2930nc2=Dataset(filename2)
    30 nc3=Dataset(filename3)
    3131
    3232lat=getvar(nc1,"latitude")
     
    3535tim=getvar(nc1,"Time")
    3636############################
     37
     38def 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
    3753
    3854def getwinds(lon,lat,vecx,vecy):
     
    4460          color='black'    # arrow color
    4561          pivot='mid'      # arrow around middle of box. Alternative : tip
    46           scale=100     # scale arrow : plus grand = fleche plus petite
     62          scale=200     # scale arrow : plus grand = fleche plus petite
    4763          width=0.002      # width arrow
    48           linewidths=0.5   # epaisseur contour arrow
     64          linewidths=0.1   # epaisseur contour arrow
    4965          edgecolors='k'   # couleur contour arrow
    5066
     
    86102def getfigvar(i):
    87103    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
    91106    CF=mpl.contourf(newlon, lat, myvarbis,lev,cmap=pal,extend='both')
    92107    yticks=[-90,-60,-30,0,30,60,90]
    93108    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)
    96111
    97112    for t in cbar.ax.get_yticklabels():
    98113        t.set_fontsize(font)
    99114
    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)
    101116    mpl.clabel(c, fmt='%2.1f',inline=1, colors='k', fontsize=23,inline_spacing=1)
    102117    #mpl.title('Local Time at Sputnik Planum='+str(i*3)+'H00',fontsize=font)
     
    105120    mpl.xticks(xticks,fontsize=font)
    106121    mpl.yticks(yticks,fontsize=font)
    107     getwinds(newlon,lat,u,v)
     122    #getwinds(newlon,lat,u,v)
    108123
    109 # def getnumalt(choicealt,alt):
    110 #     numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
    111 #     numalt=numalt[0][0]
    112 #     return numalt
     124def getnumalt(choicealt,alt):
     125    numalt=np.where(abs(alt-choicealt)==min(abs(alt-choicealt)))
     126    numalt=numalt[0][0]
     127    return numalt
    113128
    114129#######################
    115 # numalt=getnumalt(30,alt)
    116 numalt=np.searchsorted(alt[...],30)
     130numalt=getnumalt(altitude,alt)
    117131print(('numalt =',numalt,'altitude=',alt[numalt]))
    118132uini=getvar(nc1,vari,tint,tim)[:,numalt]
    119133vini=getvar(nc1,varj,tint,tim)[:,numalt]
    120134myvar=getvar(nc2,var,tint,tim)
    121 phisinit2=getvar(nc3,phisinit) # phisinitmyvar2=phisinit2/0.6169/1000.  # altitude km
     135myvar2=getvar(nc2,var2) # phisinit
     136myvar2=myvar2/0.6169/1000.  # altitude km
    122137nbfig=uini.shape[0]
    123138print(("nbfig=",nbfig))
    124 phisinit2bis=np.copy(phisinit2)
     139myvar2bis=swinglon(myvar2)
    125140
    126141for i in range(nbfig):
     
    129144   v2=vini[i,:,:]
    130145   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)
    135150   getfigvar(i)
    136151   mpl.savefig('mapwinds'+str('{0:03}'.format(i))+'.eps',dpi=200)
     
    146161#mpl.subplots_adjust(hspace = .1)
    147162
     163
    148164#mpl.show()
    149165
Note: See TracChangeset for help on using the changeset viewer.