; Get GPC GCM  daily predicted rainfall from C3S
;ZTS: Apr 2023: Started modifying Code for CFSv2. CFS daily data
;;;;            from C3S have issues to write to netCDF. 
;;;; 6Jul 2023: Restarted working on this code for CFS.    
;;;; 2Aug 2023: Restarted working on this code for CFS.    
;;;; 8Aug 2023: CFS code finalized for both hindcast and forecats
;;;;  note that ens is greated from fdtaes and HHz (ddhhall) but 
;;;; attached current year (yyyy) for descriptive purposes. For hindcast
;;;; ens is created for each year bit has the same ddhhens. Also, lead times
;;;;  are created from valid dates of nstep (lead) added to intial dates
;;;;  (initime). YYYY are again attached to initial time (xtime0/yyyymmdd)
;;;;  which vary for hindcast but are current year for forecast 
;;;; 
   
begin

 TYPE       = "Forecast" 

 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
 xMon = iMon - 1
 if (xMon .eq. 0) then
   iMon = 12
 else
   iMon = xMon
 end if
 xMon = iMon + 1

 if (iMon .eq. 1 .or. iMon .eq. 5) then 
   fdates    = (/1, 6, 11, 16, 21, 26,31 /)
 end if

 if (iMon .eq. 2 ) then 
   fdates    = (/5, 10, 15, 20, 25 /)
 end if

 if (iMon .eq. 3 .or. iMon .eq. 11 .or. iMon .eq. 12) then 
   fdates    = (/2, 7, 12, 17, 22, 27 /)
 end if

 if (iMon .eq. 4 ) then 
   fdates    = (/1, 6, 11, 16, 21, 26 /)
 end if

 if (iMon .eq. 6 .or. iMon .eq. 7) then 
   fdates    = (/5, 10, 15, 20, 25, 30 /)
 end if

 if (iMon .eq. 8 ) then 
   fdates    = (/4, 9 , 14, 19, 24, 29 /)
 end if

 if (iMon .eq. 9 .or. iMon .eq. 10) then 
   fdates    = (/3, 8 , 13, 18, 23, 28 /)
 end if

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

 filemonth = allmonths(xMon)
 
 HHz       = (/0,6,12,18/)
 ndd       = dimsizes(fdates)
 nitd      = dimsizes(fdates)
 nihr      = 4  ; always 4 times daily
 
 ddhhall  =  new(nitd*nihr,integer)
 ddhhens  =  new(nitd*nihr,integer)
 i = 0
 do idd = 0, nitd - 1
   do ihh = 0,nihr - 1
      ddhhall(i) = toint(sprinti("%0.2i",fdates(idd)) + \
                         sprinti("%0.2i",HHz(ihh)))
      ddhhens(i) = toint(sprinti("%0.4i",toint(fileyr)) + \
                         sprinti("%0.2i",fdates(idd)) + \
                         sprinti("%0.2i",HHz(ihh)))
      i = i + 1
   end do 
 end do


 nmodelStrt = 0
 nmodelLast = 0
 
 modelNames := (/"cfs" /)
 version    := (/2 /)

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

  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

