Changeset 225
- Timestamp:
- Jul 15, 2011, 2:47:21 PM (13 years ago)
- Location:
- trunk
- Files:
-
- 1 added
- 8 edited
-
MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90 (modified) (2 diffs)
-
MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic (modified) (1 diff)
-
MESOSCALE_DEV/PLOT/PYTHON/mylib/api_g95.sh (modified) (1 diff)
-
MESOSCALE_DEV/PLOT/PYTHON/mylib/api_pgf90.sh (modified) (2 diffs)
-
MESOSCALE_DEV/PLOT/PYTHON/mylib/api_wrapper.py (modified) (1 diff)
-
MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py (modified) (5 diffs)
-
MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py (modified) (5 diffs)
-
MESOSCALE_DEV/PLOT/PYTHON/scripts/makefig.sh (added)
-
MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/MESOSCALE/LMD_MM_MARS/SRC/POSTPROC/api.F90
r202 r225 274 274 ENDIF 275 275 END DO 276 CLOSE ( 11 2 )276 CLOSE ( 111 ) !!AS: don't know why but it was written 112... fixed the bug. 277 277 print*, " " 278 278 ! We clean up our own messes. … … 1166 1166 ENDDO 1167 1167 1168 !!!! forgotten in the initial program 1169 DEALLOCATE ( input_file_names ) 1170 DEALLOCATE ( output_file_names ) 1171 1168 1172 ! 1169 1173 ! WELL... -
trunk/MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic
r191 r225 5 5 # ehouarn v2.0 6 6 7 tar xvf netcdf-3.6.1.tar 7 tar xzvf netcdf-3.6.1.tar.gz 8 tar xvf netcdf-3.6.1.tar 8 9 9 10 cd netcdf-3.6.1 -
trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_g95.sh
r195 r225 8 8 \rm logc 9 9 \rm loge 10 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> loge10 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 11 11 f2py -c -m timestuff time.F --fcompiler=g95 >> logc 2>> loge -
trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_pgf90.sh
r195 r225 5 5 #NETCDF=/distrib/local/netcdf-4.0.1/pgi_7.1-6_64/ 6 6 #NETCDF=/distrib/local/netcdf-3.6.0-p1/pgi_64bits/ 7 NETCDF=/donnees/aslmd/MODELES/MESOSCALE/NETCDF/pgf90_64_fpic/netcdf-3.6.1/ 7 NETCDF=/donnees/aslmd/MODELES/MESOSCALE_DEV/NETCDF/pgf90_64_netcdf_fpic-3.6.1/ 8 8 9 echo $NETCDF 9 10 … … 12 13 \rm logc 13 14 \rm loge 14 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> loge15 f2py -c -m api api.F90 --fcompiler=pg -L$NETCDF/lib -lnetcdf -lm -I$NETCDF/include --f90flags="-Mfree" > logc 2> loge 15 16 f2py -c -m timestuff time.F --fcompiler=pg >> logc 2>> loge -
trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/api_wrapper.py
r207 r225 26 26 if interp_method == 4: output_name = input_name+'_zabg' 27 27 28 print input_name, output_name28 #print input_name, output_name 29 29 30 30 if nocall: pass -
trunk/MESOSCALE_DEV/PLOT/PYTHON/mylib/myplot.py
r207 r225 44 44 from timestuff import sol2ls 45 45 nc = Dataset(namefile) 46 if 'Times' in nc.variables and 'vert' not in nc.variables: 46 if namefile[0]+namefile[1]+namefile[2] != "geo" \ 47 and 'Times' in nc.variables \ 48 and 'vert' not in nc.variables: 47 49 zetime = nc.variables['Times'][0] 48 50 zetimestart = getattr(nc, 'START_DATE') … … 100 102 nx = len(lon2d[0,:])-1 101 103 ny = len(lon2d[:,0])-1 102 return [[lon2d[0,0],lon2d[nx,ny]],[lat2d[0,0],lat2d[nx,ny]]] 104 lon1 = lon2d[0,0] 105 lon2 = lon2d[nx,ny] 106 lat1 = lat2d[0,0] 107 lat2 = lat2d[nx,ny] 108 wider = 0.5 * (abs(lon1)+abs(lon2)) * 0.1 109 if lon1 < lon2: wlon = [lon1, lon2 + wider] ## a tester en normal 110 else: wlon = [lon2, lon1 + wider] 111 if lat1 < lat2: wlat = [lat1, lat2] 112 else: wlat = [lat2, lat1] 113 return [wlon,wlat] 103 114 104 115 def makeplotpngres (filename,res,pad_inches_value=0.25,folder='',disp=True): … … 237 248 if char in ["cyl","lcc","merc","nsper","laea"]: step = findstep(wlon) 238 249 else: step = 10. 250 print step 239 251 m.drawmeridians(np.r_[-180.:180.:step*2.], labels=[0,0,0,1], color='grey', fontsize=fontsizemer) 240 252 m.drawparallels(np.r_[-90.:90.:step], labels=[1,0,0,0], color='grey', fontsize=fontsizemer) … … 252 264 "HGT": "%.1e",\ 253 265 "USTM": "%.2f",\ 266 "HFX": "%.0f",\ 254 267 } 255 268 if whichvar not in fmtvar: … … 264 277 "QH2O": "PuBu",\ 265 278 "USTM": "YlOrRd",\ 266 #"RdPu",\279 "HFX": "RdYlBu",\ 267 280 } 268 281 if whichone not in whichcolorb: -
trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/domain.py
r207 r225 3 3 ### A. Spiga -- LMD -- 30/06/2011 -- slight modif early 07/2011 4 4 5 def domain (namefile,proj= "ortho",back="vishires",target=None,var='HGT'):5 def domain (namefile,proj=None,back="vishires",target=None): 6 6 from netCDF4 import Dataset 7 from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,getprefix,dumpbdy 7 from myplot import getcoord2d,define_proj,makeplotpng,simplinterv,getprefix,dumpbdy,getproj,latinterv,wrfinterv,simplinterv 8 8 from matplotlib.pyplot import contourf,rcParams 9 9 from numpy.core.defchararray import find … … 11 11 nc = Dataset(namefile) 12 12 ### 13 [lon2d,lat2d] = getcoord2d(nc) 13 if proj == None: proj = "ortho" #proj = getproj(nc) 14 ### 15 prefix = namefile[0] + namefile[1] + namefile[2] 16 if prefix == "geo": 17 [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') 18 var = 'HGT_M' 19 zeplot = "domain" 20 else: 21 [lon2d,lat2d] = getcoord2d(nc) 22 var = "HGT" 23 zeplot = getprefix(nc) + "domain" 24 ### 14 25 lon2d = dumpbdy(lon2d) 15 26 lat2d = dumpbdy(lat2d) 16 [wlon,wlat] = simplinterv(lon2d,lat2d) 27 if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") 28 elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) 29 else: [wlon,wlat] = simplinterv(lon2d,lat2d) 17 30 ### 18 31 m = define_proj(proj,wlon,wlat,back=back) … … 22 35 contourf(x, y, what_I_plot, 40) 23 36 ### 24 zeplot = getprefix(nc) + "domain"25 37 if not target: zeplot = namefile[0:find(namefile,'wrfout')] + zeplot 26 38 else: zeplot = target + "/" + zeplot … … 36 48 parser = OptionParser() 37 49 parser.add_option('-f', action='store', dest='namefile', type="string", default=None, help='name of WRF file [NEEDED]') 38 parser.add_option('-p', action='store', dest='proj', type="string", default= "ortho",help='projection')50 parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection') 39 51 parser.add_option('-b', action='store', dest='back', type="string", default="vishires", help='background') 40 52 parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder') 41 parser.add_option('-v', action='store', dest='var', type="string", default='HGT', help='variable contoured')42 53 (opt,args) = parser.parse_args() 43 54 if opt.namefile is None: … … 45 56 exit() 46 57 print "Options:", opt 47 domain (opt.namefile,proj=opt.proj,back=opt.back,target=opt.target ,var=opt.var)58 domain (opt.namefile,proj=opt.proj,back=opt.back,target=opt.target) -
trunk/MESOSCALE_DEV/PLOT/PYTHON/scripts/winds.py
r207 r225 14 14 target=None, 15 15 stride=3,\ 16 numplot= 4,\16 numplot=2,\ 17 17 var=None,\ 18 18 colorb=True,\ … … 52 52 elif 'vert' in nc.variables: typefile = 'mesoapi' 53 53 elif 'U' in nc.variables: typefile = 'meso' 54 elif 'HGT_M' in nc.variables: typefile = 'geo' 54 55 else: 55 56 print "typefile not supported." … … 59 60 ############################################################## 60 61 ### Try to guess the projection from wrfout if not set by user 61 if typefile in ['mesoapi','meso' ]:62 if typefile in ['mesoapi','meso','geo']: 62 63 if proj == None: proj = getproj(nc) 63 64 ### (il faudrait passer CEN_LON dans la projection ?) … … 70 71 if not back: 71 72 if not var: back = "mola" ## if no var: draw mola 72 elif typefile in ['mesoapi','meso' ] \73 and proj not in ['merc','lcc','nsper','laea']:back = "molabw" ## if var but meso: draw molabw73 elif typefile in ['mesoapi','meso','geo'] \ 74 and proj not in ['merc','lcc','nsper','laea']: back = "molabw" ## if var but meso: draw molabw 74 75 else: pass ## else: draw None 75 76 … … 82 83 elif typefile in ['gcm']: 83 84 [lon2d,lat2d] = getcoord2d(nc,nlat="latitude",nlon="longitude",is1d=True) 85 elif typefile in ['geo']: 86 [lon2d,lat2d] = getcoord2d(nc,nlat='XLAT_M',nlon='XLONG_M') 84 87 if proj == "npstere": [wlon,wlat] = latinterv("North_Pole") 85 88 elif proj in ["lcc","laea"]: [wlon,wlat] = wrfinterv(lon2d,lat2d) 86 89 else: [wlon,wlat] = simplinterv(lon2d,lat2d) 90 print wlon 91 print wlat 87 92 if zoom: 88 93 dlon = abs(wlon[1]-wlon[0])/2. … … 105 110 metwind = False ## geometrical (wrt grid) 106 111 print "Beware ! Not using meteorological winds. You trust numerical grid as being (x,y)." 112 elif typefile is 'geo': 113 winds = None 107 114 108 115 ##################################################### … … 122 129 if vmax is None: zevmax = mean(field) + dev 123 130 else: zevmax = vmax 131 if vmin == vmax: 132 zevmin = mean(field) - dev ### for continuity 133 zevmax = mean(field) + dev ### for continuity 134 print "field ", min(field), max(field) 124 135 print "bounds ", zevmin, zevmax 125 136 ### some already defined colormaps … … 137 148 ######################################### 138 149 ### Name for title and graphics save file 139 if winds: basename = "UV_" 140 else: basename = "" 141 if var: basename = basename + var 142 ### 143 if dimension == 4: 150 if var and winds: basename = var + '_UV' 151 elif var: basename = var 152 elif winds: basename = 'UV' 153 else: exit() 154 ### 155 if dimension == 4 or winds: 144 156 if typefile is 'meso': stralt = "_lvl" + str(nvert) 145 157 elif typefile is 'mesoapi': … … 158 170 ### Open a figure and set subplots 159 171 fig = figure() 172 rcParams['font.size'] = 12. 173 if typefile in ['geo']: numplot = 1 160 174 if numplot > 0: 161 175 if numplot == 4: … … 228 242 if var: 229 243 if typefile in ['mesoapi','meso']: what_I_plot = dumpbdy(field[i,:,:]) 244 elif typefile in ['geo']: what_I_plot = field[0,:,:] 230 245 elif typefile in ['gcm']: 231 246 if dimension == 2: what_I_plot = field[:,:] … … 233 248 palette = get_cmap(name=colorb) 234 249 #palette.set_over('b', 1.0) 250 #print np.array(x).shape 251 #print np.array(y).shape 252 #print np.array(what_I_plot).shape 235 253 if not tile: 236 254 zelevels = np.linspace(zevmin,zevmax) 255 what_I_plot[ what_I_plot < zevmin ] = zevmin*(1. + 1.e-7) 256 what_I_plot[ what_I_plot > zevmax ] = zevmax*(1. - 1.e-7) 237 257 contourf( x, y, what_I_plot, 10, cmap = palette, levels = zelevels ) 238 258 else: 239 259 pcolor( x, y, what_I_plot, cmap = palette, vmin=zevmin, vmax=zevmax ) 240 if var in ['HGT']: 260 if var in ['HGT']: pass 241 261 elif colorb: 242 243 244 245 246 262 ndiv = 10 263 colorbar(fraction=0.05,pad=0.1,format=fmtvar(var),\ 264 ticks=np.linspace(zevmin,zevmax,ndiv+1),\ 265 extend='max',spacing='proportional') 266 # both min max neither 247 267 ### Vector plot 248 268 if winds: … … 257 277 else: colorvec = definecolorvec(colorb) 258 278 vectorfield(vecx, vecy,\ 259 x, y, stride=stride, csmooth= stride,\279 x, y, stride=stride, csmooth=2,\ 260 280 scale=15., factor=300., color=colorvec, key=key) 261 #200. ## or csmooth= 2281 #200. ## or csmooth=stride 262 282 263 283 ### Next subplot 264 284 plottitle = basename 265 if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) 266 else: plottitle = plottitle + "_LT"+str(ltst) 285 if typefile in ['mesoapi','meso']: 286 if addchar: plottitle = plottitle + addchar + "_LT"+str(ltst) 287 else: plottitle = plottitle + "_LT"+str(ltst) 267 288 ptitle( plottitle ) 268 289 sub += 1 … … 297 318 from netCDF4 import Dataset 298 319 from myplot import getlschar 320 from numpy import ones 299 321 300 322 ############################# 301 323 ### Get options and variables 302 324 parser = OptionParser() 303 parser.add_option('-f', action=' store', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file')325 parser.add_option('-f', action='append', dest='namefile', type="string", default=None, help='[NEEDED] name of WRF file (append)') 304 326 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)') 305 327 parser.add_option('-p', action='store', dest='proj', type="string", default=None, help='projection') … … 307 329 parser.add_option('-t', action='store', dest='target', type="string", default=None, help='destination folder') 308 330 parser.add_option('-s', action='store', dest='stride', type="int", default=3, help='stride vectors (def=3)') 309 parser.add_option('-v', action=' store', dest='var', type="string", default=None, help='variable contoured')310 parser.add_option('-n', action='store', dest='numplot', type="int", default= 4, help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)')331 parser.add_option('-v', action='append', dest='var', type="string", default=None, help='variable contoured (append)') 332 parser.add_option('-n', action='store', dest='numplot', type="int", default=2, help='number of plots (def=1)(<0: 1 plot of LT -*numplot*)') 311 333 parser.add_option('-i', action='store', dest='interp', type="int", default=None, help='interpolation method (2: press, 3: z-amr, 4:z-als)') 312 334 parser.add_option('-c', action='store', dest='colorb', type="string", default=True, help='change colormap') 313 335 parser.add_option('-x', action='store_false', dest='winds', default=True, help='no wind vectors') 314 parser.add_option('-m', action=' store', dest='vmin', type="float", default=None, help='bounding minimum value for color plot')315 parser.add_option('-M', action=' store', dest='vmax', type="float", default=None, help='bounding maximum value for color plot')336 parser.add_option('-m', action='append', dest='vmin', type="float", default=None, help='bounding minimum value (append)') 337 parser.add_option('-M', action='append', dest='vmax', type="float", default=None, help='bounding maximum value (append)') 316 338 parser.add_option('-T', action='store_true', dest='tile', default=False, help='draw a tiled plot (no blank zone)') 317 339 parser.add_option('-z', action='store', dest='zoom', type="float", default=None, help='zoom factor in %') … … 323 345 exit() 324 346 print "Options:", opt 325 zefile = opt.namefile 326 zelevel = opt.nvert 327 stralt = None 328 [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) 329 330 ##################################################### 331 ### Call Fortran routines for vertical interpolations 332 if opt.interp is not None: 333 if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 334 ### winds or no winds 335 if opt.winds : zefields = 'uvmet' 336 else : zefields = '' 337 ### var or no var 338 if opt.var is None : pass 339 elif zefields == '' : zefields = opt.var 340 else : zefields = zefields + "," + opt.var 341 print zefields 342 zefile = api_onelevel ( path_to_input = '', \ 343 input_name = zefile, \ 344 path_to_output = opt.target, \ 345 fields = zefields, \ 346 interp_method = opt.interp, \ 347 onelevel = zelevel, \ 348 nocall = opt.nocall ) 349 print zefile 350 zelevel = 0 ## so that zelevel could play again the role of nvert 351 352 ############# 353 ### Main call 354 name = winds (zefile,int(zelevel),\ 355 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ 356 addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom) 357 print 'Done: '+name 358 359 ######################################################### 360 ### Generate a .sh file with the used command saved in it 361 command = "" 362 for arg in sys.argv: command = command + arg + ' ' 363 f = open(name+'.sh', 'w') 364 f.write(command) 347 348 if opt.var is None: 349 pass 350 else: 351 listvar = '' 352 zelen = len(opt.var) 353 if zelen == 1: listvar = opt.var[0] + ',' 354 else : 355 for jjj in range(zelen): listvar += opt.var[jjj] + ',' 356 listvar = listvar[0:len(listvar)-1] 357 if opt.vmin: 358 if zelen != len(opt.vmin): 359 print "not enough or too much vmin values... setting same values all variables" 360 vmintab = ones(zelen) * opt.vmin[0] 361 else: 362 vmintab = opt.vmin 363 else: 364 vmintab = None 365 if opt.vmax: 366 if zelen != len(opt.vmax): 367 print "not enough or too much vmax values... setting same values all variables" 368 vmaxtab = ones(zelen) * opt.vmax[0] 369 else: 370 vmaxtab = opt.vmax 371 else: 372 vmaxtab = None 373 374 for i in range(len(opt.namefile)): 375 376 zefile = opt.namefile[i] 377 print zefile 378 zelevel = opt.nvert 379 stralt = None 380 [lschar,zehour,zehourin] = getlschar ( zefile ) ## getlschar from wrfout (or simply return "" if another file) 381 382 ##################################################### 383 ### Call Fortran routines for vertical interpolations 384 if opt.interp is not None: 385 if opt.nvert is 0 and opt.interp is 4: zelevel = 0.010 386 ### winds or no winds 387 if opt.winds : zefields = 'uvmet' 388 else : zefields = '' 389 ### var or no var 390 if opt.var is None : pass 391 elif zefields == '' : zefields = listvar 392 else : zefields = zefields + "," + listvar 393 print zefields 394 zefile = api_onelevel ( path_to_input = '', \ 395 input_name = zefile, \ 396 #path_to_output = opt.target, \ 397 fields = zefields, \ 398 interp_method = opt.interp, \ 399 onelevel = zelevel, \ 400 nocall = opt.nocall ) 401 print zefile 402 zelevel = 0 ## so that zelevel could play again the role of nvert 403 404 if opt.var is None: 405 ############# 406 ### Main call 407 name = winds (zefile,int(zelevel),\ 408 proj=opt.proj,back=opt.back,target=opt.target,stride=opt.stride,var=opt.var,numplot=opt.numplot,colorb=opt.colorb,winds=opt.winds,\ 409 addchar=lschar,interv=[zehour,zehourin],vmin=opt.vmin,vmax=opt.vmax,tile=opt.tile,zoom=opt.zoom) 410 print 'Done: '+name 411 else: 412 for jjj in range(len(opt.var)): 413 if vmintab: argvmin = vmintab[jjj] 414 else: argvmin = None 415 if vmaxtab: argvmax = vmaxtab[jjj] 416 else: argvmax = None 417 ############# 418 ### Main call 419 name = winds (zefile,int(zelevel),\ 420 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,\ 421 addchar=lschar,interv=[zehour,zehourin],vmin=argvmin,vmax=argvmax,tile=opt.tile,zoom=opt.zoom) 422 print 'Done: '+name 423 424 ######################################################### 425 ### Generate a .sh file with the used command saved in it 426 command = "" 427 for arg in sys.argv: command = command + arg + ' ' 428 f = open(name+'.sh', 'w') 429 f.write(command)
Note: See TracChangeset
for help on using the changeset viewer.
