diff --git a/openquake/cat/completeness/analysis.py b/openquake/cat/completeness/analysis.py
index 77b61e0ff..d0dc80154 100644
--- a/openquake/cat/completeness/analysis.py
+++ b/openquake/cat/completeness/analysis.py
@@ -310,7 +310,7 @@ def _completeness_analysis(fname, years, mags, binw, ref_mag, ref_upp_mag,
         print(f'Iteration: {iper:05d} norm: {norm:12.6e}', end="\r")
 
         ctab = _make_ctab(prm, years, mags)
-        print(ctab)
+        #print(ctab)
         if isinstance(ctab, str):
             continue
 
diff --git a/openquake/ghm/mosaic.py b/openquake/ghm/mosaic.py
index 1560fe5e7..3a3f16557 100644
--- a/openquake/ghm/mosaic.py
+++ b/openquake/ghm/mosaic.py
@@ -39,6 +39,22 @@
            'nea': {'RUS': ['75 0 180 0 180 88 75 88',
                            '-180 0 -160 0 -160 88 -180 88']},
            'nwa': {'RUS': ['75 0 25 0 25 10 25 20 25 40 25 50 25 88 75 88']},
+           'pan': {'ATF': ['49.5 -50.1 71 -50.1 71 -45.6 49.5 -45.6']},
+           'oat': {'BRA': ['-33.34 -21.5 -27.18 -21.5 -27.18 -2.44 -33.34 -2.44'],
+                   'ATA': ['-63.76 -64.82 -43.7 -64.82 -43.7 -59.91 -63.76 -59.91'],
+                   'ATF': ['46.48 -16.71 55.2 -16.71 55.2 -11.35 46.48 -11.35']},
+           'oin': {'IND': ['71.32 7.686 74.224 7.686 74.224 12.10 71.32 12.10'],
+                   'MUS': ['63.27 -19.83 63.53 -19.83 63.53 -19.62 63.27 -19.62'],
+                   'ATF': ['77.4 -38.76 77.64 -38.76 77.64 -37.73 77.4 -37.73']},
+           'opa': {'AUS': ['158 -55.5 160 -55.5 160 -31 158 -31'],
+                   'JPN': ['128.58 23.5 155 23.5 155 28.55 128.58 28.55'],
+                   'CHL': ['-110.1 -33.9 -77.7 -33.9 -77.7 -24.2 -110.1 -24.2'],
+                   'ECU': ['-94.14 -2.50 -88.37 -2.50 -88.37 2.54 -94.14 2.54'],
+                   'CRI': ['-87.78 5.07 -86.16 5.07 -86.16 6.10 -87.78 5.07'],
+                   'MEX': ['-114.99 17.10 -109.37 19.59 -117.75 29.84 -121.38 26.69'],
+                   'NZL': ['169.38 -53.29 179 -50.38 178.94 -46.77 166.61 -47.61 165.37 -50.81',
+                           '171.8 -34.44 172.3 -34.44 172.3 -33.85 171.8 -33.85',
+                           '-177.63 -44.60 -175.44 -44.60 -175.44 -43.12 -177.63 -43.12']},
            'pac': {'NZL': ['-180.00 -32.00 -180.00 -3.47 -174.35 -5.49 -169.62 -8.96 -166.25 -13.69 -166.83 -22.66 -172.90 -32.00 -180.00 -32.00',
                            '156.22 -6.28 155.95 -6.81 155.88 -6.81 155.83 -6.84 155.49 -6.98 154.58 -7.29 159.72 -11.44 165.23 -23.85 167.73 -29.60 178.01 -32.55 180.00 -33.01 180.00 -3.47 177.66 -3.26 169.67 0.20 160.26 -0.64 156.22 -6.28'],
                    'FJI': ['-180.00 -32.00 -180.00 -3.47 -174.35 -5.49 -169.62 -8.96 -166.25 -13.69 -166.83 -22.66 -172.90 -32.00 -180.00 -32.00',
@@ -178,11 +194,18 @@
                'MAR', 'MLI', 'MRT', 'NER', 'PSA',
                'SDN', 'TCD', 'TUN'],
        'nzl': ['NZL'],
-        #'oat': ['']
-        #'oin': ['']
+       'oan': ['ATF', 'HMD'],
+       'oat': ['CPV', 'STP', 'SHN', 'SYC','MUS',
+               'REU', 'GNQ', 'MYT', 'COM','BMU',
+               'FLK', 'SGS', 'BVT', 'ATA', 'BRA'],
+       'oin': ['MDV', 'CCK', 'CXR', 'MUS', 
+               'IOT', 'IND', 'ATF'], 
         # KIR - Kiribati
         # COK - Cook Islands
-       'opa': ['COK', 'KIR'],
+       'opa': ['KIR', 'PLW', 'FSM', 'PCN', 'COK', 
+               'PYF', 'XCL','UMI', 'CRI', 'MNP', 
+               'GUM', 'MHL', 'NRU', 'CHL', 'ECU', 
+               'JPN', 'MEX', 'NZL', 'AUS'],
         # ASM - American Samoa
         # FJI - Fiji
         # KIR - Kiribati
diff --git a/openquake/man/source_tests.py b/openquake/man/source_tests.py
new file mode 100644
index 000000000..201537880
--- /dev/null
+++ b/openquake/man/source_tests.py
@@ -0,0 +1,264 @@
+import os
+import glob
+import numpy as np
+import pandas as pd
+import geopandas as gpd
+import pygmt
+from pathlib import Path
+import pathlib
+from scipy.stats import poisson
+from tabulate import tabulate
+import matplotlib.pyplot as plt
+from openquake.hazardlib.nrml import to_python
+from openquake.hazardlib.sourceconverter import SourceConverter
+
+
+DEFAULT = SourceConverter(
+	investigation_time=1.0,
+	rupture_mesh_spacing=10.,
+	area_source_discretization=10.)
+
+def get_params(ssm, mmin = 0, total_for_zone = False):
+	'''
+	Retrieve source parameters from xml source file for multipoint source
+	Specify total_for_zone to retain only total a value summed over all points
+
+	:param ssm:
+	    seismic source model retreived from xml file with to_python
+	:param mmin:
+	    minimum magnitude to consider
+	:param total_for_zone:
+	    if True, returns the total a-value for the zone. 
+	    if False, returns a-values at point locations
+	'''
+	total_rate_above_zero = 0
+	data = []
+	for ps in ssm[0][0]:
+		agr = ps.mfd.a_val
+		bgr = ps.mfd.b_val
+		mmax = ps.mfd.max_mag
+		lo = ps.location.longitude
+		la = ps.location.latitude
+		log_rate_above = agr - bgr * mmin
+		total_rate_above_zero += 10**(agr)
+			
+		data.append([lo, la, log_rate_above, bgr, mmax])
+
+	if total_for_zone == True:  
+		#total_rate_above += 10**(agr - bgr * m_min)
+		return np.log10(total_rate_above_zero), bgr, mmax
+	else: 
+		return np.array(data)
+
+def plot_sources(folder_name, region, mmin, fig_folder,
+				 rate_range = [-9, 3, 0.2], mmax_range = [6.0, 9.0, 0.2], 
+				 sconv = DEFAULT, poly_name = '', plot_poly = False):
+	'''
+	Create plots of source rates and mmax for each source in folder_name
+
+	:param folder_name:
+		folder containing source xml files
+	:param region:
+	    region for pygmt plotting
+	:param mmin:
+		minimum magnitude to be used in model
+	:param fig_folder:
+		folder in which to store plots
+	:param rate_range:
+		range of rate to use in colour scale,
+		specified as [min, max, interval]
+	:param mmax_range:
+	    range of mmax for colour scale, 
+	    specified as [min, max, interval]
+	:param poly_name:
+		location of file containing the model source polygon (optional)
+	:param plot_poly:
+	    boolean specifying if polygon outline should be plotted
+	'''
+	path = pathlib.Path(folder_name)
+
+	# make rate and mmax folders if they do not exist
+	if os.path.exists(os.path.join(os.path.join(fig_folder, 'rate'))):
+		print("found folders at ", os.path.join(os.path.join(fig_folder, 'rate')))
+	else:
+		os.mkdir(os.path.join(fig_folder, 'rate'))
+		os.mkdir(os.path.join(fig_folder, 'mmax'))
+
+	# set up colour scales
+	cpt_rate = os.path.join(fig_folder, 'rate.cpt')
+	pygmt.makecpt(cmap="turbo", series=rate_range, output=cpt_rate)
+	cpt_mmax = os.path.join(fig_folder, 'mmax.cpt')
+	pygmt.makecpt(cmap="rainbow", series=mmax_range, output=cpt_mmax)
+    
+	if poly_name:
+    	# set plotting polygon
+		poly = gpd.read_file(poly_name)
+		poly["x"] = poly.representative_point().x
+		poly["y"] = poly.representative_point().y
+
+	for fname in sorted(path.glob('src*.xml')):
+		ssm = to_python(fname, sconv)
+		fig_a = pygmt.Figure()
+		fig_a.basemap(region=region, projection="M15c", frame=True)
+		fig_a.coast(land="grey", water="white")
+
+		fig_b = pygmt.Figure()
+		fig_b.basemap(region=region, projection="M15c", frame=True)
+		fig_b.coast(land="grey", water="white")    
+    
+		vmin = +1e10
+		vmax = -1e10
+    
+		for grp in ssm:
+			for src in grp:
+				name = src.name    
+        
+		data = get_params(ssm, mmin=mmin, total_for_zone = False)
+		vmin = np.min([vmin, min(data[:,2])])
+		vmax = np.max([vmin, max(data[:,2])])
+
+		fig_a.plot(x=data[:,0], 
+			y=data[:,1], 
+			style="h0.2", 
+			color=data[:, 2],
+			cmap=cpt_rate)
+ 
+		fig_b.plot(x=data[:,0], 
+			y=data[:,1], 
+			style="h0.2", 
+			color=data[:, 4],
+			cmap=cpt_mmax)
+        
+		if plot_poly == True:
+			fig_a.plot(data=poly, pen=".5p,black")
+			fig_b.plot(data=poly, pen=".5p,black")
+    
+		fig_a.colorbar(frame=f'af+l"Log((N(m)>{mmin}))"', cmap=cpt_rate)    
+		fig_b.colorbar(frame='af+l"Mmax"', cmap=cpt_mmax)
+    
+		out = os.path.join(fig_folder, 'rate', name+'rate.png')
+		fig_a.savefig(out)
+    
+		out = os.path.join(fig_folder, 'mmax', name+'_mmax.png')
+		fig_b.savefig(out)
+
+def simulate_occurrence(agr, bgr, rate, minmag, mmax, time_span, N=2000):
+	'''
+	Simulate number of occurrences from a Poisson distribution given the FMD parameters
+
+	:param agr:
+		a value for source
+	:param bgr: 
+	    b value for source
+	:param minmag:
+	    minimum magnitude to be considered 
+	:param mmax: 
+        maximum magnitude
+    :param time_span:
+        time span (in years) over which to simulate occurrences
+    :param N:
+        Number of simulations. Default to 2000
+	'''
+	num_occ = np.random.poisson(rate*time_span, N)
+	return(num_occ)
+
+def occurence_table(path_oq_input, path_to_subcatalogues, minmag, minyear, maxyear, N, src_ids, sconv = DEFAULT):
+	'''
+	Check number of events expected from the source model against the number of observations.
+	Uses N samples from a Poisson distribution with rate from source a and b value.
+	Returns a table summarising the catalogue vs source model for zones in src_ids.
+
+	:param path_oq_input: 
+	    path to location of xml source models
+	:param path_to_subcatalogues:
+	    path to subcatalogues to compare source with
+	:param minmag:
+	    minimum magnitude to consider
+	:param minyear: 
+	    year to start analysis
+	:param maxyear: 
+	    end year for analysis
+	:param N:
+	    number of Poisson samples to use for comparison
+	:param src_ids:
+	    list of sources to use
+	:param sconv:
+	    source converter object specifying model setup
+	'''
+	table = []
+	time_span = maxyear - minyear
+
+	for src_id in sorted(src_ids):
+
+		fname_src = os.path.join(path_oq_input, "src_{:s}.xml".format(src_id))
+		ssm = to_python(fname_src, sconv)
+
+		fname_cat = os.path.join(path_to_subcatalogues, "subcatalogue_zone_{:s}.csv".format(src_id))
+		df = pd.read_csv(fname_cat)
+		df = df.loc[df.magnitude > minmag]
+		df = df.loc[(df.year >= minyear) & (df.year <= maxyear)]
+		obs = len(df)
+
+		agr, bgr, mmax = get_params(ssm, minmag, total_for_zone = True) 
+		rate = 10.0**(agr-minmag*bgr)-10.0**(agr-mmax*bgr)
+		num_occ_per_time_span = simulate_occurrence(agr, bgr, rate, minmag, mmax, time_span, N)
+    	
+		mioc = min(num_occ_per_time_span)
+		maoc = max(num_occ_per_time_span)
+    
+		perc_16 = np.percentile(num_occ_per_time_span, 16)
+		perc_84 = np.percentile(num_occ_per_time_span, 84)
+    
+		perc_16 = poisson.ppf(0.16, rate*time_span)
+		perc_84 = poisson.ppf(0.84, rate*time_span)
+    
+		agr_cat = np.nan
+		if obs > 1e-10:
+			agr_cat = np.log10(obs/time_span) + bgr * minmag
+    
+		pss = ""
+		if obs >= perc_16 and obs <= perc_84:
+			pss = "="
+		elif obs < perc_16:
+			pss = "<"
+		else:
+			pss = ">"
+
+		table.append([src_id, agr_cat, agr, bgr, mioc, maoc, perc_16, perc_84, obs, pss])
+	
+	heads = ["Zone", "agr_cat", "agr", "bgr", "min", "max", "%16", "%84", "observed", "obs Vs. pred"]
+	print(tabulate(table, headers=heads))
+
+def source_info_table(folder_name, sconv = DEFAULT):
+	'''
+	Print a table describing the sources in this model, inclduing their ID, upper and lower depth limits,
+	tectonic region, magnitude scaling relation, magnitude limits and depth distributions.
+
+	:param folder_name:
+	    folder to find soure xmls
+	:param sconv:
+	    source converter object
+	'''
+	columns = ['ID', 'Name', 'TR', 'USD', 'LSD', 'MSR', 'Mmin', 'Mmax', 'h_depths', 'b_vals']
+	sdata = pd.DataFrame(columns=columns)
+	path = pathlib.Path(folder_name)
+	for fname in sorted(path.glob('src*.xml')):
+		ssm = to_python(fname, sconv)
+		for grp in ssm:
+			for src in grp:
+				mmin, mmax = src.get_min_max_mag()
+				hdeps = [d[1] for d in src.hypocenter_distribution.data]
+				row = [src.source_id,
+					src.name,
+					src.tectonic_region_type,
+					src.upper_seismogenic_depth,
+					src.lower_seismogenic_depth,
+					src.magnitude_scaling_relationship,
+					mmin,
+					mmax,
+					hdeps,
+					src.mfd.kwargs['b_val']
+					]
+				sdata.loc[len(sdata.index)] = row
+	
+	print(tabulate(sdata, headers="keys", tablefmt="psql"))
\ No newline at end of file
diff --git a/openquake/mbi/wkf/analysis_nodal_plane.py b/openquake/mbi/wkf/analysis_nodal_plane.py
index c4a3fedee..686e79bb4 100755
--- a/openquake/mbi/wkf/analysis_nodal_plane.py
+++ b/openquake/mbi/wkf/analysis_nodal_plane.py
@@ -5,8 +5,8 @@
 from openquake.wkf.seismicity.nodal_plane import process_gcmt_datafames
 
 
-def main(fname_folder, folder_out):
-    process_gcmt_datafames(fname_folder, folder_out)
+def main(fname_folder, folder_out, save_csv = False):
+    process_gcmt_datafames(fname_folder, folder_out, save_csv)
 
 
 main.fname_folder = 'Name of the folder with input files'
diff --git a/openquake/mbi/wkf/check_toml.py b/openquake/mbi/wkf/check_toml.py
new file mode 100644
index 000000000..dbd1a3652
--- /dev/null
+++ b/openquake/mbi/wkf/check_toml.py
@@ -0,0 +1,122 @@
+#!/usr/bin/env python
+# coding: utf-8
+# ------------------- The Model Building Toolkit ------------------------------
+#
+# Copyright (C) 2022 GEM Foundation
+#
+# This program is free software: you can redistribute it and/or modify it under
+# the terms of the GNU Affero General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# This program is distributed in the hope that it will be useful, but WITHOUT
+# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+# FOR A PARTICULAR PURPOSE.  See the GNU Affero General Public License for more
+# details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+# -----------------------------------------------------------------------------
+
+import toml
+import shutil
+import os
+from openquake.wkf.utils import get_list
+
+
+HERE = os.path.dirname(__file__)
+
+def main(fname_config: str, copy_loc: str, use: str = []):
+	""" Check config file contains all necessary parameters and make a working copy
+	
+	:param fname_config:
+		Location of config file describing model
+	:param copy_loc:
+		location to store an editable copy of the config file
+
+	"""
+	model = toml.load(fname_config)
+
+	if all(i in model for i in ('name', 'mmin', 'bin_width', 'rupture_mesh_spacing')):
+	    print("Checking model ", model.get('name'))
+	else: 
+		print("missing default data (name, mmin, bin_width or rupture_mesh_spacing not found)")
+		
+	# Check for declustering parameters
+	if 'declustering' in model:
+		declust_params = model['declustering']
+		if declust_params['method'] == 'Zaliapin':
+			if all(i in declust_params for i in ('fractal_dim', 'b_value', 'threshold', 'depth', 'output_nearest_neighbor_distances')):
+				print("All parameters specified for Zaliapin declustering")
+			else:
+				print("Check Zaliapin declustering parameters. Should contain: 'fractal_dim', 'b_value', 'threshold', 'depth', 'output_nearest_neighbor_distances'.")
+
+		elif declust_params['method'] == 'Reasenberg':
+			if all(i in declust_params for i in ('taumin', 'taumax', 'p', 'xk', 'xmeff', 'rfact', 'horiz_error', 'depth_error', 'interaction_formula', 'max_interaction_dist')):
+				print("All parameters specified for Reasenberg declustering")
+			else:
+				print(declust_params)
+				print("Check Reasenberg parameters. Should contain: 'taumin', 'taumax', 'p', 'xk', 'xmeff', 'rfact', 'horiz_error', 'depth_error', 'interaction_formula', 'max_interaction_dist'.")
+
+		elif declust_params['method'] == 'windowing':
+			if all(i in declust_params for i in ('time_distance_window', 'fs_time_prop')):
+				print("All parameters specified for windowing declustering")
+			else:
+				print(declust_params)
+				print("Check windowing parameters. Should contain: 'time_distance_window', 'fs_time_prop'.")
+
+		else: 
+			print("unrecognised declustering algorithm. Please choose from 'Zaliapin', 'Reasenberg' or 'windowing'" )
+
+
+	else:
+		print("No declustering algorithm found")
+
+	# Check smoothing parameters
+	if 'smoothing' in model: 
+		gauss = False
+		adap = False
+		smoothing_params = model['smoothing']
+		if all(i in smoothing_params for i in ('n_v', 'kernel', 'd_i_min', 'maxdist', 'h3res')):
+			print("Found parameters for adaptive smoothing")
+			adap = True
+		if all(i in smoothing_params for i in ('kernel_maximum_distance', 'kernel_smoothing')):
+			print("Found parameters for Gaussian smoothing")
+			gauss = True
+
+		if adap == False and gauss == False:
+			print("Smoothing paramaters missing. 'kernel_maximum_distance', 'kernel_smoothing' needed for Gaussian smoothing. 'n_v', 'kernel', 'd_i_min', 'maxdist', 'h3res' needed for adaptive smoothing. ")
+
+	else:
+		print("No smoothing parameters found")
+		
+	# check sources
+	if len(use) > 0:
+		use = get_list(use)
+		for src_id in use:
+		# Add sources specified by use
+			if src_id not in model['sources']:
+				model['sources'][src_id] = {}
+		# Should remove sources not in use, but for some reason this is stupidly difficult to find out how to do!
+		print("Configured ", len(use), "sources as specified by use argument")
+	    
+	elif 'sources' in model:
+		print("Found ", len(model['sources']), " sources in model config")
+
+	else:
+		print("No sources found. Please add some!")
+
+	# Check for msr section
+	if 'msr' not in model: 
+		print("No magnitude scaling relationships defined")
+
+	print("copying toml to ", copy_loc)
+	source = fname_config
+	destination = copy_loc
+	shutil.copy(source, destination)
+	
+main.fname_config = 'location of toml file for source model'
+main.copy_loc = 'location to store an editable copy of the config'
+
+if __name__ == '__main__':
+    sap.run(main)
diff --git a/openquake/mbi/wkf/create_smoothing_per_zone.py b/openquake/mbi/wkf/create_smoothing_per_zone.py
index f5d9c9960..ef9d79c9b 100755
--- a/openquake/mbi/wkf/create_smoothing_per_zone.py
+++ b/openquake/mbi/wkf/create_smoothing_per_zone.py
@@ -5,9 +5,9 @@
 from openquake.wkf.seismicity.smoothing import create_smoothing_per_zone
 
 
-def main(fname_points: str, fname_polygons: str, folder_out: str='/tmp',
-         skip=[]):
-    create_smoothing_per_zone(fname_points, fname_polygons, folder_out, skip)
+def main(fname_points: str, fname_polygons: str, folder_out: str='/tmp',*,
+         skip=[], use = []):
+    create_smoothing_per_zone(fname_points, fname_polygons, folder_out, skip, use)
 
 
 main.fname_points = ".csv file created by the smoothing code"
diff --git a/openquake/mbi/wkf/create_subcatalogues_per_zone.py b/openquake/mbi/wkf/create_subcatalogues_per_zone.py
index 618341696..43ae0bac0 100755
--- a/openquake/mbi/wkf/create_subcatalogues_per_zone.py
+++ b/openquake/mbi/wkf/create_subcatalogues_per_zone.py
@@ -6,7 +6,7 @@
 
 
 def main(fname_polygons: str, fname_cat: str, folder_out: str, *,
-         source_ids: list=[]):
+         source_ids: str=[]):
     """
     Given a file (e.g. a shapefile) with a set of polygons and an earthquake
     catalogue, it creates a set of .csv files each one containing the
diff --git a/openquake/mbi/wkf/focal_mech_loc_plots.py b/openquake/mbi/wkf/focal_mech_loc_plots.py
new file mode 100644
index 000000000..f3fb5d22e
--- /dev/null
+++ b/openquake/mbi/wkf/focal_mech_loc_plots.py
@@ -0,0 +1,67 @@
+from obspy.imaging.beachball import beach
+from openquake.hmtk.parsers.catalogue.gcmt_ndk_parser import ParseNDKtoGCMT
+import pandas as pd
+import matplotlib.pyplot as plt
+from openquake.sub.utils import  mecclass
+from matplotlib.gridspec import GridSpec
+import numpy as np
+
+KAVERINA = {'N': 'blue',
+            'SS': 'green',
+            'R': 'red',
+            'N-SS': 'turquoise',
+            'SS-N': 'palegreen',
+            'R-SS': 'goldenrod',
+            'SS-R': 'yellow'}
+
+def focal_mech_loc_plots(fname, figsize = (15, 10), width = 0.5, size = 0.1):
+    """
+    Produce a figure consisting of:
+    	1) nodal planes plotted in space (lat/Lon) with Kaverina classification colours
+    	2) scatterplot of event Kaverina classificatons and magnitudes
+    	3) scatterplot of event strike vs rake, coloured by Kaverina classification
+    Please note that the 'width' parameter might need to be adjusted for different models 
+    """
+    
+    cmt_cat_zone = pd.read_csv(fname)
+    plungeb = cmt_cat_zone['plunge_b']
+    plungep = cmt_cat_zone['plunge_p']
+    plunget = cmt_cat_zone['plunge_t']
+    mclass = ['']*len(plunget)
+
+    for i in range(0, len(plungeb)):
+        mclass[i] = mecclass(plunget[i], plungeb[i], plungep[i])
+
+    cmt_cat_zone['class'] = mclass    
+    
+    mts = np.column_stack([cmt_cat_zone.strike1, cmt_cat_zone.dip1, cmt_cat_zone.rake1])
+    
+    fig = plt.figure(layout="constrained", figsize = figsize)
+    gs = GridSpec( 2, 3, figure=fig)
+    a0 = fig.add_subplot(gs[0:, :-1])
+    
+    a0.set_xlim(np.min(cmt_cat_zone['longitude']) - 0.1, np.max(cmt_cat_zone['longitude'])+ 0.1)
+    a0.set_ylim(np.min(cmt_cat_zone['latitude']) - 0.1, np.max(cmt_cat_zone['latitude']) + 0.1)
+
+    a0.margins(0.05)
+
+    idx = 0
+    for i in range(0, len(plungeb)):
+        bcc = beach(mts[idx],xy=(cmt_cat_zone['longitude'][idx], cmt_cat_zone['latitude'][idx]), width=width, linewidth=1, zorder=20, size=size, facecolor=KAVERINA[mclass[idx]])
+        bcc.set_alpha(0.5)
+        a0.add_collection(bcc)
+        idx += 1
+    
+    a1 = fig.add_subplot(gs[0, -1])
+
+    a1.scatter(cmt_cat_zone['class'], cmt_cat_zone['magnitude'], c=cmt_cat_zone['class'].map(KAVERINA))
+    a1.set_xlabel("Kaverina classification")
+    a1.set_ylabel("magnitude")
+    
+    a2 = fig.add_subplot(gs[1, -1])
+    a2.scatter(cmt_cat_zone['strike1'], cmt_cat_zone['rake1'], c=cmt_cat_zone['class'].map(KAVERINA), s = 1, alpha = 0.5)
+    a2.set_xlabel("strike")
+    a2.set_ylabel("rake")
+
+    fig.suptitle("Zone nodal plane distribution")
+    plt.show()
diff --git a/openquake/mbi/wkf/set_defaults.py b/openquake/mbi/wkf/set_defaults.py
index 02c146949..d72e7cc5a 100755
--- a/openquake/mbi/wkf/set_defaults.py
+++ b/openquake/mbi/wkf/set_defaults.py
@@ -7,8 +7,13 @@
 from openquake.wkf.utils import get_list
 
 
-def main(fname_conf: str, fname_defaults: str):
+def main(fname_conf: str, fname_defaults: str, *, use: str = [], skip: str = []):
 
+    if len(use) > 0:
+        use = get_list(use)
+    if len(skip) > 0:
+        skip = get_list(skip)
+        
     # Parsing config
     model = toml.load(fname_conf)
     output = copy.copy(model)
@@ -17,10 +22,12 @@ def main(fname_conf: str, fname_defaults: str):
     defaults = toml.load(fname_defaults)
 
     # Adding fields
-    if 'sources' not in output:
-        output['sources'] = {}
+    #if 'sources' not in output:
+    #    output['sources'] = {}
 
     for src_id in defaults['sources']:
+        if (len(use) and src_id not in use) or (src_id in skip):
+            continue
         if src_id not in output['sources']:
             output['sources'][src_id] = {}
         for key in defaults['sources'][src_id]:
diff --git a/openquake/mbi/wkf/set_gr_params.py b/openquake/mbi/wkf/set_gr_params.py
index 2348c7d32..b8bf7a30e 100755
--- a/openquake/mbi/wkf/set_gr_params.py
+++ b/openquake/mbi/wkf/set_gr_params.py
@@ -45,6 +45,8 @@ def set_gr_params(fname_conf: str, use: str = "*", method: str = "weichert",
     for src_id in model['sources']:
         if exclude is not None and src_id in exclude:
             continue
+        if use != "*" and src_id not in use:
+            continue
         else:
             print("src_id:", src_id, " ", method)
             if use == "*" or src_id in use:
diff --git a/openquake/mbi/wkf/set_h3_to_zones.py b/openquake/mbi/wkf/set_h3_to_zones.py
index 94bca7f0a..cb9cf1cbd 100755
--- a/openquake/mbi/wkf/set_h3_to_zones.py
+++ b/openquake/mbi/wkf/set_h3_to_zones.py
@@ -5,12 +5,13 @@
 from openquake.wkf.h3.zones import discretize_zones_with_h3_grid
 
 
-def main(h3_level: str, fname_poly: str, folder_out: str, *, use=[]):
+
+def main(h3_level: str, fname_poly: str, folder_out: str,*, use: str = []):
     """
     Given a set of polygons, using H3 this creates grids representing in a
     'discrete' sense the original geometries
     """
-    discretize_zones_with_h3_grid(h3_level, fname_poly, folder_out, use)
+    discretize_zones_with_h3_grid(h3_level, fname_poly, folder_out, use = use)
 
 
 descr = 'The level of the H3 grid'
diff --git a/openquake/mbt/tools/mfd.py b/openquake/mbt/tools/mfd.py
index e20aa4c0a..9b5bf93a7 100644
--- a/openquake/mbt/tools/mfd.py
+++ b/openquake/mbt/tools/mfd.py
@@ -4,7 +4,6 @@
 
 import scipy
 import numpy as np
-
 from scipy.stats import truncnorm
 
 from openquake.hazardlib.mfd import (TruncatedGRMFD, EvenlyDiscretizedMFD,
@@ -13,7 +12,7 @@
 from openquake.hazardlib.mfd.youngs_coppersmith_1985 import (
         YoungsCoppersmith1985MFD)
 
-log = True
+#log = True
 log = False
 
 
@@ -226,11 +225,40 @@ def get_evenlyDiscretizedMFD_from_multiMFD(mfd, bin_width=None):
             min_m = min_mag[0] if len(min_mag) == 1 else min_mag[i]
             if i == 0:
                 emfd = EEvenlyDiscretizedMFD(min_m, binw, occ)
+                
             else:
                 tmfd = EEvenlyDiscretizedMFD(min_m, binw, occ)
                 emfd.stack(tmfd)
+                 
+    elif mfd.kind == 'truncGutenbergRichterMFD':
+        aval = mfd.kwargs['a_val']
+        bval = mfd.kwargs['b_val']
+        min_mag = mfd.kwargs['min_mag']
+        max_mag = mfd.kwargs['max_mag']
+        binw = mfd.kwargs['bin_width'][0]
+        # take max mag here so we don't have to rescale the FMD later
+        max_m = np.max(max_mag) 
+        min_m = min_mag[0] if len(min_mag) == 1 else np.min(min_mag) 
+
+        for i in range(mfd.size):
+            bgr = bval[0] if len(bval) == 1 else bval[i]
+            agr = aval[0] if len(aval) == 1 else aval[i]
+            
+            left = np.arange(min_m, max_m, binw)
+            rates = 10.**(agr-bgr*left)-10.**(agr-bgr*(left+binw))
+            
+            if i == 0:
+                emfd = EEvenlyDiscretizedMFD(min_m+binw/2.,
+                                binw,
+                                list(rates))
+            else:
+                tmfd = EEvenlyDiscretizedMFD(min_m+binw/2.,
+                                binw,
+                                list(rates))
+                emfd.stack(tmfd)
+    
     else:
-        raise ValueError('Unsupported MFD type')
+        raise ValueError('Unsupported MFD type ', mfd.kind)
 
     return emfd
 
@@ -275,6 +303,8 @@ def stack(self, imfd):
 
         mfd1 = self
         bin_width = self.bin_width
+        
+        
         #
         # check bin width of the MFD to be added
         if (isinstance(mfd2, EvenlyDiscretizedMFD) and
@@ -309,6 +339,7 @@ def stack(self, imfd):
         #
         # mfd1 MUST be the one with the mininum minimum magnitude
         if mfd1.min_mag > mfd2.min_mag:
+            print("minimum magnitudes are different")
             if log:
                 print('SWAPPING')
             tmp = mfd2
@@ -324,8 +355,8 @@ def stack(self, imfd):
             tmpmag += bin_width
 
         rates = list(np.zeros(len(mfd1.occurrence_rates)))
-        mags = list(mfd1.min_mag+np.arange(len(rates))*bin_width)
-
+        mags = list(mfd1.min_mag+np.arange(len(rates))*bin_width)	
+	
         # Add to the rates list the occurrences included in the mfd with the
         # lowest minimum magnitude
         for idx, occ in enumerate(mfd1.occurrence_rates):
@@ -352,12 +383,13 @@ def stack(self, imfd):
             try:
                 if len(rates) > idx+delta:
                     assert abs(mag - mags[idx+delta]) < 1e-5
+                    
             except:
                 print('mag:     :', mag)
                 print('mag rates:', mags[idx+delta])
                 print('delta    :', delta)
                 print('diff     :', abs(mag - mags[idx+delta]))
-                raise ValueError('Staking wrong bins')
+                raise ValueError('Stacking wrong bins')
 
             if log:
                 print(idx, idx+delta, len(mfd2.occurrence_rates), len(rates))
@@ -552,7 +584,6 @@ def mfd_upsample(bin_width, mfd):
     # assigning occurrences
     dlt = bin_width * 1e-5
     for mag, occ in mfd.get_annual_occurrence_rates():
-        print(mag, occ)
         #
         # find indexes of lower bin limits lower than mag
         idx = np.nonzero(mag+dlt-mfd.bin_width/2 > nocc[:, 1])[0]
@@ -571,7 +602,6 @@ def mfd_upsample(bin_width, mfd):
             idxb = np.amin(idx)
         #
         #
-        print(idxa, idxb)
         if idxb is not None and idxa == idxb:
             nocc[idxa, 3] += occ
         else:
diff --git a/openquake/mbt/tools/model_building/plt_mfd.py b/openquake/mbt/tools/model_building/plt_mfd.py
index 229b431a6..280e29446 100755
--- a/openquake/mbt/tools/model_building/plt_mfd.py
+++ b/openquake/mbt/tools/model_building/plt_mfd.py
@@ -25,7 +25,7 @@ def _compute_mfd(cat, compl_table, mwid):
     weichert_config = {'magnitude_interval': mwid,
                        'reference_magnitude': 0.0}
     weichert = Weichert()
-    bval_wei, sigmab, aval_wei, sigmaa = weichert.calculate(cat,
+    bval_wei, sigmab, aval_wei, sigmaa = weichert.calc(cat,
                                                             weichert_config,
                                                             compl_table)
     #
diff --git a/openquake/sub/misc/edge.py b/openquake/sub/misc/edge.py
index cce37b520..8d846ed65 100644
--- a/openquake/sub/misc/edge.py
+++ b/openquake/sub/misc/edge.py
@@ -22,7 +22,7 @@
 from mpl_toolkits.mplot3d import Axes3D
 
 TOL = 0.7
-
+#TOL = 1.5
 
 def line_between_two_points(pnt1, pnt2):
     """
