Possibility to read daily climate forcing for ML

[giraffffeeee]

Currently, the python code for ML is adapted for reading 6 hourly CRUJRA forcings. It could be extended to read daily or even monthly climate forcings. A flag to predefine the time step of the forcing used for ML could be helpful.

[gaillardr]

A solution could be to add a 'timestep' index in varlist.json in the "climate" dict. For instance :

"climate":

  {
    "sourcepath":"/data/rgaillard/SPINacc/test_v3/forcing_files/climate/",
    "filename":"CM61-LR-pi-03g_",
    "year_start":3800,
    "year_end":3801,
    "timestep":"3H",
    "variables":["Tair","Rainf","Snowf","Qair","PSurf","SWdown","LWdown"]
  }

Then in readvar.py, one could add if-elif statement to select the forcing timestep. A solution could be :

# 0.1 read climate variables
  climvar=varlist['climate']
  varname_clim=climvar['variables']
  days=np.array(calendar.mdays[1:])
  nyear=climvar['year_end']-climvar['year_start']+1
  tstep=climvar['timestep']
  for index in range(len(varname_clim)):
    var_month_year=np.full((nyear,12,auxil.nlat,auxil.nlon),np.nan)
    for year in range(climvar['year_start'],climvar['year_end']+1):
      check.display('reading %s from year %i'%(varname_clim[index],year),logfile)
      # DSG bugfix_start: remove hardcode
      #filename='crujra_twodeg_v2_'+str(year)+'.nc'
      #f=Dataset(climvar['sourcepath']+filename,'r')
      f=Dataset(climvar['sourcepath']+climvar['filename']+str(year)+'.nc','r')
      #DSG bugfix_end
      da=f[varname_clim[index]][:]
      #DSG: fix to read in compressed netCDF files
      if "land" in f[varname_clim[index]].dimensions:
          land=f["land"][:]-1
          ntime=len(da)
          nlat=len(f.dimensions["y"])
          nlon=len(f.dimensions["x"])
          uncomp=np.ma.masked_all((ntime,nlat*nlon))
          uncomp[:,land]=da
          da=uncomp.reshape((ntime,nlat,nlon))
      #DSG: end

      # calculate the monthly value from 6h data
      zstart=1
      var_month=np.full((12,auxil.nlat,auxil.nlon),np.nan)
      count=0
      for month in range(1,13):
        count=np.nansum(days[:month])
        if tstep=="3H":
         if month==0: 
           print("readvar : using 3H climate forcings.")
          mkk=np.mean(da[8*(zstart-1):8*count][:][:],axis=0)
        elif tstep=="6H":
          if month==0:
            print("readvar : Using 6H climate forcings.")
          mkk=np.mean(da[4*(zstart-1):4*count][:][:],axis=0)
        elif tstep=="1D":
          if month==0:
            print("readvar : Using daily climate forcings.")
          mkk=np.mean(da[(zstart-1):count][:][:],axis=0)
        elif tstep=="1M":
          if month==0:
            print("readvar : Using monthly climate forcings.")
          mkk=np.mean(da[(month-1):month][:][:],axis=0)
        else:
          raise ValueError("Climate timestep should be in [3H, 6H, 1D, 1M].")
        #mkk[da[0][:][:]==mask]=np.nan
        var_month[month-1][:][:]=mkk.filled(np.nan)
        zstart=count+1

      var_month_year[year-climvar['year_start']][:][:][:]=var_month
      #eval('MY'+varname_clim[index]+'=var_month_year')
    adict['MY%s'%varname_clim[index]]=var_month_year[:]

Here only 3H, 6H, daily and monthly timesteps are available. More can be added.

A match - case syntax would be cleaner but only works with Python >v3.10 while if - elif works with older versions.