;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; LEC_AGG_NETCDF ; This program has to be run with Pvwave. ; It reads the netcdf files given for the Rhone Aggreg experiment, ; and plot 1D and 2D instantaneous or averaged fields. ; User defined variables have to be modified in section 0. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; 0. USER DEFINED PARAMETERS ; These Variables should be modified ;************************************* ; Folder where to find the Rhone-Agg masks file: REPMASK='/data0/mc2/lemoigne/boone/rhoneagg/Sfc_dat/Exp1/' ; Folder where to find the file to plot REP='/data0/mc2/lemoigne/boone/rhoneagg/sfc_dat/Exp1/' ; Name of the file to plot file='rhone_agg_exp1_veg.nc' ; Nombers of color (or class) used to plot NB_CLASS=10 ; Flag used only if atmospheric forcing is plotted AVG_FLAG=1 ; 1: will plot time averaged value ; 0: will plot instantaneous values ; 1. modify the wave_startup : check for HDF Interface ; This is not necessary for my platform, but maybe you'll need it ;****************************************************************** ;if (N_elements(deb) eq 0) then deb=1 ; IF ((deb ne 0) or STRMATCH ( !Path, 'hdf-|HDF-' ) NE 0 ) THEN BEGIN & $ ; DEFSYSV, "!Hdf_debug", 0 & $ ; @hdf_common & $ ; HDF_INIT & $ ; LOAD_OPTION, 'hdf', 3.3 & $ ; IF !Err EQ 0 THEN PRINT, " PV-WAVE:HDF 3.30 Module Initialized" & $ ; confirmspace = 1 & $ ; ENDIF ;deb=0 flag=-999. zinf=[0.]/[0.] ;2. READ MASKS ;************* ; Exept for experience 1, Two masks are needed. ; The program read the corresponding mask according ; to the number of the experiment ; 2.1 Read the 8x8 km2 mask ;========================== file_mask=REPMASK+'rhone_agg_exp1_masks.nc' iloop=2 ncid=ncopen(file_mask,nc_nowrite) status=ncinquire(ncid,ndims,nvars,natts,recdim) repeat begin i=ncvarinq(ncid,iloop,titre,datatype,ndims,dim,natts) start =[0] count =[1471] ztab=lonarr(count(0)) j=ncvarget(ncid,iloop,start,count,ztab) if (iloop eq 2) then maskx=ztab else masky=ztab iloop=iloop+1 endrep until (iloop eq 4) ; 2.2 Read the mask associated to the number of the experiment ;============================================================= case 1 of (strpos(file,'exp1' ) ge 0) : exp=1 (strpos(file,'exp2a') ge 0) : exp=2 (strpos(file,'exp2b') ge 0) : exp=3 (strpos(file,'exp2c') ge 0) : exp=2 (strpos(file,'exp3' ) ge 0) : exp=4 endcase print,'experience is ',exp if (exp eq 1) then begin mask_exp=indgen(1471) endif else begin iloop=exp+2 i=ncvarinq(ncid,iloop,titre,datatype,ndims,dim,natts) print,titre start =[0] count =[1471] ztab=lonarr(count(0)) j=ncvarget(ncid,iloop,start,count,ztab) mask_exp=ztab-1 endelse iret=ncclose(ncid) INBTOT=max(mask_exp+1) IXMAX=max(maskx) IYMAX=max(masky) ; 3 initiliaze graphic properties ;******************************* yn='' device,/close set_plot,'x' loadct,4 ; 4 OPEN netCDF file ;******************************* filename=REP+file ncid=ncopen(filename,nc_nowrite) print,'file is ',filename,' nid is ',ncid ; 5 INQUIRE numbers of parameters ;********************************* status=ncinquire(ncid,ndims,nvars,natts,recdim) ; 6 READ EACH VARIABLE ;********************* !p.multi=[0,1,1] iloop=0 repeat begin ; 6.1 inquire name, numbers of dimension, and type of data of the variable ;========================================================= i=ncvarinq(ncid,iloop,titre,datatype,ndims,dim,natts) name=titre ; 6.2 inquire dimension of the variable ;========================================================= start =intarr(ndims) & start(*)=0 length=intarr(ndims) sdim =strarr(ndims) jdim=0 ndtot=1 repeat begin diminq=dim(jdim) j=ncdiminq(ncid,diminq,titre,len) length(jdim)=len sdim (jdim)=titre jdim=jdim+1 ndtot=ndtot*len endrep until (jdim eq ndims) count=length ; 6.3 define variable size, starting point, and numbers of data to read ;===================================================================== s=[ndims,length,datatype-1,ndtot] ztab=make_array(size=s,value=datatype) ; 6.4 read variable ;=================== j=ncvarget(ncid,iloop,start,count,ztab) ; 6.5 Plot ;=========== if ((strpos(name,'time') lt 0) and (max(s) ge INBTOT)) then begin zori=ztab nbfield=1 i=where(sdim eq 'tstep',ic) if ((ic ne 0) and (AVG_FLAG eq 1)) then begin ztab=sum(zori,i) zle=length(i)*1.0 ztab(*)=ztab(*)/zle(0) nbfield=1 endif if ((ic ne 0) and (AVG_FLAG eq 0)) then begin nbfield=length(i(0)) endif if (nbfield gt 1) then begin icol=min([3,nbfield/4.]) irow=min([4,nbfield/3.]) !p.multi=[0,icol,irow] endif for ifield=0,min([100,nbfield-1]) do begin s=[2,IXMAX,IYMAX,datatype-1,2] ztab_rhone=make_array(size=s,value=datatype) ztab_rhone(*)=flag indflag=where(mask_exp ge 0) if (nbfield eq 1) then ztab_rhone(maskx(indflag)-1,masky(indflag)-1)=ztab(mask_exp(indflag)) $ else ztab_rhone(maskx(indflag)-1,masky(indflag)-1)=ztab(mask_exp(indflag),ifield) imask=where(ztab_rhone ne flag,icf) fo='(f6.1)' if( icf gt 0) then begin ext=' n='+string(min(ztab_rhone(imask)),form=fo)+ $ ' x='+string(max(ztab_rhone(imask)),form=fo)+ $ ' a='+string(avg(ztab_rhone(imask)),form=fo) if(nbfield gt 0) then ext=ext+'t='+string(ifield,form=fo) fact=500./max([s]) ztabplot=congrid(ztab_rhone,s(1)*fact,s(2)*fact) mini=min(ztab) maxi=max(ztab) tlevels=findgen(NB_CLASS)*(maxi-mini)/(NB_CLASS-1.)+mini tcolor =findgen(NB_CLASS+1)*250./(NB_CLASS) if(strpos(!d.name,'X') ge 0 ) then window,/free,xs=s(1)*fact,ys=s(2)*fact,title=name+ext contour,ztabplot,levels=tlevels,c_color=tcolor,C_charsize=1,Path='path.dat',/follow contourfill,'path.dat',ztabplot,color_index=tcolor contour,ztabplot,levels=tlevels,C_charsize=1.5,/noerase,/follow,title=name+ext endif endfor endif else begin if (strpos(name,'Dis') lt 0) then begin plot,ztab,title=name;,xr=xr,yr=yr endif else begin i=where(sdim eq 'x',ic) zori=ztab plot,ztab,/nodata,title='streamflow' for ins=0,ic do begin zdis=ztab(ins,*) oplot,ztab(ins,*),thick=2,linestyle=ins endfor endelse endelse iloop=iloop+1 endrep until (iloop eq nvars) ; 7 CLOSE netCDF file ;******************** iret=ncclose(ncid) end