diff --git a/openquake/sub/slab/rupture.py b/openquake/sub/slab/rupture.py
index 34e561a36..5131eb9a8 100755
--- a/openquake/sub/slab/rupture.py
+++ b/openquake/sub/slab/rupture.py
@@ -41,8 +41,8 @@
 from openquake.man.checks.catalogue import load_catalogue
 from openquake.wkf.utils import create_folder
 
-#PLOTTING = True
-PLOTTING = False
+PLOTTING = True
+#PLOTTING = False
 
 
 def get_catalogue(cat_pickle_fname, treg_filename=None, label='', 
diff --git a/openquake/wkf/catalogue.py b/openquake/wkf/catalogue.py
index f5417c8d4..86fe8447f 100644
--- a/openquake/wkf/catalogue.py
+++ b/openquake/wkf/catalogue.py
@@ -36,7 +36,7 @@
 from typing import Type
 from shapely.geometry import Point
 from openquake.hmtk.seismicity.catalogue import Catalogue
-from openquake.wkf.utils import create_folder
+from openquake.wkf.utils import create_folder, get_list
 from openquake.hmtk.parsers.catalogue.gcmt_ndk_parser import ParseNDKtoGCMT
 
 
@@ -121,7 +121,7 @@ def from_df(df, end_year=None) -> Type[Catalogue]:
 
 
 def create_subcatalogues(fname_polygons: str, fname_cat: str, folder_out: str,
-                         source_ids: list = []):
+                         source_ids: str = []):
     """
     Given a catalogue and a gis-file with polygons (e.g. shapefile or
     .geojson), this code creates for each polygon a subcatalogue with the
@@ -141,7 +141,10 @@ def create_subcatalogues(fname_polygons: str, fname_cat: str, folder_out: str,
 
     # Create output folder
     create_folder(folder_out)
-
+    
+    if len(source_ids) > 0:
+        source_ids = get_list(source_ids)
+        
     # Create geodataframe with the catalogue
     df = pd.read_csv(fname_cat)
     gdf = gpd.GeoDataFrame(df, crs='epsg:4326', geometry=[Point(xy) for xy
@@ -149,21 +152,27 @@ def create_subcatalogues(fname_polygons: str, fname_cat: str, folder_out: str,
 
     # Read polygons
     polygons_gdf = gpd.read_file(fname_polygons)
-
-    # Select point in polygon
-    columns = ['eventID', 'year', 'month', 'day', 'magnitude', 'longitude',
-               'latitude', 'depth']
-
+    
+    # check there are no duplicate ids
+    assert(polygons_gdf['id'].is_unique)
+    
+    # explode for any multipolygons
+    polygons_gdf = polygons_gdf.explode()
+    
     # Iterate over sources
+    # actually this should maybe iterate over *sources* rather than the polygons, because then we could 
+    # handle the multipolygons more smoothly
     out_fnames = []
     for idx, poly in polygons_gdf.iterrows():
 
         if len(source_ids) > 0 and poly.id not in source_ids:
             continue
-
+            
         df = pd.DataFrame({'Name': [poly.id], 'Polygon': [poly.geometry]})
         gdf_poly = gpd.GeoDataFrame(df, geometry='Polygon', crs='epsg:4326')
         within = gpd.sjoin(gdf, gdf_poly, op='within')
+        
+        # check if there are other polygons for this source
         # Create output file
         if isinstance(poly.id, int):
             fname = f'subcatalogue_zone_{poly.id:d}.csv'
@@ -171,7 +180,7 @@ def create_subcatalogues(fname_polygons: str, fname_cat: str, folder_out: str,
             fname = f'subcatalogue_zone_{poly.id}.csv'
         out_fname = os.path.join(folder_out, fname)
         out_fnames.append(out_fname)
-        within.to_csv(out_fname, index=False, columns=columns)
+        within.to_csv(out_fname, index=False)
 
     return out_fnames
 
diff --git a/openquake/wkf/h3/zones.py b/openquake/wkf/h3/zones.py
index bc0011654..5cce90cde 100755
--- a/openquake/wkf/h3/zones.py
+++ b/openquake/wkf/h3/zones.py
@@ -34,8 +34,8 @@
 
 
 def discretize_zones_with_h3_grid(h3_level: str, fname_poly: str,
-                                  folder_out: str, use: str =[]):
-
+                                  folder_out: str, *, use: str = []):
+    
     h3_level = int(h3_level)
     create_folder(folder_out)
     tmp = "mapping_h{:d}.csv".format(h3_level)
@@ -51,21 +51,28 @@ def discretize_zones_with_h3_grid(h3_level: str, fname_poly: str,
     # Select point in polygon
     fout = open(fname_out, 'w')
     for idx, poly in polygons_gdf.iterrows():
+    
+        if len(use) > 0 and poly.id not in use:
+            continue
+
         tmps = shapely.geometry.mapping(poly.geometry)
         geojson_poly = eval(json.dumps(tmps))
         if geojson_poly['type'] == 'MultiPolygon':
             from shapely.geometry import shape, mapping
             # Check that there are no polygons inside
             multipoly = shape(geojson_poly)
-            assert len(multipoly.geoms) == 1
-            geojson_poly = mapping(multipoly.geoms[0])
+            try:
+                assert len(multipoly.geoms) == 1
+                geojson_poly = mapping(multipoly.geoms[0])
+            except:
+                print(poly.id)
 
         # Revert the positions of lons and lats
-        coo = [[c[1], c[0]] for c in geojson_poly['coordinates'][0]]
-        geojson_poly['coordinates'] = [coo]
+        #coo = [[c[1], c[0]] for c in geojson_poly['coordinates'][0]]
+        #geojson_poly['coordinates'] = [coo]
 
         # Discretizing
-        hexagons = list(h3.polyfill(geojson_poly, h3_level))
+        hexagons = list(h3.polyfill(geojson_poly, h3_level, geo_json_conformant=True))
         for hxg in hexagons:
             if isinstance(poly.id, str):
                 fout.write("{:s},{:s}\n".format(hxg, poly.id))
diff --git a/openquake/wkf/mfd.py b/openquake/wkf/mfd.py
index 81ba089a7..a2e8543c2 100644
--- a/openquake/wkf/mfd.py
+++ b/openquake/wkf/mfd.py
@@ -37,7 +37,7 @@
     get_evenlyDiscretizedMFD_from_truncatedGRMFD)
 
 
-def check_mfds(fname_input_pattern: str, fname_config: str, *,
+def check_mfds(fname_input_pattern: str, fname_config: str = None,  *,
                src_id: str = None):
     """
     Given a set of .xml files and a configuration file with GR params, this