;;;;;;;;

  imodel = 0 
  ;VAR_228_GDS0_SFC(initial_time0_hours,forecast_time1,g0_lat_2, g0_lon_3 )
  ;  a = addfile("./"+modelNames(imodel)+ "Daily"+ TYPE + ".grib","r") ;
  
  files = systemfunc(" ls -1 "+ "cfs*"+TYPE+".grib")
  a  = addfiles(files, "r")

  ListSetType (a, "cat") 

  print(files)

  x    := a[:]->VAR_228_GDS0_SFC
  ; CFSv2 is dimensioned: 
  ;   (initial_time0_hours, forecast_time1, g0_lat_2, g0_lon_3 )

  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)

  initial_time0_hours  = a[:]->initial_time0_hours
     

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

  forecast_time1 := a[0]->forecast_time1  ;ftime1==nlead doesnt change
  initial_time := a[:]->initial_time0 ; initime does change
  initime      := grib_stime2itime(initial_time) ; YYYYMMDDHH

  iyear         = initime/1000000
  yrs           = get_unique_values(iyear)

  glead  := forecast_time1 ; does not change for any of the grib files
  glat   := a[0]->g0_lat_2 ; does not change for any of the grib files
  glon   := a[0]->g0_lon_3 ; does not change for any of the grib files
  glat   := glat(::-1)

  nlead  := dimsizes(forecast_time1)
  nlat   := dimsizes(glat)
  nlon   := dimsizes(glon)
  nyears := dimsizes(yrs)
         
         
  y :=  gprec(initial_time0_hours|:,forecast_time1|:,g0_lat_2|::-1,g0_lon_3|:)

  x := y
  x = (/ y * 1000 /) ; daily precip is in m; change to mm
  x@units  =    "mm/day"
   
  nens   := nitd*nihr ; ensembles are created by intital date and time
  ntime0 := nlead + nens + 4  ; this is the maximum time length 
                                    ; need to find out why i=20 yyyymmdd is
                                    ; not long enough 
  yyyymmdd1 := new(nyears*ntime0,integer)
           
  do i = 1, nyears 
     yr0                = yrs(i-1)
     yyyymmdd          := yyyymmdd_time(yr0, yr0+1, "integer")
     yri               := toint(yr0+sprinti("%0.2i",iMon)+"01")
     yyyymmdd          := yyyymmdd({yri:})
     ixs                = (i-1)*ntime0     
     ixe                = ixs + ntime0 - 1
     ;print(ixs+ " " + ixe +" " + ntime0 )
     ;zz2 = dimsizes(ispan(0,(ntime0-1),1))
     ;zz1 = dimsizes(ispan(ixs,ixe,1))
     ;zz3 = dimsizes(yyyymmdd1)
     ;zz4= dimsizes(yyyymmdd)
     ;print(zz1 + " " +zz2 + " " + zz3 + " " + zz4)
     yyyymmdd1(ixs:(ixe)) = yyyymmdd(0:(ntime0-1))  
     ;print(yyyymmdd(0:(ntime0-1)) + "  " + yyyymmdd1(ixs:ixe))
  end do

  yr0                = yrs(0)
  yyyymmdd          := yyyymmdd_time(yr0, yr0+1, "integer")
  yri               := toint(yr0+sprinti("%0.2i",iMon)+"01")
  yyyymmdd          := yyyymmdd({yri:})
  xtime0            := yyyymmdd(0:(ntime0-1))
  xtime0!0           = "initime"
  xtime0@units       = "yyyymmdd"
  mmddhh            := xtime0  ; this is the actual time a lead will be valid 

 ; print(xtime0)
  
  yyyymmdd := yyyymmdd1
  ;print(mmddhh)
  xshp := new((/ntime0*nyears,nens,nlat,nlon/),typeof(gprec),gprec@_FillValue) 

  xshp!0      = "time0"
  xshp!1      = "ens"
  xshp!2      = "lat"
  xshp!3      = "lon"
  xshp&time0  = yyyymmdd
  xshp&ens    = ispan(0,nens - 1,1)
  xshp&lat    = glat
  xshp&lon    = glon
       
  xi := initial_time0_hours
      
  do i=0,dimsizes(xi) - 1
             
     xt    := xi(i)
     xto   := cd_calendar(xt,0)
     xhr   := toint(xto(0,3)/24. + 0.49) ; 06:1200hrs will be the same day
                                    ;but 1300:2300hrs considered next day
     xyr   := xto(0,0)
     xyri  = get1Dindex(yrs,xyr)

     ddhh := toint(sprinti("%0.2i",toint(xto(0,2))) + \
                   sprinti("%0.2i",toint(xto(0,3)))) ; use to locate ensi
 
     xensi  = get1Dindex(ddhhall,ddhh)
   
     xt1    := xt + forecast_time1
     ftind  := ispan(0,dimsizes(forecast_time1) - 1, 1)
     copy_VarAtts(xt,xt1)
     
     xt2    := cd_calendar(xt1,-3) ; yyyymmddHH
           
     xt3    := xt2/100    ;here time of the day is discarded
           
          
     if (xhr .eq. 1) then  
         xt3   := xt3(1:(dimsizes(xt3)-1)) ; nextday for 12/18Z
         ftind := ispan(1,dimsizes(forecast_time1) - 1, 1)
     end if
     zind := get1Dindex(yyyymmdd,xt3)
 
     if(any(ismissing(zind))) then
           print("zind contains missing values for i="+i)
           igood := ind(.not.ismissing(zind))
           zind  := zind(igood)
           ftind := ftind(igood)
     end if

     xshp(zind,xensi,:,:) = (/ x(i,ftind,:,:) /)

  end do ; for each initial time

  x       := reshape(xshp, (/nyears,ntime0,nens,nlat,nlon/))
  xlead   := ispan(0,ntime0-1,1) 
  xens    := ispan(0,nens - 1,1)
  nlead   := ntime0
  x!0     = "time"
  x!1     = "lead"
  x!2     = "ens"
  x!3     = "lat"
  x!4     = "lon"
  x&time  =  yrs
  x&lead  = xlead
  x&ens   = xens
  x&lat   = glat
  x&lon   = glon

  nlds = nlead - 1 
  x(:,1:nlds,:,:,:) = (/ x(:,1:nlds,:,:,:) - x(:,0:(nlds-1),:,:,:)  /)
       
  iyears  = yrs
  imonth  = yrs
  iday    = yrs
  ihr     = yrs

  imonth  = (/iMon /)
  iday    = (/ 1 /)
  ihr     = (/ 0 /)
  imn     = ihr
  isc     = ihr  
 
  TIME      := cd_inv_calendar(iyears,imonth,iday,ihr,imn,isc,\
                            "days since 1990-01-01 00:00:0",0)
  TIME!0     = "time"

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

  ;print(time)
  y          := x
  delete(x)

  glead := xlead ; xtime0
  gens  := xens 

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

  y&time = time

  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    := ddhhens   ;ispan(1,nens, 1) 
  lead   := xlead  ;ispan(0,nlead-1,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 dates"

  lead!0          = "lead"
  lead@units      = "oftset from initial date "
  lead@long_name  = "Forecast offset from initial time"

  xlead          := xtime0
  xlead!0         = "lead"
  xlead@units     = "yyyymmdd"
  xlead@long_name = "initial time for each year"

  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"

  ; print(initime)
  prec&time       = time
  prec&lead       = lead
  prec&ens        = ddhhens ;ens created from initial dates & hrs of year
  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 " + 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 /)
  ;print(dimSizes)
  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(time), getvardims(time) )
  filevarattdef(fout,"initime",time)

  filevardef(fout,"timelead",typeof(xlead), getvardims(xlead) )
  filevarattdef(fout,"timelead",xlead)

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


