Index: /trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90 =================================================================== --- /trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90 (revision 224) +++ /trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90 (revision 225) @@ -274,5 +274,5 @@ ENDIF END DO - CLOSE ( 112 ) + CLOSE ( 111 ) !!AS: don't know why but it was written 112... fixed the bug. print*, " " ! We clean up our own messes. @@ -1166,4 +1166,8 @@ ENDDO + !!!! forgotten in the initial program + DEALLOCATE ( input_file_names ) + DEALLOCATE ( output_file_names ) + ! ! WELL... Index: /trunk/MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic =================================================================== --- /trunk/MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic (revision 224) +++ /trunk/MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic (revision 225) @@ -5,5 +5,6 @@ # ehouarn v2.0 -tar xvf netcdf-3.6.1.tar +tar xzvf netcdf-3.6.1.tar.gz +tar xvf netcdf-3.6.1.tar cd netcdf-3.6.1 Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_g95.sh =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_g95.sh (revision 224) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_g95.sh (revision 225) @@ -8,4 +8,4 @@ \rm logc \rm loge -f2py -c -m api ../../../LMD_MM_MARS/SRC/POSTPROC/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 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_DEV/PLOT/PYTHON/mylib/api_pgf90.sh =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_pgf90.sh (revision 224) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_pgf90.sh (revision 225) @@ -5,5 +5,6 @@ #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/NETCDF/pgf90_64_fpic/netcdf-3.6.1/ +NETCDF=/donnees/aslmd/MODELES/MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic-3.6.1/ + echo $NETCDF @@ -12,4 +13,4 @@ \rm logc \rm loge -f2py -c -m api ../../../LMD_MM_MARS/SRC/POSTPROC/api.F90 --fcompiler=pg -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include --f90flags="-Mfree" > logc 2> 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_DEV/PLOT/PYTHON/mylib/api_wrapper.py =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_wrapper.py (revision 224) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_wrapper.py (revision 225) @@ -26,5 +26,5 @@ if interp_method == 4: output_name = input_name+'_zabg' - print input_name, output_name + #print input_name, output_name if nocall: pass Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py (revision 224) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py (revision 225) @@ -44,5 +44,7 @@ from timestuff import sol2ls nc = Dataset(namefile) - if 'Times' in nc.variables and 'vert' not in nc.variables: + if namefile[0]+namefile[1]+namefile[2] != "geo" \ + and 'Times' in nc.variables \ + and 'vert' not in nc.variables: zetime = nc.variables['Times'][0] zetimestart = getattr(nc, 'START_DATE') @@ -100,5 +102,14 @@ nx = len(lon2d[0,:])-1 ny = len(lon2d[:,0])-1 - return [[lon2d[0,0],lon2d[nx,ny]],[lat2d[0,0],lat2d[nx,ny]]] + lon1 = lon2d[0,0] + lon2 = lon2d[nx,ny] + lat1 = lat2d[0,0] + lat2 = lat2d[nx,ny] + wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 + if lon1 < lon2: wlon = [lon1, lon2 + wider] ## a tester en normal + else: wlon = [lon2, lon1 + wider] + if lat1 < lat2: wlat = [lat1, lat2] + else: wlat = [lat2, lat1] + return [wlon,wlat] def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True): @@ -237,4 +248,5 @@ if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon) else: step = 10. + print step m.drawmeridians(np.r_[-180.:180.:step*2.], 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) @@ -252,4 +264,5 @@ "HGT": "%.1e",\ "USTM": "%.2f",\ + "HFX": "%.0f",\ } if whichvar not in fmtvar: @@ -264,5 +277,5 @@ "QH2O": "PuBu",\ "USTM": "YlOrRd",\ -#"RdPu",\ + "HFX": "RdYlBu",\ } if whichone not in whichcolorb: Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py (revision 224) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py (revision 225) @@ -3,7 +3,7 @@ ### A. Spiga -- LMD -- 30/06/2011 -- slight modif early 07/2011 -def domain (namefile,proj="ortho",back="vishires",target=None,var='HGT'): +def domain (namefile,proj=None,back="vishires",target=None): from netCDF4 import Dataset - from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,getprefix,dumpbdy + from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,getprefix,dumpbdy,getproj,latinterv,wrfinterv,simplinterv from matplotlib.pyplot import contourf,rcParams from numpy.core.defchararray import find @@ -11,8 +11,21 @@ nc = Dataset(namefile) ### - [lon2d,lat2d] = getcoord2d(nc) + 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) lat2d = dumpbdy(lat2d) - [wlon,wlat] = simplinterv(lon2d,lat2d) + 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) @@ -22,5 +35,4 @@ contourf(x, y, what_I_plot, 40) ### - zeplot = getprefix(nc) + "domain" if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot else: zeplot = target + "/" + zeplot @@ -36,8 +48,7 @@ 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="ortho", help='projection') + 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') - parser.add_option('-v', action='store', dest='var', type="string", default='HGT', help='variable contoured') (opt,args) = parser.parse_args() if opt.namefile is None: @@ -45,3 +56,3 @@ exit() print "Options:", opt - domain (opt.namefile,proj=opt.proj,back=opt.back,target=opt.target,var=opt.var) + domain (opt.namefile,proj=opt.proj,back=opt.back,target=opt.target) Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/makefig.sh =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/makefig.sh (revision 225) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/makefig.sh (revision 225) @@ -0,0 +1,14 @@ +#! /bin/bash + + fold="/home/aslmd/POLAR/" + file1="$fold/POLAR_THOMAS_ls50/wrfout_d03_2024-02-45_06:00:00" + file2="$fold/POLAR_THOMAS_ls20/wrfout_d03_2024-01-42_06:00:00" + file3="$fold/POLAR_THOMAS_ls0/wrfout_d03_2024-01-06_06:00:00" + dest="/u/aslmd/WWW/antichambre/thomas/" + + domain.py -t $dest -f $file1 -p npstere + winds.py -t $dest -p lcc -n 2 -z 25 -i 4 -s 2 \ + -v HFX -m -15. -M 0.0 \ + -v USTM -m 0. -M 0.8 \ + -v HGT -m 0 -M 0 \ + -f $file1 -f $file2 -f $file3 Index: /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py =================================================================== --- /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py (revision 224) +++ /trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py (revision 225) @@ -14,5 +14,5 @@ target=None, stride=3,\ - numplot=4,\ + numplot=2,\ var=None,\ colorb=True,\ @@ -52,4 +52,5 @@ elif 'vert' in nc.variables: typefile = 'mesoapi' elif 'U' in nc.variables: typefile = 'meso' + elif 'HGT_M' in nc.variables: typefile = 'geo' else: print "typefile not supported." @@ -59,5 +60,5 @@ ############################################################## ### Try to guess the projection from wrfout if not set by user - if typefile in ['mesoapi','meso']: + if typefile in ['mesoapi','meso','geo']: if proj == None: proj = getproj(nc) ### (il faudrait passer CEN_LON dans la projection ?) @@ -70,6 +71,6 @@ if not back: if not var: back = "mola" ## if no var: draw mola - elif typefile in ['mesoapi','meso'] \ - and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw + 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 @@ -82,7 +83,11 @@ elif typefile in ['gcm']: [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) + elif typefile in ['geo']: + [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) else: [wlon,wlat] = simplinterv(lon2d,lat2d) + print wlon + print wlat if zoom: dlon = abs(wlon[1]-wlon[0])/2. @@ -105,4 +110,6 @@ metwind = False ## geometrical (wrt grid) print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)." + elif typefile is 'geo': + winds = None ##################################################### @@ -122,4 +129,8 @@ if vmax is None: zevmax = mean(field) + dev else: zevmax = vmax + if vmin == vmax: + zevmin = mean(field) - dev ### for continuity + zevmax = mean(field) + dev ### for continuity + print "field ", min(field), max(field) print "bounds ", zevmin, zevmax ### some already defined colormaps @@ -137,9 +148,10 @@ ######################################### ### Name for title and graphics save file - if winds: basename = "UV_" - else: basename = "" - if var: basename = basename + var - ### - if dimension == 4: + if var and winds: basename = var + '_UV' + elif var: basename = var + elif winds: basename = 'UV' + else: exit() + ### + if dimension == 4 or winds: if typefile is 'meso': stralt = "_lvl" + str(nvert) elif typefile is 'mesoapi': @@ -158,4 +170,6 @@ ### Open a figure and set subplots fig = figure() + rcParams['font.size'] = 12. + if typefile in ['geo']: numplot = 1 if numplot > 0: if numplot == 4: @@ -228,4 +242,5 @@ if var: if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(field[i,:,:]) + elif typefile in ['geo']: what_I_plot = field[0,:,:] elif typefile in ['gcm']: if dimension == 2: what_I_plot = field[:,:] @@ -233,16 +248,21 @@ palette = get_cmap(name=colorb) #palette.set_over('b', 1.0) + #print np.array(x).shape + #print np.array(y).shape + #print np.array(what_I_plot).shape if not tile: zelevels = np.linspace(zevmin,zevmax) + what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7) + what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7) contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels ) else: pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax ) - if var in ['HGT']: pass + if var in ['HGT']: pass elif colorb: - ndiv = 10 - colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\ - ticks=np.linspace(zevmin,zevmax,ndiv+1),\ - extend='max',spacing='proportional') - # both min max neither + ndiv = 10 + colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\ + ticks=np.linspace(zevmin,zevmax,ndiv+1),\ + extend='max',spacing='proportional') + # both min max neither ### Vector plot if winds: @@ -257,12 +277,13 @@ else: colorvec = definecolorvec(colorb) vectorfield(vecx, vecy,\ - x, y, stride=stride, csmooth=stride,\ + x, y, stride=stride, csmooth=2,\ scale=15., factor=300., color=colorvec, key=key) - #200. ## or csmooth=2 + #200. ## or csmooth=stride ### Next subplot plottitle = basename - if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) - else: plottitle = plottitle + "_LT"+str(ltst) + if typefile in ['mesoapi','meso']: + if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) + else: plottitle = plottitle + "_LT"+str(ltst) ptitle( plottitle ) sub += 1 @@ -297,9 +318,10 @@ from netCDF4 import Dataset from myplot import getlschar + from numpy import ones ############################# ### Get options and variables parser = OptionParser() - parser.add_option('-f', action='store', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file') + parser.add_option('-f', action='append', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)') parser.add_option('-l', action='store', dest='nvert', type="float", default=0, help='vertical level (def=0)(-i 2: p,mbar)(-i 3,4: z,km)') parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection') @@ -307,11 +329,11 @@ parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder') parser.add_option('-s', action='store', dest='stride', type="int", default=3, help='stride vectors (def=3)') - parser.add_option('-v', action='store', dest='var', type="string", default=None, help='variable contoured') - parser.add_option('-n', action='store', dest='numplot', type="int", default=4, help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)') + parser.add_option('-v', action='append', dest='var', type="string", default=None, help='variable contoured (append)') + parser.add_option('-n', action='store', dest='numplot', type="int", default=2, help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)') parser.add_option('-i', action='store', dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)') parser.add_option('-c', action='store', dest='colorb', type="string", default=True, help='change colormap') parser.add_option('-x', action='store_false', dest='winds', default=True, help='no wind vectors') - parser.add_option('-m', action='store', dest='vmin', type="float", default=None, help='bounding minimum value for color plot') - parser.add_option('-M', action='store', dest='vmax', type="float", default=None, help='bounding maximum value for color plot') + parser.add_option('-m', action='append', dest='vmin', type="float", default=None, help='bounding minimum value (append)') + parser.add_option('-M', action='append', dest='vmax', type="float", default=None, help='bounding maximum value (append)') parser.add_option('-T', action='store_true', dest='tile', default=False, help='draw a tiled plot (no blank zone)') parser.add_option('-z', action='store', dest='zoom', type="float", default=None, help='zoom factor in %') @@ -323,42 +345,85 @@ exit() print "Options:", opt - zefile = opt.namefile - 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 - elif zefields == '' : zefields = opt.var - else : zefields = zefields + "," + opt.var - print zefields - zefile = api_onelevel ( path_to_input = '', \ - input_name = zefile, \ - path_to_output = opt.target, \ - 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 - - ############# - ### Main call - name = winds (zefile,int(zelevel),\ - proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ - addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom) - print 'Done: '+name - - ######################################################### - ### 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) + + if opt.var is None: + pass + else: + listvar = '' + zelen = len(opt.var) + if zelen == 1: listvar = opt.var[0] + ',' + else : + for jjj in range(zelen): listvar += opt.var[jjj] + ',' + listvar = listvar[0:len(listvar)-1] + if opt.vmin: + if zelen != len(opt.vmin): + print "not enough or too much vmin values... setting same values all variables" + vmintab = ones(zelen) * opt.vmin[0] + else: + vmintab = opt.vmin + else: + vmintab = None + if opt.vmax: + if zelen != len(opt.vmax): + print "not enough or too much vmax values... setting same values all variables" + vmaxtab = ones(zelen) * opt.vmax[0] + else: + vmaxtab = opt.vmax + else: + vmaxtab = None + + 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 + elif zefields == '' : zefields = listvar + else : zefields = zefields + "," + listvar + print zefields + zefile = api_onelevel ( path_to_input = '', \ + input_name = zefile, \ + #path_to_output = opt.target, \ + 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: + ############# + ### Main call + name = winds (zefile,int(zelevel),\ + proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ + addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom) + print 'Done: '+name + else: + for jjj in range(len(opt.var)): + if vmintab: argvmin = vmintab[jjj] + else: argvmin = None + if vmaxtab: 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=opt.var[jjj],numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ + addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,tile=opt.tile,zoom=opt.zoom) + print 'Done: '+name + + ######################################################### + ### 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)