@@ -45,12 +45,12 @@ def check_mfds(fname_input_pattern: str, fname_config: str, *,
     configuration file. The ID of the source if not provided is taken from the
     name of the files (i.e., last label preceded by `_`)
     """
-
     for fname in sorted(glob(fname_input_pattern)):
 
-        if src_id is None:
-            src_id = _get_src_id(fname)
-        model = toml.load(fname_config)
+        srcs = _get_src_id(fname)
+        if src_id is not None and srcs not in src_id:
+        	continue
+        
 
         binw = 0.1
         sourceconv = SourceConverter(investigation_time=1.0,
@@ -60,7 +60,6 @@ def check_mfds(fname_input_pattern: str, fname_config: str, *,
         ssm = to_python(fname, sourceconv)
 
         for grp in ssm:
-
             for i, src in enumerate(grp):
                 if i == 0:
                     nmfd = EEvenlyDiscretizedMFD.from_mfd(src.mfd, binw)
@@ -68,20 +67,25 @@ def check_mfds(fname_input_pattern: str, fname_config: str, *,
                     ged = get_evenlyDiscretizedMFD_from_truncatedGRMFD
                     tmfd = ged(src.mfd, nmfd.bin_width)
                     nmfd.stack(tmfd)
-
+            
             occ = np.array(nmfd.get_annual_occurrence_rates())
 
-            bgr = model["sources"][src_id]["bgr_weichert"]
-            agr = model["sources"][src_id]["agr_weichert"]
-
-            tmp = occ[:, 0] - binw
-            mfd = 10.0**(agr-bgr*tmp[:-1])-10.0**(agr-bgr*(tmp[:-1]+binw))
+            if fname_config:
+                model = toml.load(fname_config)
+                bgr = model["sources"][srcs]["bgr"]
+                agr = model["sources"][srcs]["agr"]
+                tmp = occ[:, 0] - binw
+                mfd = 10.0**(agr-bgr*tmp[:-1])-10.0**(agr-bgr*(tmp[:-1]+binw))
 
             _ = plt.figure(figsize=(8, 6))
-            plt.plot(occ[:, 0], occ[:, 1], 'o')
-            plt.plot(tmp[:-1]+binw/2, mfd, 'x')
+            plt.plot(occ[:, 0], occ[:, 1], 'o', 
+                label = 'model occurrence rate for source')
+            if fname_config:
+                plt.plot(tmp[:-1]+binw/2, mfd, 'x', 
+                    label = ('config: a = ', agr, ", b = ", bgr))
             plt.title(fname)
             plt.xlabel('Magnitude')
             plt.ylabel('Annual occurrence rate')
             plt.yscale('log')
+            plt.legend()
             plt.show()
diff --git a/openquake/wkf/seismicity/smoothing.py b/openquake/wkf/seismicity/smoothing.py
index deaa87890..bd15fccbe 100755
--- a/openquake/wkf/seismicity/smoothing.py
+++ b/openquake/wkf/seismicity/smoothing.py
@@ -5,7 +5,7 @@
 import pandas as pd
 import geopandas as gpd
 import matplotlib.pyplot as plt
-from openquake.wkf.utils import create_folder
+from openquake.wkf.utils import create_folder, get_list
 from openquake.baselib import sap
 from shapely.geometry import Point
 from shapely import wkt
@@ -17,7 +17,7 @@ def main(fname_points: str, fname_polygons: str, folder_out: str='/tmp',
 
 
 def create_smoothing_per_zone(fname_points: str, fname_polygons: str,
-                              folder_out: str='/tmp', skip=[]):
+                              folder_out: str='/tmp', *, skip=[], use = []):
     """
     Creates subsets of points, one for each of the polygons included in the
     `fname_polygons` shapefile. The attibute table must have an 'id' attribute.
@@ -33,6 +33,9 @@ def create_smoothing_per_zone(fname_points: str, fname_polygons: str,
     """
 
     create_folder(folder_out)
+    
+    if len(use) > 0:
+        use=get_list(use)
 
     # Create a geodataframe with the point sources
     df = pd.read_csv(fname_points)
@@ -47,6 +50,9 @@ def create_smoothing_per_zone(fname_points: str, fname_polygons: str,
 
         if poly.id in skip:
             continue
+            
+        if len(use) > 0 and poly.id not in use:
+            continue
 
         # Create a geodataframe with the polygon in question
         df = pd.DataFrame({'Name': [poly.id], 'Polygon': [poly.geometry]})