;;;;;;
;;; ZTS 25 Aug 2023: RFE2 subseted from Africa domain
;;; Africa domain is from 40S-40N, 20W to 55E
;;;;;;;;

begin

 startday = 1
 startmon = 7  ; start month for  onset/cessation
 endmon   = 9  ; last month of observation
 endday   = 30 
 SEASON   = "JAS" ;e. "JJA"
 YearS    = 1998  ;This is for observation year
 YearE    = 2020   ;This is for observation year
 latS     = -22.0 
 latN     = 9.0 
 lonW     = 132.0 
 lonE     = 205.0 
 obsT     = "CMORPH"
 resl     = 0.5
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;; END OF CHANGEABLE PARAMETERS
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 if (resl .eq. 0.5) then
   obsFiles  = (/"cmorph.v2.0.daily.cPAp5."/) ; 
 else if (resl .eq. 0.1) then
   obsFiles  = (/"cmorph.v2.0.daily.cPAp1."/) ; 
 end if
 end if

 if (lonW .lt. 0.) then
    lonW = lonW + 360.
 end if

 if (lonE .lt. 0.) then
    lonE = lonE + 360.
 end if

 obsOut    = (/"CMORPH"/)
 if (startmon .gt. endmon) then
   nextYr  = True
 else
   nextYr  = False
 end if

 datdir   = "./"
 obsdir   = "../precip/cmorph/"

 monthdef = (/"","Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep", \
                    "Oct","Nov","Dec"/)

 initMonth = monthdef(startmon)

 writeoutany  = True

 iobsStrt = 0
 iobsLast = 0

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; this assignment use lat-lon defined earlier near the top
 mlat  = toint((latN-latS)/resl + 1)
 mlon  = toint((lonE-lonW)/resl + 1)
 REGlat = fspan(latS, latN, mlat)
 REGlon = fspan(lonW, lonE, mlon)
 REGlat@units = "degrees_north"
 REGlon@units  = "degrees_east"
 REGlat!0      = "lat"
 REGlon!0      = "lon"
  print(mlat+ " " +mlon)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;   Define hindcast time variable to save computed values

 yrs   = ispan(1998,2024,1) 
 mm    = yrs
 mm(:) = startmon
 dd    = mm
 dd(:) = startday
 hh    = dd
 hh(:) = 00
 mn    = hh
 ss    = hh
 units = "days since 1981-01-01 00:00"

 TIME      = cd_inv_calendar(yrs,mm,dd,hh,mn,ss,units,0)
 nyears    = dimsizes(TIME)
 TIME!0    = "time"

 time   = cd_calendar(TIME,-2)
 time!0 = "time"
 time@units = "yyyymmdd"

 ntype       = dimsizes(obsFiles)
 otype       = ispan(0,ntype - 1,1)
 otype!0     = "otype"
 otype@obs1  = obsOut
; otype@obs2  = obsOut(1)

;;;;;;

 jlatStrt   = 0
 jlatLast   = mlat - 1
 ilonStrt   = 0
 ilonLast   = mlon - 1

;;;; save onset for all years (nyears can be beyond obs data)
;;;;  and all observation tipes
 
 precipAll     = new((/nyears, dimsizes(obsFiles),mlat, mlon /),float,-9999)
 precipAll!0   = "time" 
 precipAll!1   = "otype"
 precipAll!2   = "lat"
 precipAll!3   = "lon"

 precipAll&time  = TIME 
 precipAll&otype = otype
 precipAll&lat   = REGlat
 precipAll&lon   = REGlon

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Get data 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 do iobs = iobsStrt, iobsLast
    
    fname  := systemfunc("ls -1 " + obsdir + obsFiles(iobs) + "*" + ".nc")
      a        := addfiles(fname,"r")
      tmp      := a[:]->prec 
      rain_tot := tmp(:,{latS:latN},{lonW:lonE})
      otime    := a[:]->time
      delete(fname) 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;; to select all indices within a specified years
    oyymmdd    := cd_calendar(otime,-2)

   ; print(dimsizes(oyymmdd))
;;;;;;
    odate     := cd_calendar(otime, 0)   ;(YYYY,MM,DD,hr, mn, ss)
    oYrStart  := (/ toint(odate(0,0)) /)
    dims       = dimsizes(odate)
    oYrEnd    := (/ toint(odate(dims(0) - 1, 0)) /)
    oAllyrs   := ispan(oYrStart, oYrEnd, 1)

    yrStrt    := oAllyrs

    if(nextYr) then
      yrLast  := yrStrt + 1
      ndimyr   = dimsizes(yrStrt) - 2
    else
      yrLast  := yrStrt
      ndimyr   = dimsizes(yrStrt) - 1
    end if

;;;;;;;;;;;;;;;
;;;;;;; Loop over all years (initial time)
;;;;;;
    do iyr = 0, ndimyr
   
       iyears  = yrStrt(iyr)
       iyeare  = yrLast(iyr)
       yrind   = ind(yrs .eq. iyears ) 

       yyyy   :=  oyymmdd/10000
   
       ii     := ind(yyyy.ge.iyears  .and. yyyy.le.iyeare)
     ;print(iyears + " " + iyeare)
     ;print(yrind + " " + ii)
     ;print(yyyy)

;;;;; Times with in the data period ffor each year
       ytime  := otime(ii)  
       ydate  := cd_calendar(ytime, 0)  ;6 cols for yr,MM,DD,HH,mm,ss
       ymonth :=  tointeger(ydate(:,1)) ;all months of data for each day
       yday   := tointeger(ydate(:,2))  ;all days of data for a given year 

       timelabelsObs := monthdef(ymonth)+"-"+yday ; all times per year

;;;; for onset specfied start and end times are
       tmpS = monthdef(startmon)+ "-" + startday  ; onset search start day
       tmpE = monthdef(endmon)  + "-" + endday    ; onset search last day

      print(iyears+" with data dates from "+ timelabelsObs(0)+ " to  " +timelabelsObs(dimsizes(timelabelsObs)-1))
      print( " Data summed over " + tmpS +" to " + tmpE)

       indxOnsetStart = get1Dindex(timelabelsObs,tmpS)
       indxOnsetEnd   = get1Dindex(timelabelsObs,tmpE)

       timelabels := timelabelsObs(indxOnsetStart:indxOnsetEnd)
      ; print(timelabels)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Create a daily rainfall variable; set it to 0 if less than 0.01mm
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

       raindaily := rain_tot(ii,:,:)
       
       raindaily := raindaily(indxOnsetStart:indxOnsetEnd,:,:)
       raindaily := dim_sum_n_Wrap(raindaily,0)

       raindaily = where(raindaily .lt. 0.1, 0.0, raindaily)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
          print("Yr: " + iyears )
       precipAll(yrind,iobs,:,:)        = (/ raindaily     /) 
    end do ; iyr    

 end do ; iobs
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 

  ;print(dimsizes(precipAll)) 
  fname := "./prec" + SEASON + "_" + YearS + "-" + YearE + obsT + ".nc"
  system("/bin/rm -f "+fname)
  ncdf = addfile(fname, "c")

       fAtt               = True    ; assign file attributes
       fAtt@title         = "Obtained from continental RFE2 data"
       fAtt@source_file   =  "RFE2"
       fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

       filedimdef(ncdf,"time",-1,True)

       ncdf->prec  = precipAll

end