source: trunk/LMDZ.PLUTO/util/script_figures/proftempNH_SP.py @ 3823

#! /usr/bin/env python
from ppclass import pp
from    netCDF4               import    Dataset
from    numpy                 import    *
import  numpy                 as        np
import  matplotlib.pyplot     as        mpl
from matplotlib.cm import get_cmap
import pylab
from matplotlib import ticker
import matplotlib.colors as colors
import datetime
from mpl_toolkits.basemap import Basemap, shiftgrid

############################
filename1="restart_ref/diagfi2015_A.nc"
var="temperature" #variable

nc1=Dataset(filename1)
lat=nc1.variables["lat"][:]
lon=nc1.variables["lon"][:]
alt=nc1.variables["altitude"][:]
tim=nc1.variables["time_counter"][:]

print(('Time = ',tim))
print(('Long = ',lon))
print(('Lat = ',lat))


############################

def findindextime(tini,lt0,lt1,p):
    lt180=(lt0+12)%24   # a t=0 et longitude=180
    lt0p1=lt180*((p[0]+360)%360)/180   # local time a t=0 a la longitude p1
    diff=lt1-lt0p1                 # diff du local time a p1 a t=0 avec celui recherche
    indp1=tini+diff/24.                 # on adapte lindice pour tomber sur le bon moment

    # on cherche dans Time l'indice le plus proche :
    indt1=np.where(abs(tim[:]-indp1)==min(abs(tim[:]-indp1)))[0][0]
    print(('Point =',p,' Time=',tim[indt1]))
    return indt1
 
def getvar(filename,var):
    myvar = pp(file=filename,var=var,compute="nothing").getf()
    return myvar

def getindex(lat,lon,mylat,mylon):
    indlat=np.where(abs(lat[:]-mylat)==min(abs(lat[:]-mylat)))[0][0]
    indlon=np.where(abs(lon[:]-mylon)==min(abs(lon[:]-mylon)))[0][0]
    print((lon[indlon],lat[indlat]))
    return indlat,indlon

def main(tini,lt0,ltchoice,p):
    indt=findindextime(tini,lt0,ltchoice,p)
    indlat,indlon=getindex(lat,lon,p[1],p[0])
    myvar=getvar(filename1,var)[indt,:,indlat,indlon]
    return myvar
############################

#points: entre -180 et 180 
p1=[180, 25] # Entry
p2=[180, 20]    # low SP
p3=[180, 15]    # Burney
p4=[180, 10]    # SP
p5=[180, 5]    # BTD
p6=[180, 0]    # TR
p7=[180, -5]    # Lowell
p8=[180, -10]    # south pole

tini=32 # choix du jour dans le diagfi
lt0=0   # a t=0 et longitude=0
ltchoice=12   # local time choisi

myvar1=main(tini,lt0,ltchoice,p1)
#myvar2=main(tini,lt0,ltchoice,p2)
myvar3=main(tini,lt0,ltchoice,p3)
#myvar4=main(tini,lt0,ltchoice,p4)
myvar5=main(tini,lt0,ltchoice,p5)
myvar6=main(tini,lt0,ltchoice,p6)
myvar7=main(tini,lt0,ltchoice,p7)
#myvar8=main(tini,lt0,ltchoice,p8)


mpl.figure(figsize=(10, 10))

font=23

#xticks=[-90,-60,-30,0,30,60,90]
#yticks=np.linspace(0,240,9)
alt=alt/1000.

mpl.plot(myvar1,alt,'r',label=r'lat=25 $^\circ$')
mpl.plot(myvar3,alt,'k',label=r'lat=15 $^\circ$')
mpl.plot(myvar5,alt,'c',label=r'lat=5 $^\circ$')
mpl.plot(myvar6,alt,'b',label=r'lat= 0 $^\circ$')
mpl.plot(myvar7,alt,'y',label=r'lat=-5 $^\circ$')

mpl.legend(prop={'size':20},loc='upper left')
#mpl.title('Latitude ='+str(tintstr[i]),fontsize=font)
mpl.ylabel('Altitude (km)',labelpad=10,fontsize=font)
mpl.xlabel('Temperature (K)',labelpad=10, fontsize=font)
#mpl.xticks(xticks,fontsize=font)
mpl.xticks(fontsize=font)
#mpl.yticks(yticks,fontsize=font)
mpl.yticks(fontsize=font)
mpl.grid()
#mpl.legend(["Ref","Alt"])

pylab.ylim([-3,0])
pylab.xlim([35,45])

mpl.savefig('proftempNH_SP180_'+str(ltchoice)+'h.eps',dpi=200)
mpl.savefig('proftempNH_SP180_'+str(ltchoice)+'h.png',dpi=200)
mpl.show()