;***********************************************
; terminate define mode: not necessary but for clarity
;*********************************************** 
  setfileoption(fout, "definemode", False)
;  print(fout)
;===================================================================
; 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   = (/ time /)
  fout->timelead  = (/ xlead /)
  fout->prec      = (/ prec /)

  delete(prec)
  delete(time)
  delete(lead)
  delete(ens)
  delete(glat)
  delete(glon)
  delete(xlead)
  delete(xtime0)
  delete(fout)
;;;;;;;;
;;;;  Do for hindcast below
;;;;;;;;
  

 TYPE       = "Hindcast"

 glatS = latS - 1
 glatN = latN + 1
 glonW = lonW - 1
 glonE = lonE + 1
 
 
  imodel = 0 
  
  files := systemfunc(" ls -1 "+ "cfs*"+TYPE+"*.grib")
  a  := addfiles(files, "r")

  ListSetType (a, "cat") 

  print(files)

  x    := a[:]->VAR_228_GDS0_SFC
  ; CFSv2 is dimensioned: 
  ;   (initial_time0_hours, forecast_time1, g0_lat_2, g0_lon_3 )

  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)

  initial_time0_hours  := a[:]->initial_time0_hours
     

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

  forecast_time1 := a[0]->forecast_time1  ;ftime1==nlead doesnt change
  initial_time := a[:]->initial_time0 ; initime does change
  initime      := grib_stime2itime(initial_time) ; YYYYMMDDHH

  iyear         := initime/1000000
  yrs           := get_unique_values(iyear)

  glead  := forecast_time1 ; does not change for any of the grib files
  glat   := a[0]->g0_lat_2 ; does not change for any of the grib files
  glon   := a[0]->g0_lon_3 ; does not change for any of the grib files
  glat   := glat(::-1)

  nlead  := dimsizes(forecast_time1)
  nlat   := dimsizes(glat)
  nlon   := dimsizes(glon)
  nyears := dimsizes(yrs)
         
         
  y :=  gprec(initial_time0_hours|:,forecast_time1|:,g0_lat_2|::-1,g0_lon_3|:)

