; Get GPC GCM  daily predicted rainfall from C3S
;ZTS: Jan 2022: This code reads forecast grib & write out netcdf 
;;;;  CFS and eccc  data from C3S have issues to write to netCDF.
;;;; 30 Apr 2022: For eccc (likely also for GloSea) use grib_to_netcdf
;;;; grib_to_netcdf -D NC_FLOAT -R 19900101 -T -o out.nc ecccDailyHindcast.grib
;;;;  This will solve conflict of writing two identical dates, which would fail 
;;;;  because netcdf does not allow writing values with two identical times
;;;; https://confluence.ecmwf.int/display/ECC/GRIB+tools+examples 
   
begin

 TYPE      = "Forecast" 

 modelNames := (/"eccc" /)
 version    := (/3/)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;REG boundaries

 latS           = -20.0
 latN           = 25.0
 lonW           = -180.0
 lonE           = 179.0

 fileyr             = CURRENTYR
 iMon               = OMONS
 latS               = LATS
 latN               = LATN
 lonW               = LONW
 lonE               = LONE

 Res                = RESL
 mlat  = toint((latN-latS)/Res + 1)
 mlon  = toint((lonE-lonW)/Res + 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"
 REGlat&lat    = REGlat
 REGlon&lon    = REGlon

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

 filemonth = allmonths(iMon)

 
 imodel = 0
 fName = "./"+modelNames(imodel)+ "Daily"+ filemonth + fileyr + TYPE + ".grib"

 system(" grib_to_netcdf -D NC_FLOAT -R 19900101 -T -o output.nc " + fName )
 a  = addfile("output.nc","r")
    ; output is: tp(date, step, number, latitude, longitude)
    ; or       : tp(step, number, latitude, longitude)
 x       := a->tp
 ndims    = dimsizes(x)
 gens := a->number
 glead := a->step

 if (any(isnan_ieee(x))) then
     if(.not.isatt(x,"_FillValue")) then
        x@_FillValue = default_fillvalue(typeof(x))
     end if
     replace_ieeenan (x, x@_FillValue, 0)
 end if

 gprec     = x
 delete(x)


 glat  := a->latitude
 glon  := a->longitude
 glat  := glat(::-1)

 nlat   = dimsizes(glat)
 nlon   = dimsizes(glon)
 nens   = dimsizes(gens)

 vNames := getfilevarnames (a)
 nNames := dimsizes(vNames)

 ii = str_match_ind(vNames,"date")
    
 if (ismissing(ii)) then
  ;; No date dimension only step
         
    gunits = glead@units
    xstr   = str_split(gunits," ")     
    xstr  := xstr(2)   ; this should contain YYYY-MM-DD
    gdate  = str_split(xstr,"-")
    iyear  = toint(gdate(0))
    imonth = toint(gdate(1))
    iday   = toint(gdate(2))

    TIME   := cd_inv_calendar(iyear,imonth,iday,00,00,00,\
                              "days since 1990-01-01 00:00:0",0)
    TIME!0  = "time"

    time       := cd_calendar(TIME,-3)
    time!0     = "time"
    time@units = "yyyymmddhh"

    initime    = time

    nyears     = dimsizes(time)
    nlead      = dimsizes(glead)
 
    y     := gprec(step|:,\
                    number|:,latitude|::-1,longitude|:)
    x := y
    x = (/ y * 1000. /) ; daily precip is in m; change to mm
    x@units  =    "mm/day"
    x(1:(nlead-1),:,:,:) = x(1:(nlead-1),:,:,:) - x(0:(nlead-2),:,:,:)
    delete(y)

    y  = new((/nyears,nlead,nens,nlat,nlon/),typeof(gprec),gprec@_FillValue)
    y(0:(nyears-1),:,:,:,:) = (/x /)
    delete(x)

 else

    xtime   := a->date
    TIME    := cd_convert(xtime,"days since 1990-01-01 00:00:0")
    TIME!0  = "time"

    time      := cd_calendar(TIME,-3)
    time!0     = "time"
    time@units = "yyyymmddhh"

    initime    := time
    nyears      = dimsizes(time)
    nlead       = dimsizes(glead) 
    


    y     := gprec(date|:,step|:,\ 
                   number|:,latitude|::-1,longitude|:)
    x := y
    x = (/ y * 1000. /) ; daily precip is in m; change to mm
    x@units  =    "mm/day"
    x(:,1:(nlead-1),:,:,:) = x(:,1:(nlead-1),:,:,:) - x(:,0:(nlead-2),:,:,:)
    y    := x
    delete(x)
 end if

 y!0    = "time"
 y!1    = "lead"
 y!2    = "ens"
 y!3    = "lat"
 y!4    = "lon"
 y&time = time

 glead           = glead/24
 glead!0         = "lead"
 glead&lead      = glead
 glead@units     = "days"
 glead@long_name = "Forecast offset from initial time"

 y@_FillValue    = gprec@_FillValue
 y@description   = "Total Forecast Precipitation "
 y&lead = glead
 y&ens  = gens
 y&lat  = glat
 y&lon  = glon

 delete(gprec)

 ens    := ispan(1,nens, 1) 
 lead   := glead  ;ispan(1,nlead,1)  ; Just like data downloaded from IRI,
                         ; except leadtimes are expressed as ceil(0.5,1.5,..)
 ens!0           = "ens"
 ens@units       = "ensemble members id"
 lead!0          = "lead"
 lead@units      = "days"
 lead@long_name  = "Forecast offset from initial time"

 x := lonFlip(y)
 y := x
 glon = (/ y&lon /)

 prec :=  linint2_Wrap(glon,glat,y, False, REGlon,REGlat, 0)
 delete(y)

 prec!0    = "time"
 prec!1    = "lead"
 prec!2    = "ens"
 prec!3    = "lat"
 prec!4    = "lon"

 prec&time = time
 prec&lead       = lead
 prec&ens        = ens
 prec&lat        = REGlat
 prec&lon        = REGlon
 prec@_FillValue = -9999

 printVarSummary(prec)
 delete(a)

;;;;;;;;;;;;;;;
 
 tmp  := tochar(tostring(Res))
 dir0 = "./regionalp" + tmp(2) + "/"
 system("if ! test -d " + dir0 +" ; then mkdir -p " + dir0 + " ; fi")

 MODEL = str_upper(modelNames(imodel))+ "v" +version(imodel)

 fname = dir0+ "precDaily"+ filemonth + fileyr + MODEL +".nc"

 system("/bin/rm -f " + fname) 
 fout  = addfile(fname ,"c")  ; netCDF out file

;***************************************************
; define file options [Version a033]
;***********************************************
 setfileoption(fout, "prefill", False)
 setfileoption(fout, "suppressclose", True)
 setfileoption(fout, "definemode", True)

;===================================================================
; create global attributes of the file (optional)
;===================================================================
 fAtt               = True            ; assign file attributes
 fAtt@title         = "Regridded Global Data from C3S"
 fAtt@source_file   = TYPE + " from CS3 portal for " + modelNames(imodel)
 fAtt@Conventions   = "None"
 fAtt@creation_date = systemfunc ("date")
 fileattdef( fout, fAtt )            ; copy file attributes

;***********************************************
; predefine coordinate information                
;***********************************************
 dimNames = (/"time", "lead","ens", "lat", "lon" /)
 dimSizes = (/ nyears ,nlead, nens, mlat,  mlon /)
 dimUnlim = (/ True ,False, False, False,  False /)
 filedimdef(fout, dimNames  , dimSizes,  dimUnlim )
 filevardef(fout, "time"  , typeof(time), getvardims(time) )
 filevarattdef(fout, "time", time)
 filevardef(fout, "lead" , typeof(lead), getvardims(lead))
 filevarattdef(fout, "lead", lead)
 filevardef(fout, "ens"  , typeof(ens) , getvardims(ens) )
 filevarattdef(fout, "ens", ens)
 filevardef(fout, "lat", typeof(REGlat), getvardims(REGlat) )
 filevarattdef(fout, "lat", REGlat)
 filevardef(fout, "lon", typeof(REGlon), getvardims(REGlon) )
 filevarattdef(fout, "lon", REGlon)

;***********************************************
; predefine variable sizes
;*********************************************** 
 filevardef(fout,"TIME",typeof(TIME), getvardims(TIME) )
 filevarattdef(fout,"TIME",TIME)

 filevardef(fout,"initime",typeof(initime), getvardims(initime) )
 filevarattdef(fout,"initime",initime)

 filevardef(fout,"prec",typeof(prec), getvardims(prec) )
 filevarattdef(fout,"prec",prec)

;***********************************************
; terminate define mode: not necessary but for clarity
;*********************************************** 
 setfileoption(fout, "definemode", False)

;===================================================================
; write data values to predefined locations  
; (/ .../) operator transfer values only        
;***********************************************
 fout->time   = (/ time /)
 fout->lead   = (/ lead /)
 fout->ens    = (/ ens /)
 fout->lat    = (/ REGlat /)
 fout->lon    = (/ REGlon /)

 fout->TIME    = (/ TIME /)
 fout->initime = (/ initime /)
 fout->prec    = (/ prec /)

 delete(prec)
 delete(time)
 delete(lead)
 delete(ens)
 delete(glat)
 delete(glon)

 delete(fout)

;;;;;;;;;;;;;
;;;;  do now for hindcast
;;;;;;;;;;;;;

 TYPE       = "Hindcast"

 glatS = latS - 1
 glatN = latN + 1
 glonW = lonW - 1
 glonE = lonE + 1

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 fName = "./"+modelNames(imodel)+ "Daily"+ filemonth + TYPE + "*.grib"

 system(" grib_to_netcdf -D NC_FLOAT -R 19900101 -T -o output.nc " + fName )
 a  := addfile("output.nc","r")
    ; output is: tp(date, step, number, latitude, longitude)
    ; or       : tp(step, number, latitude, longitude)
 x       := a->tp
 ndims    := dimsizes(x)
 gens := a->number
 glead := a->step

 if (any(isnan_ieee(x))) then
     if(.not.isatt(x,"_FillValue")) then
        x@_FillValue = default_fillvalue(typeof(x))
     end if
     replace_ieeenan (x, x@_FillValue, 0)
 end if

 gprec     = x
 delete(x)

 glat  := a->latitude
 glon  := a->longitude
 glat  := glat(::-1)

 nlat   = dimsizes(glat)
 nlon   = dimsizes(glon)
 nens   = dimsizes(gens)

 vNames := getfilevarnames (a)
 nNames := dimsizes(vNames)

 ii = str_match_ind(vNames,"date")
    
 if (ismissing(ii)) then
  ;; No date dimension only step
         
    gunits = glead@units
    xstr   = str_split(gunits," ")     
    xstr  := xstr(2)   ; this should contain YYYY-MM-DD
    gdate  = str_split(xstr,"-")
    iyear  = toint(gdate(0))
    imonth = toint(gdate(1))
    iday   = toint(gdate(2))

    TIME   := cd_inv_calendar(iyear,imonth,iday,00,00,00,\
                              "days since 1990-01-01 00:00:0",0)
    TIME!0  = "time"

    time       := cd_calendar(TIME,-3)
    time!0     = "time"
    time@units = "yyyymmddhh"

    initime    = time

    nyears     = dimsizes(time)
    nlead      = dimsizes(glead)
 
    y     := gprec(step|:,\
                    number|:,latitude|::-1,longitude|:)
    x := y
    x = (/ y * 1000. /) ; daily precip is in m; change to mm
    x@units  =    "mm/day"
    x(1:(nlead-1),:,:,:) = x(1:(nlead-1),:,:,:) - x(0:(nlead-2),:,:,:)
    delete(y)

    y  = new((/nyears,nlead,nens,nlat,nlon/),typeof(gprec),gprec@_FillValue)
    y(0:(nyears-1),:,:,:,:) = (/x /)
    delete(x)

 else

    xtime   := a->date
    TIME    := cd_convert(xtime,"days since 1990-01-01 00:00:0")
    TIME!0  = "time"

    time      := cd_calendar(TIME,-3)
    time!0     = "time"
    time@units = "yyyymmddhh"

    initime    := time
    nyears      = dimsizes(time)
    nlead       = dimsizes(glead) 
    


    y     := gprec(date|:,step|:,\ 
                   number|:,latitude|::-1,longitude|:)
    x := y
    x = (/ y * 1000. /) ; daily precip is in m; change to mm
    x@units  =    "mm/day"
    x(:,1:(nlead-1),:,:,:) = x(:,1:(nlead-1),:,:,:) - x(:,0:(nlead-2),:,:,:)
    y    := x
    delete(x)
 end if

 y!0    = "time"
 y!1    = "lead"
 y!2    = "ens"
 y!3    = "lat"
 y!4    = "lon"
 y&time = time

 glead           = glead/24
 glead!0         = "lead"
 glead&lead      = glead
 glead@units     = "days"
 glead@long_name = "Forecast offset from initial time"

 y@_FillValue    = gprec@_FillValue
 y@description   = "Total Forecast Precipitation "
 y&lead = glead
 y&ens  = gens
 y&lat  = glat
 y&lon  = glon

 delete(gprec)

 ;y := y(:,:,:,{glatS:glatN}, {glonW:glonE})
 ;glat  := glat({glatS:glatN})
 ;glon  := glon({glonW:glonE})


 ens    := ispan(1,nens, 1) 
 lead   := glead  ;ispan(1,nlead,1)  ; Just like data downloaded from IRI,
                         ; except leadtimes are expressed as ceil(0.5,1.5,..)
 ens!0           = "ens"
 ens@units       = "ensemble members id"
 lead!0          = "lead"
 lead@units      = "days"
 lead@long_name  = "Forecast offset from initial time"

 x    := lonFlip(y)
 y := x
 y := y(:,:,:,{glatS:glatN}, {glonW:glonE})
 glon :=  y&lon 
 glat := y&lat 
 ;nlat   := dimsizes(glat)
 ;nlon   := dimsizes(glon)

 prec :=  linint2_Wrap(glon,glat,y, False, REGlon,REGlat, 0)
 delete(y)

 prec!0    = "time"
 prec!1    = "lead"
 prec!2    = "ens"
 prec!3    = "lat"
 prec!4    = "lon"

 prec&time = time
 prec&lead       = lead
 prec&ens        = ens
 prec&lat        = REGlat
 prec&lon        = REGlon
 prec@_FillValue = -9999

 printVarSummary(prec)
 delete(a)

;;;;;;;;;;;;;;;
; Save data on local disk. The text file is for CPT pre-preparation
 
 tmp  := tochar(tostring(Res))
 dir0 = "./regionalp" + tmp(2) + "/"
 system("if ! test -d " + dir0 +" ; then mkdir -p " + dir0 + " ; fi")

 MODEL = str_upper(modelNames(imodel))+ "v" +version(imodel)

 fname = dir0 + "precDaily"+filemonth + TYPE + MODEL +".nc"

 system("/bin/rm -f " + fname) 
 fout  = addfile(fname ,"c")  ; netCDF out file

;***************************************************
; define file options [Version a033]
;***********************************************
 setfileoption(fout, "prefill", False)
 setfileoption(fout, "suppressclose", True)
 setfileoption(fout, "definemode", True)

;===================================================================
; create global attributes of the file (optional)
;===================================================================
 fAtt               = True            ; assign file attributes
 fAtt@title         = "Regridded Global Data from C3S"
 fAtt@source_file   = TYPE + " from CS3 portal for " + modelNames(imodel)
 fAtt@Conventions   = "None"
 fAtt@creation_date = systemfunc ("date")
 fileattdef( fout, fAtt )            ; copy file attributes

;***********************************************
; predefine coordinate information                
;***********************************************
 dimNames = (/"time", "lead","ens", "lat", "lon" /)
 dimSizes = (/ nyears ,nlead, nens, mlat,  mlon /)
 dimUnlim = (/ True ,False, False, False,  False /)
 filedimdef(fout, dimNames  , dimSizes,  dimUnlim )
 filevardef(fout, "time"  , typeof(time), getvardims(time) )
 filevarattdef(fout, "time", time)
 filevardef(fout, "lead" , typeof(lead), getvardims(lead))
 filevarattdef(fout, "lead", lead)
 filevardef(fout, "ens"  , typeof(ens) , getvardims(ens) )
 filevarattdef(fout, "ens", ens)
 filevardef(fout, "lat", typeof(REGlat), getvardims(REGlat) )
 filevarattdef(fout, "lat", REGlat)
 filevardef(fout, "lon", typeof(REGlon), getvardims(REGlon) )
 filevarattdef(fout, "lon", REGlon)

;***********************************************
; predefine variable sizes
;*********************************************** 
 filevardef(fout,"TIME",typeof(TIME), getvardims(TIME) )
 filevarattdef(fout,"TIME",TIME)

 filevardef(fout,"initime",typeof(initime), getvardims(initime) )
 filevarattdef(fout,"initime",initime)

 filevardef(fout,"prec",typeof(prec), getvardims(prec) )
 filevarattdef(fout,"prec",prec)

;***********************************************
; terminate define mode: not necessary but for clarity
;*********************************************** 
 setfileoption(fout, "definemode", False)

;===================================================================
; write data values to predefined locations  
; (/ .../) operator transfer values only        
;***********************************************
 fout->time   = (/ time /)
 fout->lead   = (/ lead /)
 fout->ens    = (/ ens /)
 fout->lat    = (/ REGlat /)
 fout->lon    = (/ REGlon /)

 fout->TIME    = (/ TIME /)
 fout->initime = (/ initime /)
 fout->prec    = (/ prec /)

 delete(prec)
 delete(time)
 delete(lead)
 delete(ens)
 delete(glat)
 delete(glon)

 delete(fout)

end