It would be very useful to have a tool in pyfesom2 to directly compare zonal averages, over large regions such the Atlantic or Pacific basins, to observational products such as the WOA (latitude / depth plots).
For now, I first create the zonal averages of FESOM outputs with the command:
cdo --pedantic -zonmean,zonal_1 -ifthen ${fmask} -setgrid,${mesh_file} -setmissval,0 ${input_file} ${zonmean_file}, then compare them to zonal averages with a plotting function.