;; gpcprn_dlysmanom_wrtofil.ncl
;; Dr. Sujata Mandke  6-11-2017 
;;;To write GPCP daily smoothed anomaly 
;;;for all 365 days of all years to nc file
;++++++++++++++++++++++++++++++++++++++++
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   
begin
                           ; time window 
  twStrt = 19960101        
  twLast = 20081231
;;-------------------------------------------------
;;GPCP daily rain anomalies are calculated as follows
;;;;;----------------------------------------------

   ymdStrt = 19960101                         ; start yyyymmdd
   ymdLast = 20081231                         ; last  
   yrStrt  = ymdStrt/10000
   yrLast  = ymdLast/10000
  
   nhar    = 4                                ; number of fourier comp      

   var     = "data"                            ; name of file                   
   vName   = "precip"                            ; name for plots

   diri    = "/iitm3/erpas-res/skm/data/cmp5/gfdl/mjoclivardata/"       ; input directory   
   fili    = "pregpcp19962008.daily.nc" ; input file 

;***********************************************************
; Read user specified time and create required yyyyddd                    
;***********************************************************
   f       = addfile (diri+fili , "r")                          

   time    = f->time                          ; all times on file
   ymd     = cd_calendar(time, -2)            ; yyyymmdd
;;;;
   iStrt   = ind(ymd.eq.ymdStrt)              ; index start
   iLast   = ind(ymd.eq.ymdLast)              ; index last 
   delete(time)
   delete(ymd)

;***********************************************************
; Read user specified time and create required yyyyddd                    
;***********************************************************
   timsuj    = f->time(iStrt:iLast)         ; time:units = "hours since"
   timnew    = cd_calendar(timsuj, 0)       ; type float 
   year    = floattointeger( timnew(:,0) )
   month   = floattointeger( timnew(:,1) )
   day     = floattointeger( timnew(:,2) ) 
   ddd     = day_of_year(year, month, day) 
   yyyyddd = year*1000 + ddd                  ; needed for input

;***********************************************************
; Read data: short2flt                                     
;*********************************************************** 
   x       =  short2flt( f->$var$(iStrt:iLast,:,:) ) ; convert to float 
   printVarSummary( x ) 

;***********************************************************
; Compute daily climatology: raw and then 'smoothed'  
;***********************************************************
   xClmDay = clmDayTLL(x, yyyyddd)     ; daily climatology at each grid point
   printVarSummary(xClmDay)            

;***********************************************************
; Compute smoothed daily climatology using 'nhar' harmonics
;***********************************************************
   xClmDay_sm = smthClmDayTLL(xClmDay, nhar)
   printVarSummary(xClmDay_sm)

;***********************************************************
; Compute daily anomalies using raw and smoothed climatologies
;***********************************************************
    xAnom      = calcDayAnomTLL (x, yyyyddd, xClmDay)     
    printVarSummary(xAnom)

    xAnom_sm   = calcDayAnomTLL (x, yyyyddd, xClmDay_sm)     
    xAnom_sm@long_name = "Anomalies from Smooth Daily Climatology"
    printVarSummary(xAnom_sm)

    delete( x )    ; no longer needed
;;;--------------------------------
;;;follows writing smoothed daily gpcp 
;;;rain anomaly  to output netdcf file
;************************************************
  diro="/moes/home/amin/cmip5/gfdl/mjoclivar_ncl/wavenofrq/data/"
  filo1   = "prgpcp19962008.daily.smanom.nc"
  system("/bin/rm -f "+diro+filo1)   ; rm pre-exist file, if any
   f1       = addfile (diro+filo1 , "c") 
      
     filAtt = True
     filAtt@title         = "GPCP rain daily smooth anomaly 1996-2008" 
     filAtt@source_file   = fili
     filAtt@creation_date = systemfunc("date")
     fileattdef( f1, filAtt )    ; copy file attributes 
;;;	 
   f1->gpcprnsman=xAnom_sm 

end