;MIR ADT CEN - 10 March 2009. ;Run by starting IDL - type: ; idl ;Then (once in IDL) type: ; .run simp2.pro ;Simple IDL program to read in wrfout files and plot some results. ;simp1.pro focused on zonal means of 3D [lon-lat-alt] variables ;simp2.pro focuses on 2D surface variables [lon-lat] ;Intended to be used mainly by those wishing to validate their ;installation of the PlanetWRF code with the default Mars test run. ;Also provides some insight into what variables are available and how ;to manipulate them to obtain some very basic diagnostics ;----------------------------------------------------------------------------- ;INPUTS ;1. Open the wrfout file(s) print,'How many wrfout files do you have so far? read,nfiles tcum=0 for n=0,nfiles-1 do begin if (n eq 0) then ncid=ncdf_open('wrfout_d01_0001-00001_00:00:00',/nowrite) if (n eq 1) then ncid=ncdf_open('wrfout_d01_0001-00101_00:00:00',/nowrite) if (n eq 2) then ncid=ncdf_open('wrfout_d01_0001-00201_00:00:00',/nowrite) if (n eq 3) then ncid=ncdf_open('wrfout_d01_0001-00301_00:00:00',/nowrite) if (n eq 4) then ncid=ncdf_open('wrfout_d01_0001-00401_00:00:00',/nowrite) if (n eq 5) then ncid=ncdf_open('wrfout_d01_0001-00501_00:00:00',/nowrite) if (n eq 6) then ncid=ncdf_open('wrfout_d01_0001-00601_00:00:00',/nowrite) ;2. Read in grid dimensions (from the first file only): if (n eq 0) then begin print,'Reading in constant data' ;Latitude and longitude in deg ncdf_varget,ncid,'XLAT',xlat lat=reform(xlat(0,*)) ncdf_varget,ncid,'XLONG',xlong lon=reform(xlong(*,0)) ;3. Read in some planetary constants (from the first file only): ;Reference potential temperature (in K) [single value, constant] ncdf_attget,ncid,/global,'T0',t0 ;Reference pressure (in Pa) [single value, constant] ncdf_attget,ncid,/global,'P0',p0 ;Specific heat capacity (in J/kg/K) at constant pressure ncdf_attget,ncid,/global,'CP',cp ;Specific gas constant (in J/kg/K) [single value, constant] ncdf_attget,ncid,/global,'R_D',rdas ;Gravity (in m/s2) [single value, constant] ncdf_attget,ncid,/global,'G',grav ;Radius of body (in m) [single value, constant] ncdf_attget,ncid,/global,'RADIUS',radius ;Find some grid dimensions and set up some variables: a=size(xlat) imax=a(1) jmax=a(2) a=size(sigma) kmax=a(1) tstot=fltarr(imax,jmax,400*nfiles) ;400 outputs per file pstot=fltarr(imax,jmax,400*nfiles) co2tot=fltarr(imax,jmax,400*nfiles) lstot=fltarr(400*nfiles) endif ;4. Read in 2D variables that change with time print,'Reading in data for file ',n+1 ;Surface pressure (in Pa) ncdf_varget,ncid,'PSFC',ps ;Surface temperature (in K) ncdf_varget,ncid,'TSK',ts ;Surface CO2 ice cover (in kg) ncdf_varget,ncid,'CO2ICE',co2 ;Planetocentric solar longitude (in degrees) giving season of year: ;Ls=0 => northern spring equinox ;Ls=90 => northern summer solstice ;Ls=180 => northern autumn equinox ;Ls=270 => northern winter solstice ncdf_varget,ncid,'L_S',ls ;Find number of times and set up some variables: a=size(ls) tmax=a(1) ;5. Store 2D variables in longer arrays for t=0,tmax-1 do begin print,'Storing variables for time ',t,' file ',n+1 lstot(tcum)=ls(t) pstot(*,*,tcum)=ps(*,*,t) tstot(*,*,tcum)=ts(*,*,t) co2tot(*,*,tcum)=co2(*,*,t) tcum=tcum+1 endfor ;6. Close the wrfout file ncdf_close,ncid endfor ;----------------------------------------------------------------------------- ;DIAGNOSTICS ;7. Output some information lstot2=fltarr(tcum) pstot2=fltarr(imax,jmax,tcum) tstot2=fltarr(imax,jmax,tcum) co2tot2=fltarr(imax,jmax,tcum) lstot2(0:tcum-1)=lstot(0:tcum-1) pstot2(*,*,0:tcum-1)=pstot(*,*,0:tcum-1) tstot2(*,*,0:tcum-1)=tstot(*,*,0:tcum-1) co2tot2(*,*,0:tcum-1)=co2tot(*,*,0:tcum-1) print,'Your run has longitudes (in deg) ',lon(0),' to ',lon(imax-1) print,'Your run has latitudes (in deg) ',lat(0),' to ',lat(jmax-1) print,'Your run has reached Ls=',lstot2(tcum-1) ;----------------------------------------------------------------------------- ;PLOTTING ;8. Plot surface pressure, temperature and CO2 ice cover for several times of ;year, beginning with the initial conditions. print,'Using window command to set window size [setting xsize and ysize]' print,'May need to edit these in simp2.pro for your machine' window,xsize=800, ysize=1000 ;Set up some things to make the plots look better: print,'NB: using device, decomposed=0 - may not be needed' print,'If colors are bad, comment it out' device,decomposed=0 loadct,39 ;Set up reasonable contour levels for the whole spin-up year: plevs=findgen(18)*75.+150. tlevs=findgen(19)*10.+130. clevs=[0.,0.1,findgen(12)*150.+150.] ;Instruct IDL to label every other contour labs=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0] labs2=[0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0] ;If a time is not yet available (because the run has not progressed ;that far) then print a message to this effect. ;INITIAL CONDITIONS: !p.multi=[0,1,3] tls=1 string1=string('Surface pressure [Pa], Ls =',+lstot2(tls)) string2=string('Surface temperature [K], Ls =',+lstot2(tls)) string3=string('Surface CO2 ice cover [kg], Ls =',+lstot2(tls)) contour,pstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string1,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=plevs,/fill,charsize=2 contour,pstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string1,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=plevs,/over,c_labels=labs,charsize=2 contour,tstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string2,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=tlevs,/fill,charsize=2 contour,tstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string2,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=tlevs,/over,c_labels=labs,charsize=2 contour,co2tot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string3,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=clevs,/fill,charsize=2 contour,co2tot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string3,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=clevs,/over,c_labels=labs2,charsize=2 print,'Hit 1 to continue' read,hit ;LS = 10, 20, 30, 60, 90, 180, 270, 355 tim=8 lstim=fltarr(tim) lstim(0)=10. lstim(1)=20. lstim(2)=30. lstim(3)=60. lstim(4)=90. lstim(5)=180. lstim(6)=270. lstim(7)=355. for t=0,tim-1 do begin if (lstot2(tcum-1) lt lstim(t)) then begin print,'Your run must reach Ls=',lstim(t),' before I can plot anything else' stop endif else begin tls=0 while (lstot2(tls) lt lstim(t)) do tls=tls+1 endelse !p.multi=[0,1,3] string1=string('Surface pressure [Pa], Ls =',+lstot2(tls)) string2=string('Surface temperature [K], Ls =',+lstot2(tls)) string3=string('Surface CO2 ice cover [kg], Ls =',+lstot2(tls)) contour,pstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string1,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=plevs,/fill,charsize=2 contour,pstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string1,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=plevs,/over,c_labels=labs,charsize=2 contour,tstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string2,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=tlevs,/fill,charsize=2 contour,tstot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string2,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=tlevs,/over,c_labels=labs,charsize=2 contour,co2tot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string3,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=clevs,/fill,charsize=2 contour,co2tot2(*,*,tls),lon,lat,xrange=[-180,180],xstyle=1,yrange=[-90,90],ystyle=1,title=string3,xtitle='Longitude (deg E)',ytitle='Latitude (deg N)',levels=clevs,/over,c_labels=labs2,charsize=2 print,'Hit 1 to continue' read,hit endfor end