; Get GPC GCM  daily predicted rainfall from C3S
;ZTS: Oct 2023: This code is revised. It reads forecast grib 
;;;;   and regrids, processes and writes out netcdf. Now adjusted 
;;;;  for automation. 
;;;;
;;;;  
;;;;
;;;; These are place holders for fileyr to start at line# 25 (Do Contr + G)   
;;;; to see the line number. Don't delete the space (line# 10) below

begin

 TYPE       = "Forecast" 

 modelNames := (/"metfr"/)
 version    := (/8 /)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;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"

 
 a   = addfile(fName,"r")


 x    := a->TP_GDS0_SFC
 gens := a->ensemble0
     ;Dimensioned: (ensemble0, forecast_time1, g0_lat_2, g0_lon_3 )
     ; initial_time :	01/01/2021 (00:00)
      ; OR
      ; (ensemble0,initial_time1_hours, forecast_time2, g0_lat_3, g0_lon_4) 

 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)

 nens = dimsizes(gens)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Identify whether the forecast has a single or multiple years
;;; In C3S variables  dimensions are either 
;;;      ensemble0, forecast_time1, g0_lat_2, g0_lon_3
;;;      ensemble0, initial_time_hour1,forecast_time2, g0_lat_3, g0_lon_4
;;; but 
;;;     initial_time0_hours, forecast_time1,g0_lat_2, g0_lon_3
;;; for ncep with initial_time0_hours==ensemble0 constructed by 4 times daily
;;;
;;; If the forecast has forecast_time1, then the initial time is an attribute
;;; of the variable.
;;; Note this allow us to use past forecasts & consitent time axis 
;;; Also nearly the same code can be used for hindcast 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

 ii = str_match_ind(vNames,"forecast_time1")
 if (ismissing(ii)) then
     glead := a->forecast_time2
     glat  := a->g0_lat_3
     glon  := a->g0_lon_4
     glat  := glat(::-1)

     y     := gprec(initial_time1_hours|:,forecast_time2|:,\ 
                     ensemble0|:,g0_lat_3|::-1,g0_lon_4|:)

     xtime   := a->initial_time1_hours
     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) 
     nlat       = dimsizes(glat)
     nlon       = dimsizes(glon)

     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)

 else
  
     initial_time := gprec@initial_time
     initime      := grib_stime2itime(initial_time) ; YYYYMMDDHH

     iyear         = initime/1000000
     mmddhh        = initime - iyear*1000000 
     imonth        = mmddhh/10000  
     ddhh          = mmddhh - imonth*10000
     iday          = ddhh/100
     ihr           = ddhh - iday*100

     glead  := a->forecast_time1
     glat   := a->g0_lat_2
     glon   := a->g0_lon_3
     glat   := glat(::-1)
     y      := gprec(forecast_time1|:,ensemble0|:,g0_lat_2|::-1,g0_lon_3|:)

     TIME   := cd_inv_calendar(iyear,imonth,iday,ihr,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"

     nyears     = dimsizes(time)
     nlead      = dimsizes(glead)
     nlat       = dimsizes(glat)
     nlon       = dimsizes(glon)

     x := y
     x = (/ y * 1000 /) ; daily precip is in m; change to mm
     x@units  =    "mm/day"
     delete(y)

     x(1:(nlead-1),:,:,:) = x(1:(nlead-1),:,:,:) - x(0:(nlead-2),:,:,:)
         
     y  = new((/nyears,nlead,nens,nlat,nlon/),typeof(gprec),gprec@_FillValue)
     y(0:(nyears-1),:,:,:,:) = (/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)

;;;;;;;;;;;;;;;
; 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 + 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)

; end do


 TYPE       = "Hindcast" 

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;

 filemonth = allmonths(iMon)

 imodel = 0
 filepre = "./"+modelNames(imodel)+ "Daily"+ filemonth  

 fName := systemfunc(" ls -1 "+ filepre +TYPE+"*.grib")
 
 
 a   = addfile(fName,"r")


 x    := a->TP_GDS0_SFC
 gens := a->ensemble0
     ;Dimensioned: (ensemble0, forecast_time1, g0_lat_2, g0_lon_3 )
     ; initial_time :	01/01/2021 (00:00)
      ; OR
      ; (ensemble0,initial_time1_hours, forecast_time2, g0_lat_3, g0_lon_4) 

 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)

 nens = dimsizes(gens)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;; Identify whether the forecast has a single or multiple years
;;; In C3S variables  dimensions are either 
;;;      ensemble0, forecast_time1, g0_lat_2, g0_lon_3
;;;      ensemble0, initial_time_hour1,forecast_time2, g0_lat_3, g0_lon_4
;;; but 
;;;     initial_time0_hours, forecast_time1,g0_lat_2, g0_lon_3
;;; for ncep with initial_time0_hours==ensemble0 constructed by 4 times daily
;;;
;;; If the forecast has forecast_time1, then the initial time is an attribute
;;; of the variable.
;;; Note this allow us to use past forecasts & consitent time axis 
;;; Also nearly the same code can be used for hindcast 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

 ii = str_match_ind(vNames,"forecast_time1")
 if (ismissing(ii)) then
     glead := a->forecast_time2
     glat  := a->g0_lat_3
     glon  := a->g0_lon_4
     glat  := glat(::-1)

     y     := gprec(initial_time1_hours|:,forecast_time2|:,\ 
                     ensemble0|:,g0_lat_3|::-1,g0_lon_4|:)

     xtime   := a->initial_time1_hours
     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) 
     nlat       = dimsizes(glat)
     nlon       = dimsizes(glon)

     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)

 else
  
     initial_time := gprec@initial_time
     initime      := grib_stime2itime(initial_time) ; YYYYMMDDHH

     iyear         = initime/1000000
     mmddhh        = initime - iyear*1000000 
     imonth        = mmddhh/10000  
     ddhh          = mmddhh - imonth*10000
     iday          = ddhh/100
     ihr           = ddhh - iday*100

     glead  := a->forecast_time1
     glat   := a->g0_lat_2
     glon   := a->g0_lon_3
     glat   := glat(::-1)
     y      := gprec(forecast_time1|:,ensemble0|:,g0_lat_2|::-1,g0_lon_3|:)

     TIME   := cd_inv_calendar(iyear,imonth,iday,ihr,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"

     nyears     = dimsizes(time)
     nlead      = dimsizes(glead)
     nlat       = dimsizes(glat)
     nlon       = dimsizes(glon)

     x := y
     x = (/ y * 1000 /) ; daily precip is in m; change to mm
     x@units  =    "mm/day"
     delete(y)

     x(1:(nlead-1),:,:,:) = x(1:(nlead-1),:,:,:) - x(0:(nlead-2),:,:,:)
         
     y  = new((/nyears,nlead,nens,nlat,nlon/),typeof(gprec),gprec@_FillValue)
     y(0:(nyears-1),:,:,:,:) = (/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})
; nlat   := dimsizes(glat)
; nlon   := dimsizes(glon)

 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 

 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 do


end