load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
;*****************************************************
; -------------   MAIN DRIVER ------------------------
;*****************************************************

;*****************************************************
; Miscellaneous initializations
;*****************************************************
setfileoption("nc", "SuppressClose", False)
	yrStrt  = 1990
        yrLast  = 2010
	yrLast  = yrStrt
	nmoStrt =    1
	nmoLast =    12
	var     = "R_GDS4_ISBL"
dRoot   = "/data1/data/idosam/Era_netcdf2/"   
fRoot   = "ei.oper.an.pl.regn128sc."                   ; file root    

   print("===================================================================")
   print("dRoot="+dRoot)
   print("fRoot="+fRoot)
   print("===================================================================")


;*****************************************************
; Missing information from the file
;*****************************************************
        lev = (/1,2,3,5,7,10,20,30,50,70,100                 \
            ,125,150,175,200,225,250,300,350,400,450,500  \
            ,550,600,650,700,750,775,800 ,825,850,875,900 \
            ,925,950,975,1000 /)
        lev!0 = "lev"
        lev@long_name = "pressure level"
        lev@units     = "hPa" 

        nlat = 256
        lat  = latGau(nlat, "lat", "latitude", "degrees_north")
        mlon = 2*nlat
        lon  = lonGlobeF(mlon, "lon", "longitude", "degrees_east")

;*****************************************************
; Read 4x-daily data
;*****************************************************

 	do year=yrStrt,yrLast
    	do nmo=nmoStrt,nmoLast
       	yyyymm = year*100+nmo

       	print(" ") 
       	print("=====> Calculating Monthly Mean: "+yyyymm+" <=====")
            
       	diri   = dRoot 
       	fili   = systemfunc("cd "+diri+" ; ls "+fRoot+yyyymm+"*")    
       	print(fili)
       	nfili  = dimsizes(fili)

       	f      = addfiles( diri+fili, "r")
	ListSetType(f, "join")
	x      = f[:]->$var$
       	printVarSummary(x)
        delete(x@initial_time)   ; not nneded for monthly mean
                                 ; create time coordinate information
        date  = toint(str_get_cols(fili, 24, 33))

        yr    = toint(str_get_cols(fili, 24, 27))
        mon   = toint(str_get_cols(fili, 28, 29))
        day   = toint(str_get_cols(fili, 30, 31))
        hour  = toint(str_get_cols(fili, 32, 33))
        mn    = new (dimsizes(hour), "integer")
        sc    = new (dimsizes(hour), "integer")
        mn    = 0
        sc    = 0

        units = "hours since 1901-1-1 00:00:0.0"
        time  = cd_inv_calendar(yr,mon,day,hour,mn,sc,units, 0)
        time!0= "time"
        time&time = time
        delete(time@_FillValue)
         
        x!0    = "time"    ; create time coordinate
        x&time =  time 
        x!1    = "lev"     ; rename to nicer dimension names
        x!2    = "lat"
        x!3    = "lon"
  
       	xMon   = dim_avg_n_Wrap(x, 0) 
       	printVarSummary(xMon)         ; [lev | 37] x [lat | 256] x [lon | 512]
     
                                      ; add a time dimension
       	dimx   = dimsizes(x)
       	rankx  = dimsizes(dimx)
       	print(dimx)
       	print("rank="+rankx)

       	dimMon = dimx
       	dimMon(0) = 1
       	print(dimMon)

       	if (rankx.eq.3) then
            xMON   = conform_dims( dimMon, xMon, (/1,2/) )
      	    copy_VarMeta(xMon, xMON(0,:,:))
	else
	    xMON   = conform_dims( dimMon, xMon, (/1,2,3/) )
	    copy_VarMeta(xMon, xMON(0,:,:,:))
       	end if
      ;;printVarSummary(xMON)

       	xMON!0    = "time"
       	xMON&time = x&$x!0$(0)       ; assign is day of month

        xMON!1    = "lev"
        xMON!2    = "lat"
        xMON!3    = "lon"
        xMON&lev  =  lev
        xMON&lat  =  lat
        xMON&lon  =  lon

      	printVarSummary(xMON)

        yyyymm    = date(0)/10000
        yyyymm!0  = "time"
        yyyymm@units= "yyyymm"
        printVarSummary(yyyymm)

;*****************************************************
; Write monthly means to netCDF
; pleas read: http://www.ncl.ucar.edu/Applications/method_1.shtml
;*****************************************************
        fout = "RHMON."+date(0)+".nc"
	print("Writing out netcdf: "+fout)
	system("/bin/rm -f "+fout)    

	ncdf = addfile(fout, "c")    
	filedimdef(ncdf,"time",-1, True)
        ncdf->date = yyyymm               ; only the 
	ncdf->RH   = xMON


     	delete([/ fili, x, xMon, xMON, dimx, dimMon /])
     	delete([/ yr, mon, day, hour, mn, sc, date, time /])
     
      end do  ; nmo
      end do      ; year