Errors around coast lines on SST plots with v1-HR model output

[golaz]

I ran e3sm_diags on output from v1-HR, regridded using map_ne120np4_to_cmip6_720x1440_nco.20190601.nc.

The SST plots show erroneous values near the coastlines, leading to a very large RMSE error:

Full E3SM Diags output: https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/ac.golaz/E3SM/v1-HR/20210112.A_WCYCL1950S_CMIP6_HR.ne120_oRRS18v3_ICG.unc06/e3sm_diags/720x1440_nco/model_vs_obs_0137-0155/viewer/

Here is my configuration to reproduce on chrysalis using zppy:

[default]
input = /lcrc/group/e3sm/ac.ndkeen/scratch/chrys/E3SM_simulations/production-unc06
input_subdir = run
output = /lcrc/group/e3sm/ac.golaz/E3SM_analysis/20210112.A_WCYCL1950S_CMIP6_HR.ne120_oRRS18v3_ICG.unc06
case = 20210112.A_WCYCL1950S_CMIP6_HR.ne120_oRRS18v3_ICG.unc06
www = /lcrc/group/e3sm/public_html/diagnostic_output/ac.golaz/E3SM/v1-HR
partition = compute
e3sm_unified = latest

[climo]
active = True
years = "137-155",

  [[ atm_monthly_720x1440_nco ]]
  input_files = "cam.h0"
  mapping_file = /home/ac.zender/data/maps/map_ne120np4_to_cmip6_720x1440_nco.20190601.nc

[ts]
active = True
years = "137:155:1",

  [[ atm_monthly_720x1440_nco ]]
  input_files = "cam.h0"
  mapping_file = /home/ac.zender/data/maps/map_ne120np4_to_cmip6_720x1440_nco.20190601.nc
  vars = "FSNTOA,FLUT,FSNT,FLNT,FSNS,FLNS,SHFLX,QFLX,TAUX,TAUY,PRECC,PRECL,PRECSC,PRECSL,TS,TREFHT,CLDTOT,CLDHGH,CLDMED,CLDLOW"

  [[ atm_monthly_glb ]]
  input_files = "cam.h0"
  mapping_file = "glb"
  vars = "FSNTOA,FLUT,FSNT,FLNT,FSNS,FLNS,SHFLX,QFLX,TAUX,TAUY,PRECC,PRECL,PRECSC,PRECSL,TS,TREFHT,CLDTOT,CLDHGH,CLDMED,CLDLOW"

[glb]
active = False

[e3sm_diags]
active = True
years = "137-155",

  [[ atm_monthly_720x1440_nco ]]
  short_name = 'unc06'
  grid = '720x1440_nco'
  reference_data_path = '/lcrc/soft/climate/e3sm_diags_data/obs_for_e3sm_diags/climatology'
  output_format_subplot = "pdf",

[e3sm_diags_vs_model]
active = False

[amwg]
active = False

[mpas_analysis]
active = True
walltime = "24:00:00"
parallelTaskCount = 6
subdir_ocean = run
subdir_ice = run
ts_years = "137-155",
enso_years = "137-155",
climo_years ="137-155",
mesh = "oRRS18to6v3"

[chengzhuzhang]

I can reproduce the problem with:

e3sm_diags lat_lon --no_viewer --case_id 'SST_CL_HadISST' --sets 'lat_lon' --run_type 'model_vs_obs' --variables 'SST' --seasons 'ANN' --regions 'global' --regrid_tool 'esmf' --regrid_method 'conservative' --main_title 'SST ANN global' --backend 'mpl' --output_format 'png' --output_format_subplot 'pdf' --canvas_size_w '1212' --canvas_size_h '1628' --figsize '8.5' '11.0' --dpi '150' --arrows --contour_levels '-1' '0' '1' '3' '6' '9' '12' '15' '18' '20' '22' '24' '26' '28' '29' --test_name '20210112.A_WCYCL1950S_CMIP6_HR.ne120_oRRS18v3_ICG.unc06' --short_test_name 'unc06' --test_colormap 'cet_rainbow.rgb' --ref_name 'HadISST_CL' --reference_name 'HadISST/OI.v2 (Climatology)' --reference_colormap 'cet_rainbow.rgb' --diff_title 'Model - Observations' --diff_colormap 'diverging_bwr.rgb' --diff_levels '-5' '-4' '-3' '-2' '-1' '-0.5' '-0.2' '0.2' '0.5' '1' '2' '3' '4' '5'  --granulate 'variables' 'seasons' 'plevs' 'regions' --selectors 'sets' 'seasons' --test_data_path '/lcrc/group/e3sm/ac.golaz/E3SM_analysis/20210112.A_WCYCL1950S_CMIP6_HR.ne120_oRRS18v3_ICG.unc06/post/atm/720x1440_nco/clim/19yr' --reference_data_path '/lcrc/soft/climate/e3sm_diags_data/obs_for_e3sm_diags/climatology' --results_dir '/lcrc/group/e3sm/public_html/diagnostic_output/zhang40' --save_netcdf