;  y := y(:,:,{glatS:glatN}, {glonW:glonE})
;  glat  := glat({glatS:glatN})
;  glon  := glon({glonW:glonE})
;  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) 
  nens   := nitd*nihr ; ensembles are created by intital date and time
  ntime0 := nlead + nens + 4  ; this is the maximum time length 
                                    ; need to find out why i=20 yyyymmdd is
                                    ; not long enough 
  yyyymmdd1 := new(nyears*ntime0,integer)
           
  do i = 1, nyears 
     yr0                = yrs(i-1)
     yyyymmdd          := yyyymmdd_time(yr0, yr0+1, "integer")
     yri               := toint(yr0+sprinti("%0.2i",iMon)+"01")
     yyyymmdd          := yyyymmdd({yri:})
     ixs                = (i-1)*ntime0     
     ixe                = ixs + ntime0 - 1
     yyyymmdd1(ixs:(ixe)) = yyyymmdd(0:(ntime0-1))  
  end do

  yr0                = yrs(0)
  yyyymmdd          := yyyymmdd_time(yr0, yr0+1, "integer")
  yri               := toint(yr0+sprinti("%0.2i",iMon)+"01")
  yyyymmdd          := yyyymmdd({yri:})
  xtime0            := yyyymmdd(0:(ntime0-1))
  xtime0!0           = "initime"
  xtime0@units       = "yyyymmdd"
  mmddhh            := xtime0  ; this is the actual time a lead will be valid 

 ; print(xtime0)
  
  yyyymmdd := yyyymmdd1
  ;print(mmddhh)
  xshp := new((/ntime0*nyears,nens,nlat,nlon/),typeof(x),x@_FillValue) 

  xshp!0      = "time0"
  xshp!1      = "ens"
  xshp!2      = "lat"
  xshp!3      = "lon"
  xshp&time0  = yyyymmdd
  xshp&ens    = ispan(0,nens - 1,1)
  xshp&lat    = glat
  xshp&lon    = glon
       

  xi := initial_time0_hours
      
  do i=0,dimsizes(xi) - 1
             
     xt    := xi(i)
     xto   := cd_calendar(xt,0)
     xhr   := toint(xto(0,3)/24. + 0.49) ; 06:1200hrs will be the same day
                                    ;but 1300:2300hrs considered next day
     xyr   := xto(0,0)
     xyri  = get1Dindex(yrs,xyr)

     ddhh := toint(sprinti("%0.2i",toint(xto(0,2))) + \
                   sprinti("%0.2i",toint(xto(0,3)))) ; use to locate ensi
 
     xensi  = get1Dindex(ddhhall,ddhh)
   
     xt1    := xt + forecast_time1
     ftind  := ispan(0,dimsizes(forecast_time1) - 1, 1)
     copy_VarAtts(xt,xt1)
     
     xt2    := cd_calendar(xt1,-3) ; yyyymmddHH
           
     xt3    := xt2/100    ;here time of the day is discarded
           
          
     if (xhr .eq. 1) then  
         xt3   := xt3(1:(dimsizes(xt3)-1)) ; nextday for 12/18Z
         ftind := ispan(1,dimsizes(forecast_time1) - 1, 1)
     end if
     zind := get1Dindex(yyyymmdd,xt3)
 
     if(any(ismissing(zind))) then
           print("zind contains missing values for i="+i)
           igood := ind(.not.ismissing(zind))
           zind  := zind(igood)
           ftind := ftind(igood)
     end if

     xshp(zind,xensi,:,:) = (/ x(i,ftind,:,:) /)

  end do ; for each initial time

  x       := reshape(xshp, (/nyears,ntime0,nens,nlat,nlon/))
  xlead   := ispan(0,ntime0-1,1) 
  xens    := ispan(0,nens - 1,1)
  nlead   := ntime0
  x!0     = "time"
  x!1     = "lead"
  x!2     = "ens"
  x!3     = "lat"
  x!4     = "lon"
  x&time  =  yrs
  x&lead  = xlead
  x&ens   = xens
  x&lat   = glat
  x&lon   = glon

  nlds = nlead - 1 
  x(:,1:nlds,:,:,:) = (/ x(:,1:nlds,:,:,:) - x(:,0:(nlds-1),:,:,:)  /)
       
  iyears  := yrs
  imonth  := yrs
  iday    := yrs
  ihr     := yrs

  imonth  = (/iMon /)
  iday    = (/ 1 /)
  ihr     = (/ 0 /)
  imn     := ihr
  isc     := ihr  
 
  TIME      := cd_inv_calendar(iyears,imonth,iday,ihr,imn,isc,\
                            "days since 1990-01-01 00:00:0",0)
  TIME!0     = "time"

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

  ;print(time)
  y          := x
  delete(x)

  glead := xlead ;xtime0
  gens  := xens 

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

  y&time = time

  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    := ddhhens   ;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 dates"

  lead!0          = "lead"
  lead@units      = "oftset from initial date "
  lead@long_name  = "Forecast offset from initial time"

  xlead          := xtime0
  xlead!0         = "lead"
  xlead@units     = "yyyymmdd"
  xlead@long_name = "initial time for each year"

  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"

  ; print(initime)
  prec&time       = time
  prec&lead       = lead
  prec&ens        = ddhhens ;ens created from initial dates & hrs of year
  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 /)
  ;print(dimSizes)
  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(time), getvardims(time) )
  filevarattdef(fout,"initime",time)

  filevardef(fout,"timelead",typeof(xlead), getvardims(xlead) )
  filevarattdef(fout,"timelead",xlead)

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

;***********************************************
; terminate define mode: not necessary but for clarity
;*********************************************** 
  setfileoption(fout, "definemode", False)
;  print(fout)
;===================================================================
; 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  = (/ time /)
  fout->timelead = (/ xlead /)
  fout->prec     = (/ prec /)

  delete(prec)
  delete(time)
  delete(lead)
  delete(ens)
  delete(glat)
  delete(glon)
  delete(xlead)
  delete(fout)
;;;;;;;
;;; End of Hindcast
;;;;;;;;
end