Index: trunk/MESOSCALE/LMD_MM_MARS/SIMU/runmeso =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SIMU/runmeso (revision 258) +++ trunk/MESOSCALE/LMD_MM_MARS/SIMU/runmeso (revision 261) @@ -22,10 +22,18 @@ ########################################### -changeregis=0 ## 1 if changed registry -#changeregis=1 ## 1 if changed registry + i_want_to_compile=1 - scenario="" #scenario="storm" + +########################################### +########################################### + +changeregis=0 +while getopts "r" options; do + case $options in + r ) changeregis=1;; + esac +done if [[ "${MESO}" = "" ]] Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/namelist.api =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/namelist.api (revision 258) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/namelist.api (revision 261) @@ -1,12 +1,22 @@ &io - path_to_input = './' - path_to_output = './' - input_name = 'wrfout_d01_9999-09-09_09:00:00' - process = 'list' !! list fields required in "fields" (available tk, tpot, GHT) - fields = 'tk,W,uvmet' + path_to_input = './' !! where input wrfout* files are located + path_to_output = './' !! where output API files will be located + input_name = 'wrfout_d01_9999-09-09_09:00:00' !! input file to API (could be wrfout*) + process = 'list' !! [do not modify] + fields = 'tk,W,uvmet' !! a list of fields to interpolate + !! - either fields in wrfout* + !! - or tk for temperature + uvmet for meteorological winds + tpot for potential temperature + debug = .TRUE. !! [add this if you want more information on screen] / &interp_in - interp_method = 4 - interp_levels = 0.050 + interp_method = 4 !! 1 --> INTERPOLATION: PRESSURE [LINEAR in p] output: wrfout*_p + !! 2 --> INTERPOLATION: PRESSURE [LINEAR in log(p)] output: wrfout*_p + !! 3 --> INTERPOLATION: ALTITUDE ABOVE MOLA AREOID output: wrfout*_z + !! 4 --> INTERPOLATION: ALTITUDE ABOVE LOCAL SURFACE output: wrfout*_zabg + interp_levels = 0.050 !! Interpolation levels: - pressure in hPa for interp_method = 1 or 2 + !! - altitude in km for interp_method = 3 or 4 + !! - [pressure shall be in decreasing order] / Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api.F90 =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api.F90 (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api.F90 (revision 261) @@ -0,0 +1,1 @@ +link ../../../../MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90 Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_g95.sh =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_g95.sh (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_g95.sh (revision 261) @@ -0,0 +1,11 @@ +#! /bin/bash + +#NETCDF=/distrib/local/netcdf-4.0.1/g95_4.0.3_64/ +#echo $NETCDF + +touch logc +touch loge +\rm logc +\rm loge +f2py -c -m api api.F90 --fcompiler=g95 -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include --f90flags="-Wall -Wno=112,141,137,155 -fno-second-underscore -ffree-form" > logc 2> loge +f2py -c -m timestuff time.F --fcompiler=g95 >> logc 2>> loge Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_ifort.sh =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_ifort.sh (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_ifort.sh (revision 261) @@ -0,0 +1,11 @@ +#! /bin/bash + + +echo $NETCDF + +touch logc +touch loge +\rm logc +\rm loge +f2py -c -m api api.F90 --fcompiler=intelem -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include > logc 2> loge +f2py -c -m timestuff time.F --fcompiler=intelem >> logc 2>> loge Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_pgf90.sh =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_pgf90.sh (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_pgf90.sh (revision 261) @@ -0,0 +1,16 @@ +#! /bin/bash + + +#NETCDF=/distrib/local/netcdf-4.0.1/pgi_7.1-6_32/ +#NETCDF=/distrib/local/netcdf-4.0.1/pgi_7.1-6_64/ +#NETCDF=/distrib/local/netcdf-3.6.0-p1/pgi_64bits/ +NETCDF=/donnees/aslmd/MODELES/MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic-3.6.1/ + +echo $NETCDF + +touch logc +touch loge +\rm logc +\rm loge +f2py -c -m api api.F90 --fcompiler=pg -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include --f90flags="-Mfree" > logc 2> loge +f2py -c -m timestuff time.F --fcompiler=pg >> logc 2>> loge Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_wrapper.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_wrapper.py (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/api_wrapper.py (revision 261) @@ -0,0 +1,35 @@ + +### A. Spiga -- LMD -- 03/07/2011 + +def api_onelevel ( path_to_input = './', \ + input_name = 'wrfout_d0?_????-??-??_??:00:00', \ + path_to_output = None, \ + output_name = None, \ + process = 'list', \ + fields = 'tk,W,uvmet,HGT', \ + debug = False, \ + bit64 = False, \ + oldvar = True, \ + interp_method = 4, \ + extrapolate = 0, \ + unstagger_grid = False, \ + onelevel = 0.020, \ + nocall = False ): + import api + import numpy as np + + if not path_to_output: path_to_output = path_to_input + + if not output_name: + if interp_method <= 2: output_name = input_name+'_p' + if interp_method == 3: output_name = input_name+'_z' + if interp_method == 4: output_name = input_name+'_zabg' + + #print input_name, output_name + + if nocall: pass + else: api.api_main ( path_to_input, input_name, path_to_output, output_name, \ + process, fields, debug, bit64, oldvar, np.arange(299), \ + interp_method, extrapolate, unstagger_grid, onelevel ) + + return output_name Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/domain.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/domain.py (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/domain.py (revision 261) @@ -0,0 +1,67 @@ +#!/usr/bin/env python + +### A. Spiga -- LMD -- 30/06/2011 -- slight modif early 07/2011 + +def domain (namefile,proj=None,back="vishires",target=None): + from netCDF4 import Dataset + from myplot import getcoord2d,define_proj,makeplotres,simplinterv,getprefix,dumpbdy,getproj,latinterv,wrfinterv,simplinterv + from mymath import max,min + from matplotlib.pyplot import contourf,rcParams,pcolor + from numpy.core.defchararray import find + from numpy import arange + ### + nc = Dataset(namefile) + ### + if proj == None: proj = "ortho" #proj = getproj(nc) + ### + prefix = namefile[0] + namefile[1] + namefile[2] + if prefix == "geo": + [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') + var = 'HGT_M' + zeplot = "domain" + else: + [lon2d,lat2d] = getcoord2d(nc) + var = "HGT" + zeplot = getprefix(nc) + "domain" + ### + lon2d = dumpbdy(lon2d,5) + lat2d = dumpbdy(lat2d,5) + if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") + elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) + else: [wlon,wlat] = simplinterv(lon2d,lat2d) + ### + m = define_proj(proj,wlon,wlat,back=back) + x, y = m(lon2d, lat2d) + ### + what_I_plot = dumpbdy(nc.variables[var][0,:,:], 5) + #levinterv = 250. + #zelevels = arange(min(what_I_plot)-levinterv,max(what_I_plot)+levinterv,levinterv) + zelevels = 30 + contourf(x, y, what_I_plot, zelevels) + #pcolor(x,y,what_I_plot) ## on voit trop les lignes ! + ### + if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot + else: zeplot = target + "/" + zeplot + ### + pad_inches_value = 0.35 + makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature + makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False) + #makeplotpng(zeplot,pad_inches_value=0.35) + #rcParams['savefig.facecolor'] = 'black' + #makeplotpng(zeplot+"b",pad_inches_value=0.35) + +if __name__ == "__main__": + import sys + ### to be replaced by argparse + from optparse import OptionParser + parser = OptionParser() + parser.add_option('-f', action='store', dest='namefile', type="string", default=None, help='name of WRF file [NEEDED]') + parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection') + parser.add_option('-b', action='store', dest='back', type="string", default="vishires", help='background') + parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder') + (opt,args) = parser.parse_args() + if opt.namefile is None: + print "I want to eat one file at least ! Use domain.py -f name_of_my_file. Or type domain.py -h" + exit() + print "Options:", opt + domain (opt.namefile,proj=opt.proj,back=opt.back,target=opt.target) Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/mymath.py (revision 261) @@ -0,0 +1,15 @@ +def min (field,axis=None): + import numpy as np + return np.array(field).min(axis=axis) + +def max (field,axis=None): + import numpy as np + return np.array(field).max(axis=axis) + +def mean (field,axis=None): + import numpy as np + return np.array(field).mean(axis=axis) + +def deg (): + return u'\u00b0' + Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/myplot.py (revision 261) @@ -0,0 +1,589 @@ +def errormess(text,printvar=None): + print text + if printvar: print printvar + exit() + return + +def getname(var=False,winds=False,anomaly=False): + if var and winds: basename = var + '_UV' + elif var: basename = var + elif winds: basename = 'UV' + else: errormess("please set at least winds or var",printvar=nc.variables) + if anomaly: basename = 'd' + basename + return basename + +def localtime(utc,lon): + ltst = utc + lon / 15. + ltst = int (ltst * 10) / 10. + ltst = ltst % 24 + return ltst + +def whatkindfile (nc): + if 'controle' in nc.variables: typefile = 'gcm' + elif 'phisinit' in nc.variables: typefile = 'gcmex' + elif 'vert' in nc.variables: typefile = 'mesoapi' + elif 'U' in nc.variables: typefile = 'meso' + elif 'HGT_M' in nc.variables: typefile = 'geo' + else: errormess("whatkindfile: typefile not supported.") + return typefile + +def getfield (nc,var): + ## this allows to get much faster (than simply referring to nc.variables[var]) + dimension = len(nc.variables[var].dimensions) + if dimension == 2: field = nc.variables[var][:,:] + elif dimension == 3: field = nc.variables[var][:,:,:] + elif dimension == 4: field = nc.variables[var][:,:,:,:] + return field + +def reducefield (input,d4=None,d3=None,d2=None,d1=None): + ### we do it the reverse way to be compliant with netcdf "t z y x" or "t y x" or "y x" + ### it would be actually better to name d4 d3 d2 d1 as t z y x + import numpy as np + dimension = np.array(input).ndim + shape = np.array(input).shape + print 'dim,shape: ',dimension,shape + output = input + error = False + if dimension == 2: + if d2 >= shape[0]: error = True + elif d1 >= shape[1]: error = True + elif d1 is not None and d2 is not None: output = input[d2,d1] + elif d1 is not None: output = input[:,d1] + elif d2 is not None: output = input[d2,:] + elif dimension == 3: + if d4 >= shape[0]: error = True + elif d2 >= shape[1]: error = True + elif d1 >= shape[2]: error = True + elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,d2,d1] + elif d4 is not None and d2 is not None: output = input[d4,d2,:] + elif d4 is not None and d1 is not None: output = input[d4,:,d1] + elif d2 is not None and d1 is not None: output = input[:,d2,d1] + elif d1 is not None: output = input[:,:,d1] + elif d2 is not None: output = input[:,d2,:] + elif d4 is not None: output = input[d4,:,:] + elif dimension == 4: + if d4 >= shape[0]: error = True + elif d3 >= shape[1]: error = True + elif d2 >= shape[2]: error = True + elif d1 >= shape[3]: error = True + elif d4 is not None and d3 is not None and d2 is not None and d1 is not None: output = input[d4,d3,d2,d1] + elif d4 is not None and d3 is not None and d2 is not None: output = input[d4,d3,d2,:] + elif d4 is not None and d3 is not None and d1 is not None: output = input[d4,d3,:,d1] + elif d4 is not None and d2 is not None and d1 is not None: output = input[d4,:,d2,d1] + elif d3 is not None and d2 is not None and d1 is not None: output = input[:,d3,d2,d1] + elif d4 is not None and d3 is not None: output = input[d4,d3,:,:] + elif d4 is not None and d2 is not None: output = input[d4,:,d2,:] + elif d4 is not None and d1 is not None: output = input[d4,:,:,d1] + elif d3 is not None and d2 is not None: output = input[:,d3,d2,:] + elif d3 is not None and d1 is not None: output = input[:,d3,:,d1] + elif d2 is not None and d1 is not None: output = input[:,:,d2,d1] + elif d1 is not None: output = input[:,:,:,d1] + elif d2 is not None: output = input[:,:,d2,:] + elif d3 is not None: output = input[:,d3,:,:] + elif d4 is not None: output = input[d4,:,:,:] + dimension = np.array(output).ndim + shape = np.array(output).shape + print 'dim,shape: ',dimension,shape + return output, error + +def definesubplot ( numplot, fig ): + from matplotlib.pyplot import rcParams + rcParams['font.size'] = 12. ## default (important for multiple calls) + if numplot == 4: + sub = 221 + fig.subplots_adjust(wspace = 0.3, hspace = 0.3) + rcParams['font.size'] = int( rcParams['font.size'] * 2. / 3. ) + elif numplot == 2: + sub = 121 + fig.subplots_adjust(wspace = 0.35) + rcParams['font.size'] = int( rcParams['font.size'] * 3. / 4. ) + elif numplot == 3: + sub = 131 + fig.subplots_adjust(wspace = 0.5) + rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) + elif numplot == 6: + sub = 231 + fig.subplots_adjust(wspace = 0.4, hspace = 0.0) + rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) + elif numplot == 8: + sub = 331 #241 + fig.subplots_adjust(wspace = 0.3, hspace = 0.3) + rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) + elif numplot == 9: + sub = 331 + fig.subplots_adjust(wspace = 0.3, hspace = 0.3) + rcParams['font.size'] = int( rcParams['font.size'] * 1. / 2. ) + elif numplot == 1: + sub = 99999 + elif numplot <= 0: + sub = 99999 + else: + print "supported: 1,2,3,4,6,8,9" + exit() + return sub + +def getstralt(nc,nvert): + typefile = whatkindfile(nc) + if typefile is 'meso': + stralt = "_lvl" + str(nvert) + elif typefile is 'mesoapi': + zelevel = int(nc.variables['vert'][nvert]) + if abs(zelevel) < 10000.: strheight=str(zelevel)+"m" + else: strheight=str(int(zelevel/1000.))+"km" + if 'altitude' in nc.dimensions: stralt = "_"+strheight+"-AMR" + elif 'altitude_abg' in nc.dimensions: stralt = "_"+strheight+"-ALS" + elif 'bottom_top' in nc.dimensions: stralt = "_"+strheight + elif 'pressure' in nc.dimensions: stralt = "_"+str(zelevel)+"Pa" + else: stralt = "" + else: + stralt = "" + return stralt + +def getlschar ( namefile ): + from netCDF4 import Dataset + from timestuff import sol2ls + from numpy import array + nc = Dataset(namefile) + zetime = None + if 'Times' in nc.variables: + zetime = nc.variables['Times'][0] + shape = array(nc.variables['Times']).shape + if shape[0] < 2: zetime = None + if zetime is not None \ + and 'vert' not in nc.variables: + #### strangely enough this does not work for api or ncrcat results! + zetimestart = getattr(nc, 'START_DATE') + zeday = int(zetime[8]+zetime[9]) - int(zetimestart[8]+zetimestart[9]) + if zeday < 0: lschar="" ## might have crossed a month... fix soon + else: lschar="_Ls"+str( int( 10. * sol2ls ( getattr( nc, 'JULDAY' ) + zeday ) ) / 10. ) + ### + zetime2 = nc.variables['Times'][1] + one = int(zetime[11]+zetime[12]) + int(zetime[14]+zetime[15])/37. + next = int(zetime2[11]+zetime2[12]) + int(zetime2[14]+zetime2[15])/37. + zehour = one + zehourin = abs ( next - one ) + else: + lschar="" + zehour = 0 + zehourin = 1 + return lschar, zehour, zehourin + +def getprefix (nc): + prefix = 'LMD_MMM_' + prefix = prefix + 'd'+str(getattr(nc,'GRID_ID'))+'_' + prefix = prefix + str(int(getattr(nc,'DX')/1000.))+'km_' + return prefix + +def getproj (nc): + typefile = whatkindfile(nc) + if typefile in ['mesoapi','meso','geo']: + ### (il faudrait passer CEN_LON dans la projection ?) + map_proj = getattr(nc, 'MAP_PROJ') + cen_lat = getattr(nc, 'CEN_LAT') + if map_proj == 2: + if cen_lat > 10.: + proj="npstere" + print "NP stereographic polar domain" + else: + proj="spstere" + print "SP stereographic polar domain" + elif map_proj == 1: + print "lambert projection domain" + proj="lcc" + elif map_proj == 3: + print "mercator projection" + proj="merc" + else: + proj="merc" + elif typefile in ['gcm']: proj="cyl" ## pb avec les autres (de trace derriere la sphere ?) + else: proj="ortho" + return proj + +def ptitle (name): + from matplotlib.pyplot import title + title(name) + print name + +def polarinterv (lon2d,lat2d): + import numpy as np + wlon = [np.min(lon2d),np.max(lon2d)] + ind = np.array(lat2d).shape[0] / 2 ## to get a good boundlat and to get the pole + wlat = [np.min(lat2d[ind,:]),np.max(lat2d[ind,:])] + return [wlon,wlat] + +def simplinterv (lon2d,lat2d): + import numpy as np + return [[np.min(lon2d),np.max(lon2d)],[np.min(lat2d),np.max(lat2d)]] + +def wrfinterv (lon2d,lat2d): + nx = len(lon2d[0,:])-1 + ny = len(lon2d[:,0])-1 + lon1 = lon2d[0,0] + lon2 = lon2d[nx,ny] + lat1 = lat2d[0,0] + lat2 = lat2d[nx,ny] + if abs(0.5*(lat1+lat2)) > 60.: wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 + else: wider = 0. + if lon1 < lon2: wlon = [lon1, lon2 + wider] + else: wlon = [lon2, lon1 + wider] + if lat1 < lat2: wlat = [lat1, lat2] + else: wlat = [lat2, lat1] + return [wlon,wlat] + +def makeplotres (filename,res=None,pad_inches_value=0.25,folder='',disp=True,ext='png',erase=False): + import matplotlib.pyplot as plt + from os import system + addstr = "" + if res is not None: + res = int(res) + addstr = "_"+str(res) + name = filename+addstr+"."+ext + if folder != '': name = folder+'/'+name + plt.savefig(name,dpi=res,bbox_inches='tight',pad_inches=pad_inches_value) + if disp: display(name) + if ext in ['eps','ps','svg']: system("tar czvf "+name+".tar.gz "+name+" ; rm -f "+name) + if erase: system("mv "+name+" to_be_erased") + return + +def dumpbdy (field,n,stag=None): + nx = len(field[0,:])-1 + ny = len(field[:,0])-1 + if stag == 'U': nx = nx-1 + if stag == 'V': ny = ny-1 + return field[n:ny-n,n:nx-n] + +def getcoorddef ( nc ): + ## getcoord2d for predefined types + typefile = whatkindfile(nc) + if typefile in ['mesoapi','meso']: + [lon2d,lat2d] = getcoord2d(nc) + lon2d = dumpbdy(lon2d,6) + lat2d = dumpbdy(lat2d,6) + elif typefile in ['gcm','gcmex']: + [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) + elif typefile in ['geo']: + [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') + return lon2d,lat2d + +def getcoord2d (nc,nlat='XLAT',nlon='XLONG',is1d=False): + import numpy as np + if is1d: + lat = nc.variables[nlat][:] + lon = nc.variables[nlon][:] + [lon2d,lat2d] = np.meshgrid(lon,lat) + else: + lat = nc.variables[nlat][0,:,:] + lon = nc.variables[nlon][0,:,:] + [lon2d,lat2d] = [lon,lat] + return lon2d,lat2d + +def smooth (field, coeff): + ## actually blur_image could work with different coeff on x and y + if coeff > 1: result = blur_image(field,int(coeff)) + else: result = field + return result + +def gauss_kern(size, sizey=None): + import numpy as np + ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth + # Returns a normalized 2D gauss kernel array for convolutions + size = int(size) + if not sizey: + sizey = size + else: + sizey = int(sizey) + x, y = np.mgrid[-size:size+1, -sizey:sizey+1] + g = np.exp(-(x**2/float(size)+y**2/float(sizey))) + return g / g.sum() + +def blur_image(im, n, ny=None) : + from scipy.signal import convolve + ## FROM COOKBOOK http://www.scipy.org/Cookbook/SignalSmooth + # blurs the image by convolving with a gaussian kernel of typical size n. + # The optional keyword argument ny allows for a different size in the y direction. + g = gauss_kern(n, sizey=ny) + improc = convolve(im, g, mode='same') + return improc + +def getwinddef (nc): + ## getwinds for predefined types + typefile = whatkindfile(nc) + ### + if typefile is 'mesoapi': [uchar,vchar] = ['Um','Vm'] + elif typefile is 'gcm': [uchar,vchar] = ['u','v'] + elif typefile is 'meso': [uchar,vchar] = ['U','V'] + else: [uchar,vchar] = ['not found','not found'] + ### + if typefile in ['meso']: metwind = False ## geometrical (wrt grid) + else: metwind = True ## meteorological (zon/mer) + if metwind is False: print "Not using meteorological winds. You trust numerical grid as being (x,y)" + ### + return uchar,vchar,metwind + +def vectorfield (u, v, x, y, stride=3, scale=15., factor=250., color='black', csmooth=1, key=True): + ## scale regle la reference du vecteur + ## factor regle toutes les longueurs (dont la reference). l'AUGMENTER pour raccourcir les vecteurs. + import matplotlib.pyplot as plt + import numpy as np + posx = np.min(x) - np.std(x) / 10. + posy = np.min(y) - np.std(y) / 10. + u = smooth(u,csmooth) + v = smooth(v,csmooth) + widthvec = 0.003 #0.005 #0.003 + q = plt.quiver( x[::stride,::stride],\ + y[::stride,::stride],\ + u[::stride,::stride],\ + v[::stride,::stride],\ + angles='xy',color=color,pivot='middle',\ + scale=factor,width=widthvec ) + if color in ['white','yellow']: kcolor='black' + else: kcolor=color + if key: p = plt.quiverkey(q,posx,posy,scale,\ + str(int(scale)),coordinates='data',color=kcolor,labelpos='S',labelsep = 0.03) + return + +def display (name): + from os import system + system("display "+name+" > /dev/null 2> /dev/null &") + return name + +def findstep (wlon): + steplon = int((wlon[1]-wlon[0])/4.) #3 + step = 120. + while step > steplon and step > 15. : step = step / 2. + if step <= 15.: + while step > steplon and step > 5. : step = step - 5. + if step <= 5.: + while step > steplon and step > 1. : step = step - 1. + if step <= 1.: + step = 1. + return step + +def define_proj (char,wlon,wlat,back=None): + from mpl_toolkits.basemap import Basemap + import numpy as np + import matplotlib as mpl + from mymath import max + meanlon = 0.5*(wlon[0]+wlon[1]) + meanlat = 0.5*(wlat[0]+wlat[1]) + if wlat[0] >= 80.: blat = 40. + elif wlat[1] <= -80.: blat = -40. + elif wlat[1] >= 0.: blat = wlat[0] + elif wlat[0] <= 0.: blat = wlat[1] + print "blat ", blat + h = 50. ## en km + radius = 3397200. + if char == "cyl": m = Basemap(rsphere=radius,projection='cyl',\ + llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) + elif char == "moll": m = Basemap(rsphere=radius,projection='moll',lon_0=meanlon) + elif char == "ortho": m = Basemap(rsphere=radius,projection='ortho',lon_0=meanlon,lat_0=meanlat) + elif char == "lcc": m = Basemap(rsphere=radius,projection='lcc',lat_1=meanlat,lat_0=meanlat,lon_0=meanlon,\ + llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) + elif char == "npstere": m = Basemap(rsphere=radius,projection='npstere', boundinglat=blat, lon_0=0.) + elif char == "spstere": m = Basemap(rsphere=radius,projection='spstere', boundinglat=blat, lon_0=0.) + elif char == "nplaea": m = Basemap(rsphere=radius,projection='nplaea', boundinglat=wlat[0], lon_0=meanlon) + elif char == "laea": m = Basemap(rsphere=radius,projection='laea',lon_0=meanlon,lat_0=meanlat,lat_ts=meanlat,\ + llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) + elif char == "nsper": m = Basemap(rsphere=radius,projection='nsper',lon_0=meanlon,lat_0=meanlat,satellite_height=h*1000.) + elif char == "merc": m = Basemap(rsphere=radius,projection='merc',lat_ts=0.,\ + llcrnrlat=wlat[0],urcrnrlat=wlat[1],llcrnrlon=wlon[0],urcrnrlon=wlon[1]) + fontsizemer = int(mpl.rcParams['font.size']*3./4.) + if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon) + else: step = 10. + steplon = step*2. + #if back in ["geolocal"]: + # step = np.min([5.,step]) + # steplon = step + print step + m.drawmeridians(np.r_[-180.:180.:steplon], labels=[0,0,0,1], color='grey', fontsize=fontsizemer) + m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer) + if back: m.warpimage(marsmap(back),scale=0.75) + #if not back: + # if not var: back = "mola" ## if no var: draw mola + # elif typefile in ['mesoapi','meso','geo'] \ + # and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw + # else: pass ## else: draw None + return m + +#### test temporaire +def putpoints (map,plot): + #### from http://www.scipy.org/Cookbook/Matplotlib/Maps + # lat/lon coordinates of five cities. + lats = [18.4] + lons = [-134.0] + points=['Olympus Mons'] + # compute the native map projection coordinates for cities. + x,y = map(lons,lats) + # plot filled circles at the locations of the cities. + map.plot(x,y,'bo') + # plot the names of those five cities. + wherept = 0 #1000 #50000 + for name,xpt,ypt in zip(points,x,y): + plot.text(xpt+wherept,ypt+wherept,name) + ## le nom ne s'affiche pas... + return + +def calculate_bounds(field,vmin=None,vmax=None): + import numpy as np + from mymath import max,min,mean + ind = np.where(field < 9e+35) + fieldcalc = field[ ind ] # la syntaxe compacte ne marche si field est un tuple + ### + dev = np.std(fieldcalc)*3.0 + ### + if vmin is None: + zevmin = mean(fieldcalc) - dev + else: zevmin = vmin + ### + if vmax is None: zevmax = mean(fieldcalc) + dev + else: zevmax = vmax + if vmin == vmax: + zevmin = mean(fieldcalc) - dev ### for continuity + zevmax = mean(fieldcalc) + dev ### for continuity + ### + if zevmin < 0. and min(fieldcalc) > 0.: zevmin = 0. + print "field ", min(fieldcalc), max(fieldcalc) + print "bounds ", zevmin, zevmax + return zevmin, zevmax + +def bounds(what_I_plot,zevmin,zevmax): + from mymath import max,min,mean + ### might be convenient to add the missing value in arguments + what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7) + print "new min ", min(what_I_plot) + what_I_plot[ what_I_plot > 9e+35 ] = -9e+35 + what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7) + print "new max ", max(what_I_plot) + return what_I_plot + +def nolow(what_I_plot): + from mymath import max,min + lim = 0.15*0.5*(abs(max(what_I_plot))+abs(min(what_I_plot))) + print "on vire en dessous de ", lim + what_I_plot [ abs(what_I_plot) < lim ] = 1.e40 + return what_I_plot + +def zoomset (wlon,wlat,zoom): + dlon = abs(wlon[1]-wlon[0])/2. + dlat = abs(wlat[1]-wlat[0])/2. + [wlon,wlat] = [ [wlon[0]+zoom*dlon/100.,wlon[1]-zoom*dlon/100.],\ + [wlat[0]+zoom*dlat/100.,wlat[1]-zoom*dlat/100.] ] + print "zoom %",zoom,wlon,wlat + return wlon,wlat + +def fmtvar (whichvar="def"): + fmtvar = { \ + "tk": "%.0f",\ + "tpot": "%.0f",\ + "def": "%.1e",\ + "PTOT": "%.0f",\ + "HGT": "%.1e",\ + "USTM": "%.2f",\ + "HFX": "%.0f",\ + "ICETOT": "%.1e",\ + "TAU_ICE": "%.2f",\ + "VMR_ICE": "%.1e",\ + "MTOT": "%.0f",\ + "anomaly": "%.1f",\ + "W": "%.1f",\ + } + if whichvar not in fmtvar: + whichvar = "def" + return fmtvar[whichvar] + +#################################################################################################################### +### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png +def defcolorb (whichone="def"): + whichcolorb = { \ + "def": "spectral",\ + "HGT": "spectral",\ + "tk": "gist_heat",\ + "QH2O": "PuBu",\ + "USTM": "YlOrRd",\ + "HFX": "RdYlBu",\ + "ICETOT": "YlGnBu",\ + "MTOT": "PuBu",\ + "TAU_ICE": "Blues",\ + "VMR_ICE": "Blues",\ + "W": "jet",\ + "anomaly": "RdBu_r",\ + } +#W --> spectral ou jet +#spectral BrBG RdBu_r + print "predefined colorbars" + if whichone not in whichcolorb: + whichone = "def" + return whichcolorb[whichone] + +def definecolorvec (whichone="def"): + whichcolor = { \ + "def": "black",\ + "vis": "yellow",\ + "vishires": "yellow",\ + "molabw": "yellow",\ + "mola": "black",\ + "gist_heat": "white",\ + "hot": "tk",\ + "gist_rainbow": "black",\ + "spectral": "black",\ + "gray": "red",\ + "PuBu": "black",\ + } + if whichone not in whichcolor: + whichone = "def" + return whichcolor[whichone] + +def marsmap (whichone="vishires"): + from os import uname + mymachine = uname()[1] + ### not sure about speed-up with this method... looks the same + if "lmd.jussieu.fr" in mymachine: domain = "/u/aslmd/WWW/maps/" + else: domain = "http://www.lmd.jussieu.fr/~aslmd/maps/" + whichlink = { \ + #"vis": "http://maps.jpl.nasa.gov/pix/mar0kuu2.jpg",\ + #"vishires": "http://www.lmd.jussieu.fr/~aslmd/maps/MarsMap_2500x1250.jpg",\ + #"geolocal": "http://dl.dropbox.com/u/11078310/geolocal.jpg",\ + #"mola": "http://www.lns.cornell.edu/~seb/celestia/mars-mola-2k.jpg",\ + #"molabw": "http://dl.dropbox.com/u/11078310/MarsElevation_2500x1250.jpg",\ + "vis": domain+"mar0kuu2.jpg",\ + "vishires": domain+"MarsMap_2500x1250.jpg",\ + "geolocal": domain+"geolocal.jpg",\ + "mola": domain+"mars-mola-2k.jpg",\ + "molabw": domain+"MarsElevation_2500x1250.jpg",\ + "clouds": "http://www.johnstonsarchive.net/spaceart/marswcloudmap.jpg",\ + "jupiter": "http://www.mmedia.is/~bjj/data/jupiter_css/jupiter_css.jpg",\ + "jupiter_voy": "http://www.mmedia.is/~bjj/data/jupiter/jupiter_vgr2.jpg",\ + } + ### see http://www.mmedia.is/~bjj/planetary_maps.html + if whichone not in whichlink: + print "marsmap: choice not defined... you'll get the default one... " + whichone = "vishires" + return whichlink[whichone] + +def earthmap (whichone): + if whichone == "contrast": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthMapAtmos_2500x1250.jpg" + elif whichone == "bw": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/EarthElevation_2500x1250.jpg" + elif whichone == "nice": whichlink="http://users.info.unicaen.fr/~karczma/TEACH/InfoGeo/Images/Planets/earthmap1k.jpg" + return whichlink + +def latinterv (area="Whole"): + list = { \ + "Europe": [[ 20., 80.],[- 50., 50.]],\ + "Central_America": [[-10., 40.],[ 230., 300.]],\ + "Africa": [[-20., 50.],[- 50., 50.]],\ + "Whole": [[-90., 90.],[-180.,-180.]],\ + "Southern_Hemisphere": [[-90., 60.],[-180.,-180.]],\ + "Northern_Hemisphere": [[-60., 90.],[-180.,-180.]],\ + "Tharsis": [[-30., 60.],[-170.,- 10.]],\ + "Whole_No_High": [[-60., 60.],[-180., 180.]],\ + "Chryse": [[-60., 60.],[- 60., 60.]],\ + "North_Pole": [[ 50., 90.],[-180., 180.]],\ + "Close_North_Pole": [[ 75., 90.],[-180., 180.]],\ + "Far_South_Pole": [[-90.,-40.],[-180., 180.]],\ + "South_Pole": [[-90.,-50.],[-180., 180.]],\ + "Close_South_Pole": [[-90.,-75.],[-180., 180.]],\ + } + if area not in list: area = "Whole" + [olat,olon] = list[area] + return olon,olat + Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/time.F =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/time.F (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/time.F (revision 261) @@ -0,0 +1,110 @@ +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + real function sol2ls(sol) + +c convert a given martian day number (sol) +c into corresponding solar longitude, Ls (in degr.), +c where sol=0=Ls=0 is the +c northern hemisphere spring equinox. + + implicit none + +c input + real sol +c output + real ls + +c local variables + double precision year_day,peri_day,timeperi,e_elips + double precision pi,radtodeg +c number of martian days (sols) in a martian year + parameter (year_day=668.6d0) +c perihelion date (in sols) + parameter (peri_day=485.35d0) +c orbital eccentricity + parameter (e_elips=0.09340d0) + parameter (pi=3.14159265358979d0) +c radtodeg: 180/pi + parameter (radtodeg=57.2957795130823d0) +c timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 + parameter (timeperi=1.90258341759902d0) + + double precision zanom,xref,zx0,zdx,zteta,zz +c xref: mean anomaly, zx0: eccentric anomaly, zteta: true anomaly + integer iter + + zz=(sol-peri_day)/year_day + zanom=2.*pi*(zz-nint(zz)) + xref=dabs(zanom) + +c The equation zx0 - e * sin (zx0) = xref, solved by Newton + zx0=xref+e_elips*dsin(xref) + do 110 iter=1,10 + zdx=-(zx0-e_elips*dsin(zx0)-xref)/(1.-e_elips*dcos(zx0)) + if(dabs(zdx).le.(1.d-7)) then ! typically, 2 or 3 iterations are enough + goto 120 + endif + zx0=zx0+zdx + 110 continue + 120 continue + zx0=zx0+zdx + if(zanom.lt.0.) zx0=-zx0 + +c compute true anomaly zteta, now that eccentric anomaly zx0 is known + zteta=2.*datan(dsqrt((1.+e_elips)/(1.-e_elips))*dtan(zx0/2.)) + +c compute Ls + ls=real(zteta-timeperi) + if(ls.lt.0.) ls=ls+2.*real(pi) + if(ls.gt.2.*pi) ls=ls-2.*real(pi) +c convert Ls in deg. + ls=real(radtodeg)*ls + + sol2ls = ls + return + end + +cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc + real function ls2sol(ls) + +c Returns solar longitude, Ls (in deg.), from day number (in sol), +c where sol=0=Ls=0 at the northern hemisphere spring equinox + + implicit none + +c Arguments: + real ls + +c Local: + double precision xref,zx0,zteta,zz +c xref: mean anomaly, zteta: true anomaly, zx0: eccentric anomaly + double precision year_day + double precision peri_day,timeperi,e_elips + double precision pi,degrad + parameter (year_day=668.6d0) ! number of sols in a amartian year +c data peri_day /485.0/ + parameter (peri_day=485.35d0) ! date (in sols) of perihelion +c timeperi: 2*pi*( 1 - Ls(perihelion)/ 360 ); Ls(perihelion)=250.99 + parameter (timeperi=1.90258341759902d0) + parameter (e_elips=0.0934d0) ! eccentricity of orbit + parameter (pi=3.14159265358979d0) + parameter (degrad=57.2957795130823d0) + + if (abs(ls).lt.1.0e-5) then + if (ls.ge.0.0) then + ls2sol = 0.0 + else + ls2sol = real(year_day) + end if + return + end if + + zteta = ls/degrad + timeperi + zx0 = 2.0*datan(dtan(0.5*zteta)/dsqrt((1.+e_elips)/(1.-e_elips))) + xref = zx0-e_elips*dsin(zx0) + zz = xref/(2.*pi) + ls2sol = real(zz*year_day + peri_day) + if (ls2sol.lt.0.0) ls2sol = ls2sol + real(year_day) + if (ls2sol.ge.year_day) ls2sol = ls2sol - real(year_day) + + return + end Index: trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/winds.py =================================================================== --- trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/winds.py (revision 261) +++ trunk/MESOSCALE/LMD_MM_MARS/SRC/PYTHON/winds.py (revision 261) @@ -0,0 +1,351 @@ +#!/usr/bin/env python + +### A. Spiga -- LMD -- 30/06/2011 to 10/07/2011 +### Thanks to A. Colaitis for the parser trick + + +#################################### +#################################### +### The main program to plot vectors +def winds (namefile,\ + nvert,\ + proj=None,\ + back=None,\ + target=None, + stride=3,\ + numplot=2,\ + var=None,\ + colorb="def",\ + winds=True,\ + addchar=None,\ + interv=[0,1],\ + vmin=None,\ + vmax=None,\ + tile=False,\ + zoom=None,\ + display=True,\ + itstep=None,\ + hole=False,\ + save="gui",\ + anomaly=False,\ + var2=None,\ + ndiv=10,\ + first=1): + + #################################################################################################################### + ### Colorbars http://www.scipy.org/Cookbook/Matplotlib/Show_colormaps?action=AttachFile&do=get&target=colormaps3.png + + ################################# + ### Load librairies and functions + from netCDF4 import Dataset + from myplot import getcoord2d,define_proj,makeplotres,simplinterv,vectorfield,ptitle,latinterv,getproj,wrfinterv,dumpbdy,\ + fmtvar,definecolorvec,defcolorb,getprefix,putpoints,calculate_bounds,errormess,definesubplot,\ + zoomset,getcoorddef,getwinddef,whatkindfile,reducefield,bounds,getstralt,getfield,smooth,nolow,\ + getname,localtime,polarinterv + from mymath import deg,max,min,mean + from matplotlib.pyplot import contour,contourf, subplot, figure, rcParams, savefig, colorbar, pcolor, show + from matplotlib.cm import get_cmap + import numpy as np + from numpy.core.defchararray import find + + #if save == 'eps': rcParams['backend'] = 'PS' + + ###################### + ### Load NETCDF object + nc = Dataset(namefile) + + ################################## + ### Initial checks and definitions + typefile = whatkindfile(nc) ## TYPEFILE + if var not in nc.variables: var = False ## VAR + if winds: ## WINDS + [uchar,vchar,metwind] = getwinddef(nc) + if uchar == 'not found': winds = False + if not var and not winds: errormess("please set at least winds or var",printvar=nc.variables) + [lon2d,lat2d] = getcoorddef(nc) ## COORDINATES, could be moved below + if proj == None: proj = getproj(nc) ## PROJECTION + + ########################## + ### Define plot boundaries + ### todo: possible areas in latinterv in argument (ex: "Far_South_Pole") + if proj in ["npstere","spstere"]: [wlon,wlat] = polarinterv(lon2d,lat2d) + elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) + else: [wlon,wlat] = simplinterv(lon2d,lat2d) + if zoom: [wlon,wlat] = zoomset(wlon,wlat,zoom) + + ######################################### + ### Name for title and graphics save file + basename = getname(var=var,winds=winds,anomaly=anomaly) + basename = basename + getstralt(nc,nvert) ## can be moved elsewhere for a more generic routine + + ################################## + ### Open a figure and set subplots + fig = figure() + sub = definesubplot( numplot, fig ) + + ################################# + ### Time loop for plotting device + found_lct = False + itime = first + nplot = 1 + error = False + if itstep is None and numplot > 0: itstep = int(24./numplot) + elif numplot <= 0: itstep = 1 + while error is False: + + ### Which local time ? + ltst = localtime ( interv[0]+itime*interv[1], 0.5*(wlon[0]+wlon[1]) ) + + ### General plot settings + #print itime, int(ltst), numplot, nplot + if numplot >= 1: + if nplot > numplot: break + if numplot > 1: + if typefile not in ['geo']: subplot(sub+nplot-1) + found_lct = True + ### If only one local time is requested (numplot < 0) + elif numplot <= 0: + if int(ltst) + numplot != 0: + itime += 1 + if found_lct is True: break ## because it means LT was found at previous iteration + else: continue ## continue to iterate to find the correct LT + else: + found_lct = True + + ### Map projection + m = define_proj(proj,wlon,wlat,back=back) + x, y = m(lon2d, lat2d) + + #### Contour plot + if var2: + what_I_contour, error = reducefield( getfield(nc,var2), d4=itime, d3=nvert ) + if not error: + if typefile in ['mesoapi','meso']: what_I_contour = dumpbdy(what_I_contour,6) + zevmin, zevmax = calculate_bounds(what_I_contour) + zelevels = np.linspace(zevmin,zevmax,num=20) + if var2 == 'HGT': zelevels = np.arange(-10000.,30000.,2000.) + contour( x, y, what_I_contour, zelevels, colors='k', linewidths = 0.33 ) #colors='w' )# , alpha=0.5) + + #### Shaded plot + if var: + what_I_plot, error = reducefield( getfield(nc,var), d4=itime, d3=nvert ) + if not error: + fvar = var + ### + if anomaly: + what_I_plot = 100. * ((what_I_plot / smooth(what_I_plot,12)) - 1.) + fvar = 'anomaly' + ### + if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(what_I_plot,6) + zevmin, zevmax = calculate_bounds(what_I_plot,vmin=vmin,vmax=vmax) + if colorb in ["def","nobar"]: palette = get_cmap(name=defcolorb(fvar)) + else: palette = get_cmap(name=colorb) + if not tile: + if not hole: what_I_plot = bounds(what_I_plot,zevmin,zevmax) + zelevels = np.linspace(zevmin*(1. + 1.e-7),zevmax*(1. - 1.e-7)) #,num=20) + contourf( x, y, what_I_plot, zelevels, cmap = palette ) + else: + if hole: what_I_plot = nolow(what_I_plot) + pcolor( x, y, what_I_plot, cmap = palette, \ + vmin=zevmin, vmax=zevmax ) + if colorb != 'nobar' and var != 'HGT': + colorbar(fraction=0.05,pad=0.03,format=fmtvar(fvar),\ + ticks=np.linspace(zevmin,zevmax,ndiv+1),\ + extend='neither',spacing='proportional') + # both min max neither + + ### Vector plot + if winds: + vecx, error = reducefield( getfield(nc,uchar), d4=itime, d3=nvert ) + vecy, error = reducefield( getfield(nc,vchar), d4=itime, d3=nvert ) + if not error: + if typefile in ['mesoapi','meso']: + [vecx,vecy] = [dumpbdy(vecx,6,stag=uchar), dumpbdy(vecy,6,stag=vchar)] + key = True + elif typefile in ['gcm']: + key = False + if metwind: [vecx,vecy] = m.rotate_vector(vecx, vecy, lon2d, lat2d) + if var == False: colorvec = definecolorvec(back) + else: colorvec = definecolorvec(colorb) + vectorfield(vecx, vecy,\ + x, y, stride=stride, csmooth=2,\ + #scale=15., factor=300., color=colorvec, key=key) + scale=20., factor=250., color=colorvec, key=key) + #200. ## or csmooth=stride + + ### Next subplot + plottitle = basename + if typefile in ['mesoapi','meso']: + if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) + else: plottitle = plottitle + "_LT"+str(ltst) + ptitle( plottitle ) + itime += itstep + nplot += 1 + + ########################################################################## + ### Save the figure in a file in the data folder or an user-defined folder + if typefile in ['meso','mesoapi']: prefix = getprefix(nc) + elif typefile in ['gcm','gcmex']: prefix = 'LMD_GCM_' + else: prefix = '' + ### + zeplot = prefix + basename + if addchar: zeplot = zeplot + addchar + if numplot <= 0: zeplot = zeplot + "_LT"+str(abs(numplot)) + ### + if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot + else: zeplot = target + "/" + zeplot + ### + if found_lct: + pad_inches_value = 0.35 + if save == 'png': + makeplotres(zeplot,res=100.,pad_inches_value=pad_inches_value) #,erase=True) ## a miniature + makeplotres(zeplot,res=200.,pad_inches_value=pad_inches_value,disp=False) + elif save in ['eps','svg','pdf']: + makeplotres(zeplot, pad_inches_value=pad_inches_value,disp=False,ext=save) + elif save == 'gui': + show() + else: + print "save mode not supported. using gui instead." + show() + else: print "Local time not found" + + ############### + ### Now the end + return zeplot + +############################## +### A specific stuff for below +def adjust_length (tab, zelen): + from numpy import ones + if tab is None: + outtab = ones(zelen) * -999999 + else: + if zelen != len(tab): + print "not enough or too much values... setting same values all variables" + outtab = ones(zelen) * tab[0] + else: + outtab = tab + return outtab + +########################################################################################### +########################################################################################### +### What is below relate to running the file as a command line executable (very convenient) +if __name__ == "__main__": + import sys + from optparse import OptionParser ### to be replaced by argparse + from api_wrapper import api_onelevel + from netCDF4 import Dataset + from myplot import getlschar + from os import system + + ############################# + ### Get options and variables + parser = OptionParser() + parser.add_option('-f', '--file', action='append',dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)') + parser.add_option('-l', '--level', action='store',dest='nvert', type="float", default=0, help='level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)') + parser.add_option('-p', '--proj', action='store',dest='proj', type="string", default=None, help='projection') + parser.add_option('-b', '--back', action='store',dest='back', type="string", default=None, help='background image (def: None)') + parser.add_option('-t', '--target', action='store',dest='target', type="string", default=None, help='destination folder') + parser.add_option('-s', '--stride', action='store',dest='stride', type="int", default=3, help='stride vectors (def=3)') + parser.add_option('-v', '--var', action='append',dest='var', type="string", default=None, help='variable color-shaded (append)') + parser.add_option('-n', '--num', action='store',dest='numplot', type="int", default=2, help='plot number (def=2)(<0: plot LT -*numplot*)') + parser.add_option('-i', '--interp', action='store',dest='interp', type="int", default=None, help='interpolation (2: p, 3: z-amr, 4:z-als)') + parser.add_option('-c', '--color', action='store',dest='colorb', type="string", default="def", help='change colormap (nobar: no colorbar)') + parser.add_option('-x', '--no-vect',action='store_false',dest='winds', default=True, help='no wind vectors') + parser.add_option('-m', '--min', action='append',dest='vmin', type="float", default=None, help='bounding minimum value (append)') + parser.add_option('-M', '--max', action='append',dest='vmax', type="float", default=None, help='bounding maximum value (append)') + parser.add_option('-T', '--tiled', action='store_true',dest='tile', default=False, help='draw a tiled plot (no blank zone)') + parser.add_option('-z', '--zoom', action='store',dest='zoom', type="float", default=None, help='zoom factor in %') + parser.add_option('-N', '--no-api', action='store_true',dest='nocall', default=False, help='do not recreate api file') + parser.add_option('-d', '--display',action='store_false',dest='display', default=True, help='do not pop up created images') + parser.add_option('-e', '--itstep', action='store',dest='itstep', type="int", default=None, help='stride time (def=4)') + parser.add_option('-H', '--hole', action='store_true',dest='hole', default=False, help='holes above max and below min') + parser.add_option('-S', '--save', action='store',dest='save', type="string", default="gui", help='save mode (png,eps,svg,pdf or gui)(def=gui)') + parser.add_option('-a', '--anomaly',action='store_true',dest='anomaly', default=False, help='compute and plot relative anomaly in %') + parser.add_option('-w', '--with', action='store',dest='var2', type="string", default=None, help='variable contoured') + parser.add_option('--div', action='store',dest='ndiv', type="int", default=10, help='number of divisions in colorbar (def: 10)') + parser.add_option('-F', '--first', action='store',dest='first', type="int", default=1, help='first subscript to plot (def: 1)') + #parser.add_option('-V', action='store', dest='comb', type="float", default=None, help='a defined combination of variables') + (opt,args) = parser.parse_args() + if opt.namefile is None: + print "I want to eat one file at least ! Use winds.py -f name_of_my_file. Or type winds.py -h" + exit() + if opt.var is None and opt.anomaly is True: + print "Cannot ask to compute anomaly if no variable is set" + exit() + print "Options:", opt + + listvar = '' + if opt.var is None: + zerange = [-999999] + else: + zelen = len(opt.var) + zerange = range(zelen) + #if zelen == 1: listvar = opt.var[0] + ',' + #else : + for jjj in zerange: listvar += opt.var[jjj] + ',' + listvar = listvar[0:len(listvar)-1] + vmintab = adjust_length (opt.vmin, zelen) + vmaxtab = adjust_length (opt.vmax, zelen) + + for i in range(len(opt.namefile)): + + zefile = opt.namefile[i] + print zefile + zelevel = opt.nvert + stralt = None + [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) + + ##################################################### + ### Call Fortran routines for vertical interpolations + if opt.interp is not None: + if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 + ### winds or no winds + if opt.winds : zefields = 'uvmet' + else : zefields = '' + ### var or no var + #if opt.var is None : pass + if zefields == '' : zefields = listvar + else : zefields = zefields + "," + listvar + if opt.var2 is not None : zefields = zefields + "," + opt.var2 + print zefields + zefile = api_onelevel ( path_to_input = '', \ + input_name = zefile, \ + fields = zefields, \ + interp_method = opt.interp, \ + onelevel = zelevel, \ + nocall = opt.nocall ) + print zefile + zelevel = 0 ## so that zelevel could play again the role of nvert + + if opt.var is None: zerange = [-999999] + else: zerange = range(zelen) + for jjj in zerange: + if jjj == -999999: + argvar = None + argvmin = None + argvmax = None + else: + argvar = opt.var[jjj] + if vmintab[jjj] != -999999: argvmin = vmintab[jjj] + else: argvmin = None + if vmaxtab[jjj] != -999999: argvmax = vmaxtab[jjj] + else: argvmax = None + ############# + ### Main call + name = winds (zefile,int(zelevel),\ + proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=argvar,\ + numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ + addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,\ + tile=opt.tile,zoom=opt.zoom,display=opt.display,\ + itstep=opt.itstep,hole=opt.hole,save=opt.save,\ + anomaly=opt.anomaly,var2=opt.var2,ndiv=opt.ndiv,first=opt.first) + print 'Done: '+name + system("rm -f to_be_erased") + + ######################################################### + ### Generate a .sh file with the used command saved in it + command = "" + for arg in sys.argv: command = command + arg + ' ' + f = open(name+'.sh', 'w') + f.write(command)