Could be a regridding issue when taking the difference between model and obs. I'm looking into it.

[chengzhuzhang]

I tracked back to an old issue which is not being addressed: #175
This issue affects evaluating high resolution model data with a mask being applied, ig. SST, T at 850mb. The problem is with conservative regriding method in cdms. At the time being, CDMS issues are not being worked on, I would propose to use bilinear interpolation as a default regriding method.
Conservative for high res SST:

Bilinear for high res SST:

For SST plots, switching to bilinear interpolation doesn't seem to change the metrics for Low res:
Conservative for low res SST:

Bilinear for low res SST:

I'm going to run through all variables to see what fields will be affected.

[chengzhuzhang]

I ran two cases: one uses bilinear interpolation the other uses conservative.

Checking the difference between ANN table: https://web.lcrc.anl.gov/public/e3sm/diagnostic_output/zhang40/tests/full_bilinear_low-res/viewer/table/index.html
All radiation fluxes fields stay identical, many other fields see minor changes. TGCLDLWP_OCN has largest difference comparing two methods.

[golaz]

@chengzhuzhang : thanks for the nice detective work here! I think we can assume that the CDAT issue #175 is not going to get fixed anytime soon, likely ever.

Switching (temporarily) to bi-linear is probably the best option we have, even though it makes me a little uncomfortable.

This brings into question again the future of CDAT in E3SM Diags, and whether we should transition to something else.

[chengzhuzhang]

Thank you, @golaz . I am exploring using xarray and xesmf for regridding instead. Again, the bi-linear method works okay. But conservative method fails due to errors with getting data bounds. I will try trouble shooting further. If it won't work, we may have to switch to bi-linear as a short-term solution.

[tangq]

@chengzhuzhang , we use nco for remapping in zppy. You may want to use nco or TempestRemap to fix this issue. The advantages are that these are tools supported by E3SM and easy to use. Just a thought.

[chengzhuzhang]

Hi @tangq, thanks for the suggestion. A solution that uses a Python interface would be more straightforward to implement in this case.

@golaz, I experimented more with xarray and xesmf regridder and found that conservative regridding looks promising. To do so, a cdms variable needs to be converted to an xarray variable, with masks and bounds reformed. Example 3 for xarray here: https://github.com/lee1043/XCDAT/blob/master/compare_cdat_xarray/Regrid_ex1.ipynb

To implement, we will need to add additional dependencies, also a full test run is needed to find out the difference between the new regridding and current cdms2 regridding method.

[tangq]

@chengzhuzhang , I made the suggestion of using nco or TempestRemap mainly echoing @golaz's comment "This brings into question again the future of CDAT in E3SM Diags". I was experiencing a similar dependency issue when I saw this discussion. Another frustrating case is the cdscan dependency when creating nudging files (You were in that conversation too).

If we have to add dependencies, why not depend on something E3SM supports when possible?

It should be straightforward to implement it. You call the following ncremap command with python before calculating the model-obs differences.
ncremap -a conserve -i input.nc -d dst.nc -m map.nc -o out.nc (if the input.nc and dst.nc are on lat-lon grids. Charlie can provide more instructions on the ncremap options for various input.nc and dst.nc file formats)

Then the python code uses regridded out.nc for the differencing.

Just a suggestion. Feel free to proceed either way. But I'd advocate E3SM supported tools.

[chengzhuzhang]

Yeah, we should use E3SM supported tools as much as it is practicable (I'm HUGE fan of NCO and Tempest). The regridding operation is one of the example to operate a variable object extracted from netcdf files, in e3sm-diags, we also do regional subsetting, vertical interpolation, statistical computation... In this particular use case, using cdms/xarray or their equivalent will be more suitable...

[czender]

@chengzhuzhang and @golaz

Large model/obs SST differences along coastlines (or land model/obs differences along coastlines) are an almost fool-proof indication that the regridding weights are being applied with a global mean-preserving, not a local mean-preserving, algorithm, as described here.

ncremap will, if invoked in MPAS mode with '-P mpas' as described here automatically use the appropriate algorithm (local mean-preserving renormalization) so that the coastline problem "disappears".

This is usually only an issue when using weights generated with a flux-preserving conservative method rather than bilinear interpolation. Let me know if I can be of any help.