diff --git a/.gitignore b/.gitignore
index f33a818cf..92c6c9253 100644
--- a/.gitignore
+++ b/.gitignore
@@ -90,4 +90,9 @@ ENV/
 # Rope project settings
 .ropeproject
 
+# debug_data
+pcmdi_metrics/monsoon_wang/*.nc
+pcmdi_metrics/monsoon_wang/debug_regions_specs.py
+pcmdi_metrics/monsoon_wang/test_param.py
+
 test/
diff --git a/pcmdi_metrics/io/default_regions_define.py b/pcmdi_metrics/io/default_regions_define.py
index e24ced93f..b360bf090 100755
--- a/pcmdi_metrics/io/default_regions_define.py
+++ b/pcmdi_metrics/io/default_regions_define.py
@@ -1,9 +1,9 @@
-import xcdat as xc
 from typing import Union
 
 import xarray as xr
+import xcdat as xc
 
-#from pcmdi_metrics.io import da_to_ds, get_longitude, select_subset
+# from pcmdi_metrics.io import da_to_ds, get_longitude, select_subset
 from .xcdat_dataset_io import da_to_ds, get_longitude, select_subset
 
 
@@ -51,7 +51,7 @@ def load_regions_specs():
         # South American Monsoon
         "SAMM": {"domain": {"latitude": (-45.0, 0.0), "longitude": (240.0, 330.0)}},
         # North African Monsoon
-        #"NAFM": {"domain": {"latitude": (0.0, 45.0), "longitude": (310.0, 60.0)}},
+        # "NAFM": {"domain": {"latitude": (0.0, 45.0), "longitude": (310.0, 60.0)}},
         "NAFM": {"domain": {"latitude": (0.0, 45.0), "longitude": (-50, 60.0)}},
         # South African Monsoon
         "SAFM": {"domain": {"latitude": (-45.0, 0.0), "longitude": (0.0, 90.0)}},
@@ -77,7 +77,7 @@ def load_regions_specs():
     return regions_specs
 
 
-#def region_subset(ds, regions_specs, region=None):
+# def region_subset(ds, regions_specs, region=None):
 #    """
 #    d: xarray.Dataset
 #    regions_specs: dict
@@ -132,6 +132,7 @@ def load_regions_specs():
 #
 #    return ds
 
+
 def region_subset(
     ds: Union[xr.Dataset, xr.DataArray],
     region: str,
diff --git a/pcmdi_metrics/io/xcdat_dataset_io.py b/pcmdi_metrics/io/xcdat_dataset_io.py
index ed2bb14c2..2cdc4a194 100644
--- a/pcmdi_metrics/io/xcdat_dataset_io.py
+++ b/pcmdi_metrics/io/xcdat_dataset_io.py
@@ -3,7 +3,6 @@
 import xarray as xr
 import xcdat as xc
 
-
 # Internal function
 
 
diff --git a/pcmdi_metrics/monsoon_sperber/driver_monsoon_sperber.py b/pcmdi_metrics/monsoon_sperber/driver_monsoon_sperber.py
index baada145e..743be7e0b 100644
--- a/pcmdi_metrics/monsoon_sperber/driver_monsoon_sperber.py
+++ b/pcmdi_metrics/monsoon_sperber/driver_monsoon_sperber.py
@@ -46,17 +46,18 @@
 from shutil import copyfile
 
 import matplotlib
-matplotlib.use('Agg')
-#import matplotlib.pyplot as plt
-from matplotlib import pyplot as plt
+
+matplotlib.use("Agg")
 import numpy as np
 import pandas as pd
 import xarray as xr
 import xcdat as xc
+# import matplotlib.pyplot as plt
+from matplotlib import pyplot as plt
 
 import pcmdi_metrics
 from pcmdi_metrics import resources
-from pcmdi_metrics.io import load_regions_specs, region_subset
+from pcmdi_metrics.io import load_regions_specs, region_subset, xcdat_open
 from pcmdi_metrics.mean_climate.lib import pmp_parser
 from pcmdi_metrics.monsoon_sperber.lib import (
     AddParserArgument,
@@ -67,7 +68,6 @@
     sperber_metrics,
 )
 from pcmdi_metrics.utils import create_land_sea_mask, fill_template
-from pcmdi_metrics.io import xcdat_open
 
 
 def tree():
@@ -92,7 +92,7 @@ def pick_year_last_day(ds):
 # =================================================
 # Hard coded options... will be moved out later
 # -------------------------------------------------
-#list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
+# list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
 list_monsoon_regions = ["AIR", "Sahel"]
 # list_monsoon_regions = ["all"]
 
@@ -293,7 +293,16 @@ def pick_year_last_day(ds):
                 modpath(model=model, exp=exp, realization=realization, variable=var)
             )
             if debug:
-                print("model: ", model, "   exp: ", exp, "  realization: ", realization, "  variable: ", var)
+                print(
+                    "model: ",
+                    model,
+                    "   exp: ",
+                    exp,
+                    "  realization: ",
+                    realization,
+                    "  variable: ",
+                    var,
+                )
                 print("debug: model_path_list: ", model_path_list)
             # land fraction
             model_lf_path = modpath_lf(model=model)
@@ -314,18 +323,18 @@ def pick_year_last_day(ds):
                     ds_lf = xcdat_open(model_lf_path)
                 except Exception:
                     ds_lf = None
-        
-        if not ds_lf: 
+
+        if not ds_lf:
             lf_array = create_land_sea_mask(ds_lf, method="pcmdi")
             ds_lf = lf_array.to_dataset().compute()
             ds_lf = ds_lf.rename_vars({"lsmask": "sftlf"})
 
-	#  use pcmdi mask
-	#  lf_array = create_land_sea_mask(ds_lf, method="pcmdi")
-	#  ds_lf = lf_array.to_dataset().compute()
-	#  ds_lf = ds_lf.rename_vars({"lsmask": "sftlf"})
+        #  use pcmdi mask
+        #  lf_array = create_land_sea_mask(ds_lf, method="pcmdi")
+        #  ds_lf = lf_array.to_dataset().compute()
+        #  ds_lf = ds_lf.rename_vars({"lsmask": "sftlf"})
 
-        if model in [ "EC-EARTH" ]: #, "BNU-ESM" ]:
+        if model in ["EC-EARTH"]:  # , "BNU-ESM" ]:
             ds_lf = ds_lf.isel(lat=slice(None, None, -1))
         lf = ds_lf.sftlf.sel(lat=slice(-90, 90))  # land frac file must be global
 
@@ -350,10 +359,9 @@ def pick_year_last_day(ds):
                 # Get time coordinate information
                 print("model_path =   ", model_path)
 
-
                 dc = xcdat_open(model_path, decode_times=True)
-                dc['time'].attrs['axis'] = 'T'
-                dc['time'].attrs['standard_name'] = 'time'
+                dc["time"].attrs["axis"] = "T"
+                dc["time"].attrs["standard_name"] = "time"
                 dc = xr.decode_cf(dc, decode_times=True)
                 dc = dc.bounds.add_missing_bounds("X")
                 dc = dc.bounds.add_missing_bounds("Y")
@@ -363,7 +371,6 @@ def pick_year_last_day(ds):
                 c = xc.center_times(dc)
                 eday = pick_year_last_day(dc)
 
-
                 # Get starting and ending year and month
                 startYear = c.time.values[0].year
                 startMonth = c.time.values[0].month
@@ -493,7 +500,6 @@ def pick_year_last_day(ds):
                         d.values = d.values * 86400.0
                         d["units"] = units
 
-
                     # variable for over land only
                     d_land = model_land_only(model, d, lf, debug=debug)
 
@@ -523,7 +529,6 @@ def pick_year_last_day(ds):
                             d_sub_pr.values = d_sub_pr.values * 86400.0
                             d_sub_pr["units"] = units
 
-
                         else:
                             # land-only rainfall
 
@@ -544,11 +549,9 @@ def pick_year_last_day(ds):
                                 model, d_sub_pr, lf_sub, debug=debug
                             )
 
-
                             d_sub_pr.values = d_sub_pr.values * 86400.0
                             d_sub_pr["units"] = units
 
-
                         # Area average
 
                         ds_sub_pr = d_sub_pr.to_dataset().compute()
@@ -557,7 +560,6 @@ def pick_year_last_day(ds):
                         ds_sub_pr = ds_sub_pr.bounds.add_missing_bounds("Y")
                         ds_sub_pr = ds_sub_pr.bounds.add_missing_bounds("T")
 
-
                         if "lat_bnds" not in ds_sub_pr.variables:
                             lat_bnds = dc["lat_bnds"].sel(lat=ds_sub_pr["lat"])
                             ds_sub_pr["lat_bnds"] = lat_bnds
@@ -567,8 +569,6 @@ def pick_year_last_day(ds):
                         ).compute()
                         d_sub_aave = ds_sub_aave.pr
 
-
-
                         if debug:
                             print("debug: region:", region)
 
@@ -603,7 +603,6 @@ def pick_year_last_day(ds):
 
                                 d_sub_aave = xr.concat([part1, part2], dim="time")
 
-
                                 if debug:
                                     print(
                                         "debug: ",
@@ -638,7 +637,6 @@ def pick_year_last_day(ds):
                             coords={"time": time_coords},
                         )
 
-
                         if debug:
                             print(
                                 "debug: pentad_time_series length: ",
@@ -648,14 +646,14 @@ def pick_year_last_day(ds):
                         # Keep pentad time series length in consistent
                         ref_length = int(365 / n)
                         if len(pentad_time_series) < ref_length:
-
                             pentad_time_series = pentad_time_series.interp(
-                                time=pd.date_range(time_coords[0], time_coords[-1], periods=ref_length)
+                                time=pd.date_range(
+                                    time_coords[0], time_coords[-1], periods=ref_length
+                                )
                             )
 
                             time_coords = pentad_time_series.coords["time"]
 
-
                         pentad_time_series_cumsum = np.cumsum(pentad_time_series)
                         pentad_time_series = xr.DataArray(
                             pentad_time_series,
@@ -672,7 +670,6 @@ def pick_year_last_day(ds):
                         pentad_time_series_cumsum.attrs["units"] = str(d.units.values)
                         pentad_time_series_cumsum.coords["time"] = time_coords
 
-
                         if nc_out:
                             # Archive individual year time series in netCDF file
                             pentad_time_series.to_netcdf(file_path, mode="a")
@@ -819,9 +816,7 @@ def pick_year_last_day(ds):
 
                         # obs
                         ax[region].plot(
-                            np.array(
-                                dict_obs_composite[reference_data_name][region]
-                            ),
+                            np.array(dict_obs_composite[reference_data_name][region]),
                             c="blue",
                             label=reference_data_name,
                         )
@@ -836,9 +831,9 @@ def pick_year_last_day(ds):
                             ax[region].axvline(
                                 x=idx,
                                 ymin=0,
-                                ymax=dict_obs_composite[reference_data_name][
-                                    region
-                                ][idx].item(),
+                                ymax=dict_obs_composite[reference_data_name][region][
+                                    idx
+                                ].item(),
                                 c="blue",
                                 ls="--",
                             )
diff --git a/pcmdi_metrics/monsoon_sperber/lib/argparse_functions.py b/pcmdi_metrics/monsoon_sperber/lib/argparse_functions.py
index 6fa42d61a..a51cc81a3 100644
--- a/pcmdi_metrics/monsoon_sperber/lib/argparse_functions.py
+++ b/pcmdi_metrics/monsoon_sperber/lib/argparse_functions.py
@@ -54,7 +54,9 @@ def AddParserArgument(P):
         "--meyear", dest="meyear", type=int, help="End year for model data set"
     )
     P.add_argument("--modnames", type=str, default=None, help="List of models")
-    P.add_argument("--list_monsoon_regions", type=str, default=None, help="List of regions")
+    P.add_argument(
+        "--list_monsoon_regions", type=str, default=None, help="List of regions"
+    )
     P.add_argument(
         "-r",
         "--realization",
diff --git a/pcmdi_metrics/monsoon_sperber/lib/divide_chunks.py b/pcmdi_metrics/monsoon_sperber/lib/divide_chunks.py
index efd031175..8a215f4f8 100644
--- a/pcmdi_metrics/monsoon_sperber/lib/divide_chunks.py
+++ b/pcmdi_metrics/monsoon_sperber/lib/divide_chunks.py
@@ -30,7 +30,7 @@ def divide_chunks_advanced(data, n, debug=False):
     day = day.values
     calendar = "gregorian"
     if debug:
-        #print("month = ", month, "day = ", day)
+        # print("month = ", month, "day = ", day)
         print("debug: first day of year is " + str(month) + "/" + str(day))
     if month not in [1, 7] or day != 1:
         sys.exit(
diff --git a/pcmdi_metrics/monsoon_sperber/lib/model_land_only.py b/pcmdi_metrics/monsoon_sperber/lib/model_land_only.py
index 2dd47cac2..b1319fb65 100644
--- a/pcmdi_metrics/monsoon_sperber/lib/model_land_only.py
+++ b/pcmdi_metrics/monsoon_sperber/lib/model_land_only.py
@@ -9,7 +9,7 @@ def model_land_only(model, model_timeseries, lf, debug=False):
     # - - - - - - - - - - - - - - - - - - - - - - - - -
 
     if debug:
-        #plot_map(model_timeseries[0], "_".join(["test", model, "beforeMask.png"]))
+        # plot_map(model_timeseries[0], "_".join(["test", model, "beforeMask.png"]))
         print("debug: plot for beforeMask done")
 
     # Check land fraction variable to see if it meet criteria
@@ -33,7 +33,7 @@ def model_land_only(model, model_timeseries, lf, debug=False):
             model_timeseries_masked = model_timeseries.where(lf > 90)
 
     if debug:
-        #plot_map(model_timeseries_masked[0], "_".join(["test", model, "afterMask.png"]))
+        # plot_map(model_timeseries_masked[0], "_".join(["test", model, "afterMask.png"]))
         print("debug: plot for afterMask done")
 
     return model_timeseries_masked
diff --git a/pcmdi_metrics/monsoon_sperber/param/Bo_param.py b/pcmdi_metrics/monsoon_sperber/param/Bo_param.py
index 5ec7d9522..ca0d9bd6e 100644
--- a/pcmdi_metrics/monsoon_sperber/param/Bo_param.py
+++ b/pcmdi_metrics/monsoon_sperber/param/Bo_param.py
@@ -14,9 +14,9 @@
 # -------------------------------------------------
 update_json = False
 debug = False
-#debug = True
+# debug = True
 
-#list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
+# list_monsoon_regions = ["AIR", "AUS", "Sahel", "GoG", "NAmo", "SAmo"]
 list_monsoon_regions = ["AUS"]
 # =================================================
 # Observation
@@ -41,7 +41,7 @@
 modpath = "/work/lee1043/ESGF/xmls/cmip5/historical/day/pr/cmip5.%(model).%(exp).%(realization).day.pr.xml"
 modpath_lf = "/work/lee1043/ESGF/xmls/cmip5/historical/fx/sftlf/cmip5.%(model).historical.r0i0p0.fx.sftlf.xml"
 
-#/p/css03/scratch/published-older/cmip5/output1/CSIRO-BOM/ACCESS1-0/historical/day/atmos/day/r1i1p1/v4/pr/pr_day_ACCESS1-0_historical_r1i1p1_19750101-19991231.nc
+# /p/css03/scratch/published-older/cmip5/output1/CSIRO-BOM/ACCESS1-0/historical/day/atmos/day/r1i1p1/v4/pr/pr_day_ACCESS1-0_historical_r1i1p1_19750101-19991231.nc
 
 # modnames = ['ACCESS1-0', 'ACCESS1-3', 'BCC-CSM1-1', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-BGC', 'CESM1-CAM5', 'CESM1-FASTCHEM', 'CMCC-CESM', 'CMCC-CM', 'CMCC-CMS', 'CNRM-CM5', 'CSIRO-Mk3-6-0', 'EC-EARTH', 'FGOALS-g2', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadGEM2-AO', 'HadGEM2-CC', 'HadGEM2-ES', 'INMCM4', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC-ESM', 'MIROC-ESM-CHEM', 'MIROC4h', 'MIROC5', 'MPI-ESM-MR', 'MPI-ESM-P', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M']  # noqa
 
@@ -60,7 +60,7 @@
 # =================================================
 # Output
 # -------------------------------------------------
-#pmprdir = "/p/user_pub/pmp/pmp_results/pmp_v1.1.2"
+# pmprdir = "/p/user_pub/pmp/pmp_results/pmp_v1.1.2"
 pmprdir = "/p/user_pub/climate_work/dong12/PMP_result/"
 case_id = "{:v%Y%m%d}".format(datetime.datetime.now())
 
diff --git a/pcmdi_metrics/monsoon_wang/basic_monsoon_wang_param.py b/pcmdi_metrics/monsoon_wang/basic_monsoon_wang_param.py
index 728afac2d..2755464b1 100644
--- a/pcmdi_metrics/monsoon_wang/basic_monsoon_wang_param.py
+++ b/pcmdi_metrics/monsoon_wang/basic_monsoon_wang_param.py
@@ -1,38 +1,57 @@
-import os
-
 #
 #  OPTIONS ARE SET BY USER IN THIS FILE AS INDICATED BELOW BY:
 #
 #
 
 # LIST OF MODEL VERSIONS TO BE TESTED
-modnames = ['CanCM4']
-modnames = ['IPSL-CM6A-LR']
-#modnames = ['ACCESS-CM2', 'ACCESS-ESM1-5', 'AWI-CM-1-1-MR', 'AWI-ESM-1-1-LR', 'BCC-CSM2-MR', 'BCC-ESM1', 'CAMS-CSM1-0', 'CanESM5-1', 'CanESM5', 'CAS-ESM2-0', 'CESM2-FV2', 'CESM2', 'CESM2-WACCM-FV2', 'CESM2-WACCM', 'CIESM', 'CMCC-CM2-HR4', 'CMCC-CM2-SR5', 'CMCC-ESM2', 'E3SM-1-0', 'E3SM-1-1-ECA', 'E3SM-1-1', 'E3SM-2-0', 'EC-Earth3-AerChem', 'EC-Earth3-CC', 'EC-Earth3', 'EC-Earth3-Veg-LR', 'EC-Earth3-Veg', 'FGOALS-f3-L', 'FIO-ESM-2-0', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G-CC', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'GISS-E2-2-H', 'ICON-ESM-LR', 'INM-CM4-8', 'INM-CM5-0', 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR-INCA', 'IPSL-CM6A-LR', 'KACE-1-0-G', 'KIOST-ESM', 'MCM-UA-1-0', 'MIROC6', 'MPI-ESM-1-2-HAM', 'MPI-ESM1-2-HR', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorCPM1', 'NorESM2-LM', 'NorESM2-MM', 'SAM0-UNICON', 'TaiESM1']
 
 #  trouble model ICON-ESM-LR'
 
-
-modnames = ['ACCESS-CM2', 'ACCESS-ESM1-5', 'AWI-ESM-1-1-LR', 'BCC-CSM2-MR', 'BCC-ESM1', 'CAMS-CSM1-0', 'CanESM5', 'CAS-ESM2-0', 'CESM2', 'CESM2-WACCM', 'CIESM', 'CMCC-ESM2', 'E3SM-1-0', 'E3SM-1-1', 'E3SM-2-0', 'EC-Earth3', 'FGOALS-f3-L', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G', 'GISS-E2-1-H', 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR', 'KIOST-ESM', 'MIROC6', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorESM2-LM', 'TaiESM1']
-
-
-#modnames = [ 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR', 'KIOST-ESM', 'MIROC6', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorESM2-LM', 'TaiESM1']
-
+modnames = [
+    "ACCESS-CM2",
+    "ACCESS-ESM1-5",
+    "AWI-ESM-1-1-LR",
+    "BCC-CSM2-MR",
+    "BCC-ESM1",
+    "CAMS-CSM1-0",
+    "CanESM5",
+    "CAS-ESM2-0",
+    "CESM2",
+    "CESM2-WACCM",
+    "CIESM",
+    "CMCC-ESM2",
+    "E3SM-1-0",
+    "E3SM-1-1",
+    "E3SM-2-0",
+    "EC-Earth3",
+    "FGOALS-f3-L",
+    "GFDL-CM4",
+    "GFDL-ESM4",
+    "GISS-E2-1-G",
+    "GISS-E2-1-H",
+    "IPSL-CM5A2-INCA",
+    "IPSL-CM6A-LR",
+    "KIOST-ESM",
+    "MIROC6",
+    "MPI-ESM1-2-LR",
+    "MRI-ESM2-0",
+    "NESM3",
+    "NorESM2-LM",
+    "TaiESM1",
+]
 
 
 # ROOT PATH FOR MODELS CLIMATOLOGIES
-#test_data_path = '$INPUT_DIR$/CMIP5_demo_clims/cmip5.historical.%(model).r1i1p1.mon.pr.198101-200512.AC.v20200426.nc'
-#test_data_path = '/p/css03/cmip5_css02/data/cmip5/output1/CCCma/CanCM4/historical/mon/atmos/Amon/r1i1p1/pr/1/pr_Amon_CanCM4_historical_r1i1p1_196101-200512.nc' 
-test_data_path = '/p/user_pub/pmp/pmp_results/pmp_v1.1.2/diagnostic_results/CMIP_CLIMS/cmip6/historical/v20230823/pr/cmip6.historical.%(model).r1i1p1f1.mon.pr.198101-200512.AC.v20230823.nc'
+test_data_path = "/p/user_pub/pmp/pmp_results/pmp_v1.1.2/diagnostic_results/CMIP_CLIMS/cmip6/historical/v20230823/pr/cmip6.historical.%(model).r1i1p1f1.mon.pr.198101-200512.AC.v20230823.nc"
 
 # ROOT PATH FOR OBSERVATIONS
-#reference_data_path = '$INPUT_DIR$/obs4MIPs_PCMDI_monthly/NOAA-NCEI/GPCP-2-3/mon/pr/gn/v20210727/pr_mon_GPCP-2-3_PCMDI_gn_197901-201907.nc'
-#reference_data_path = '/p/user_pub/PCMDIobs/obs4MIPs/NASA-GSFC/GPCP-Monthly-3-2/mon/pr/gn/v20240408/pr_mon_GPCP-Monthly-3-2_RSS_gn_198301-198312.nc'
-reference_data_path = '/p/user_pub/PCMDIobs/obs4MIPs_legacy/PCMDIobs2_clims/atmos/pr/TRMM-3B43v-7/pr_mon_TRMM-3B43v-7_BE_gn_199801-201712.v20200421.AC.nc'
+reference_data_path = "/p/user_pub/PCMDIobs/obs4MIPs_legacy/PCMDIobs2_clims/atmos/pr/TRMM-3B43v-7/pr_mon_TRMM-3B43v-7_BE_gn_199801-201712.v20200421.AC.nc"
 
 # DIRECTORY WHERE TO PUT RESULTS
-#results_dir = '$OUTPUT_DIR$/monsoon_wang'
-results_dir = '/home/dong12/PMP_240131/pcmdi_metrics/pcmdi_metrics/monsoon_wang/output'
+results_dir = "$OUTPUT_DIR$/monsoon_wang"
 
 # Threshold
 threshold = 2.5 / 86400
+
+# monsoon domain mask based on observations
+obs_mask = True
diff --git a/pcmdi_metrics/monsoon_wang/debug_regions_specs.py b/pcmdi_metrics/monsoon_wang/debug_regions_specs.py
index 6058620fa..1bff30764 100644
--- a/pcmdi_metrics/monsoon_wang/debug_regions_specs.py
+++ b/pcmdi_metrics/monsoon_wang/debug_regions_specs.py
@@ -1,10 +1,13 @@
-import xcdat as xc
+from typing import Union
+
 import xarray as xr
+import xcdat as xc
+
 from pcmdi_metrics.io import load_regions_specs, region_subset
-from typing import Union
-import sys 
-#from pcmdi_metrics.io import da_to_ds, get_longitude, select_subset
-#from .xcdat_dataset_io import da_to_ds
+
+# from pcmdi_metrics.io import da_to_ds, get_longitude, select_subset
+# from .xcdat_dataset_io import da_to_ds
+
 
 def da_to_ds(d: Union[xr.Dataset, xr.DataArray], var: str = "variable") -> xr.Dataset:
     """Convert xarray DataArray to Dataset
@@ -36,8 +39,8 @@ def da_to_ds(d: Union[xr.Dataset, xr.DataArray], var: str = "variable") -> xr.Da
         )
 
 
-#xr.show_versions()
-#sys.exit()
+# xr.show_versions()
+# sys.exit()
 
 
 regions_specs = load_regions_specs()
@@ -46,14 +49,12 @@ def da_to_ds(d: Union[xr.Dataset, xr.DataArray], var: str = "variable") -> xr.Da
 if isinstance(ds, xr.DataArray):
     is_dataArray = True
     print("True")
-    ds = da_to_ds(ds, data_var)
+    ds = da_to_ds(ds)
     print("da_to_ds converted")
 
 
-ds = ds.drop_vars(['time'])
+ds = ds.drop_vars(["time"])
 ds = xc.swap_lon_axis(ds, to=(-180, 180))
 
 
-
-
 mpi_obs_reg = region_subset(ds, "NAFM", data_var="pr", regions_specs=regions_specs)
diff --git a/pcmdi_metrics/monsoon_wang/monsoon_precip_index_fncs.py b/pcmdi_metrics/monsoon_wang/monsoon_precip_index_fncs.py
index 1c84daedf..b36ebdb1f 100644
--- a/pcmdi_metrics/monsoon_wang/monsoon_precip_index_fncs.py
+++ b/pcmdi_metrics/monsoon_wang/monsoon_precip_index_fncs.py
@@ -1,9 +1,7 @@
-import MV2
+from typing import Union
+
 import numpy as np
-import sys
 import xarray as xr
-from scipy.interpolate import griddata
-from typing import Union
 
 #  SEASONAL RANGE - USING ANNUAL CYCLE CLIMATOLGIES 0=Jan, 11=Dec
 
@@ -38,10 +36,7 @@ def da_to_ds(d: Union[xr.Dataset, xr.DataArray], var: str = "variable") -> xr.Da
         )
 
 
-
-
 def regrid(da_in, da_grid, data_var="pr"):
-
     ds_in = da_to_ds(da_in, data_var)
     ds_grid = da_to_ds(da_grid, data_var)
 
@@ -57,21 +52,14 @@ def compute_season(data, season_indices, weights):
     for i in season_indices:
         out += data[i] * weights[i]
         N += weights[i]
-    #out = MV2.array(out)
-    #out.id = data.id
-    #out.setAxisList(data.getAxisList()[1:])
     return out / N
 
 
 def mpd(data):
     """Monsoon precipitation intensity and annual range calculation
-
     .. describe:: Input
-
         *  data
-
             * Assumes climatology array with 12 times step first one January
-
     """
     months_length = [
         31.0,
@@ -90,79 +78,18 @@ def mpd(data):
     mjjas = compute_season(data, [4, 5, 6, 7, 8], months_length)
     ndjfm = compute_season(data, [10, 11, 0, 1, 2], months_length)
     ann = compute_season(data, list(range(12)), months_length)
-    #print("mjjas = ", mjjas)
-    #print("data  = ", data)
-    #print("data.dims = ", data.dims)
-    #print("data.coords = ", data.coords)
 
-    #data_map = data.drop("time")
     data_map = data.isel(time=0)
-    #print("data_map =  ", data_map.dims)
-
-    #sys.exit()
-
-    #cc = xr.DataArray(data.values, coords = data.coords, dims = data.dims)
-    #print("cc =  ", cc.dims)
 
-    #annrange = MV2.subtract(mjjas, ndjfm)
     annrange = mjjas - ndjfm
 
-    #print('annrange = ', annrange)
-
-
-    #lat = annrange.getAxis(0)
-    #lat = annrange['lat']
-
-
-#    i, e = lat.mapInterval((-91, 0, "con"))
-#    print('i = ', i, '   e = ', e)
-#    if i > e:  # reveresedlats
-#        tmp = i + 1
-#        i = e + 1
-#        e = tmp
-
-    #print('annrange = ', annrange)
-    #print('annrange.shape = ', annrange.shape)
-    #print('lat = ', lat)
-    #print('i = ', i, '   e = ', e)
-
-    #annrange[slice(i, e)] = -annrange[slice(i, e)]
+    da_annrange = xr.DataArray(annrange, coords=data_map.coords, dims=data_map.dims)
+    da_annrange = da_annrange.where(da_annrange.lat >= 0, da_annrange * -1)
 
-    #annrange = np.absolute(annrange)
-
-    da_annrange = xr.DataArray(annrange, coords = data_map.coords, dims = data_map.dims)
-    da_annrange = da_annrange.where(da_annrange.lat >= 0, da_annrange*-1)
-
-    #print('slice = ', slice(i, e))
-    #print('annrange = ', annrange)
-    #annrange.id = data.id + "_ar"
-    #annrange.longname = "annual range"
-
-#    sys.exit()
-
-    #mpi = np.divide(annrange, ann, where=ann.astype(bool))
     mpi = np.divide(da_annrange.values, ann, where=ann.astype(bool))
 
-    #mpi = MV2.divide(annrange, ann)
-    #mpi.id = data.id + "_int"
-    #mpi.longname = "intensity"
-
-    #print('mpi.shape = ', mpi.shape)
-    #print('annrange.shape = ', annrange.shape)
-
-    #da_annrange = xr.DataArray(annrange, coords = data_map.coords, dims = data_map.dims)
-    da_mpi = xr.DataArray(mpi, coords = data_map.coords, dims = data_map.dims)
-
-    #print('da_annrange.dims = ', da_annrange.dims)
-    #print('da_mpi_dims = ', da_mpi.dims)
-    #print('da_annrange.coords = ', da_annrange.coords)
-    #print('da_mpi_coords = ', da_mpi.coords)
-
-    #print("da_mpi = ", da_mpi)
-
-    #sys.exit()
+    da_mpi = xr.DataArray(mpi, coords=data_map.coords, dims=data_map.dims)
 
-    #return annrange, mpi
     return da_annrange, da_mpi
 
 
@@ -184,62 +111,43 @@ def mpi_skill_scores(annrange_mod_dom, annrange_obs_dom, threshold=2.5 / 86400.0
 
               * threshold in same units as inputs
     """
-    #print('annrange_mod_dom = ', annrange_mod_dom)
-#    print('annrange_mod_dom.shape = ', annrange_mod_dom.shape)
-    #print('threshold = ', threshold)
     mt = np.ma.greater(annrange_mod_dom, threshold)
     ot = np.ma.greater(annrange_obs_dom, threshold)
 
-    #print('mt = ', mt)
-    #print('type mt = ', type(mt))
-    #print('mt.shape = ', mt.shape)
-
     hitmap = mt * ot  # only where both  mt and ot are True
     hit = float(hitmap.sum())
-    #print("hitmap = ", hitmap)
-    #print("hit = ", hit)
 
     xor = np.ma.logical_xor(mt, ot)
     missmap = xor * ot
-    #missed = float(MV2.sum(missmap))
-    #print("xor = ", xor)
-    #print("missmap = ", missmap)
-    #print("missed = ", missed)
     missed = float(missmap.sum())
-    #print("missed = ", missed)
 
     falarmmap = xor * mt
-    #falarm = float(MV2.sum(falarmmap))
-    #print("falarm = ", falarm)
     falarm = float(falarmmap.sum())
-    #print("falarm = ", falarm)
 
     if (hit + missed + falarm) > 0.0:
         score = hit / (hit + missed + falarm)
     else:
         score = 1.0e20
 
-    #hitmap.id = "hit"
-    #missmap.id = "miss"
-    #falarmmap.id = "false_alarm"
-
-    #for a in [hitmap, missmap, falarmmap]:
-    #    a.setAxisList(annrange_mod_dom.getAxisList())
-
-
-    xr_hitmap = xr.DataArray(hitmap, name='hitmap',\
-        coords={ 'lat': annrange_mod_dom.lat, 'lon': annrange_mod_dom.lon},\
-        dims=[ 'lat', 'lon'])
-
-    xr_missmap = xr.DataArray(missmap, name='missmap',\
-        coords={ 'lat': annrange_mod_dom.lat, 'lon': annrange_mod_dom.lon},\
-        dims=[ 'lat', 'lon'])
-
-    xr_falarmmap = xr.DataArray(falarmmap, name='falarmmap',\
-        coords={ 'lat': annrange_mod_dom.lat, 'lon': annrange_mod_dom.lon},\
-        dims=[ 'lat', 'lon'])
-
-
+    xr_hitmap = xr.DataArray(
+        hitmap,
+        name="hitmap",
+        coords={"lat": annrange_mod_dom.lat, "lon": annrange_mod_dom.lon},
+        dims=["lat", "lon"],
+    )
+
+    xr_missmap = xr.DataArray(
+        missmap,
+        name="missmap",
+        coords={"lat": annrange_mod_dom.lat, "lon": annrange_mod_dom.lon},
+        dims=["lat", "lon"],
+    )
+
+    xr_falarmmap = xr.DataArray(
+        falarmmap,
+        name="falarmmap",
+        coords={"lat": annrange_mod_dom.lat, "lon": annrange_mod_dom.lon},
+        dims=["lat", "lon"],
+    )
 
-    #return hit, missed, falarm, score, hitmap, missmap, falarmmap
     return hit, missed, falarm, score, xr_hitmap, xr_missmap, xr_falarmmap
diff --git a/pcmdi_metrics/monsoon_wang/monsoon_wang_driver.py b/pcmdi_metrics/monsoon_wang/monsoon_wang_driver.py
index dc9503c95..623764a1a 100644
--- a/pcmdi_metrics/monsoon_wang/monsoon_wang_driver.py
+++ b/pcmdi_metrics/monsoon_wang/monsoon_wang_driver.py
@@ -2,23 +2,19 @@
 
 import collections
 import os
+import sys
 
 import cdms2
 import numpy
 import numpy as np
-from genutil import statistics
+import xarray as xr
+from monsoon_precip_index_fncs import mpd, mpi_skill_scores, regrid
 
 import pcmdi_metrics
 from pcmdi_metrics import resources
+from pcmdi_metrics.io import load_regions_specs, region_subset
 from pcmdi_metrics.mean_climate.lib.pmp_parser import PMPParser
-from monsoon_precip_index_fncs import mpd, mpi_skill_scores
-from monsoon_precip_index_fncs import regrid, da_to_ds
 from pcmdi_metrics.utils import StringConstructor
-import sys
-import xcdat as xc
-import xarray as xr
-from pcmdi_metrics.io import load_regions_specs, region_subset
-#import da_to_ds
 
 
 def create_monsoon_wang_parser():
@@ -28,6 +24,7 @@ def create_monsoon_wang_parser():
     P.use("--results_dir")
     P.use("--reference_data_path")
     P.use("--test_data_path")
+    P.use("--obs_mask")
 
     P.add_argument(
         "--outnj",
@@ -110,17 +107,10 @@ def monsoon_wang_runner(args):
     # PMP monthly default PR obs
     cdms2.axis.longitude_aliases.append("longitude_prclim_mpd")
     cdms2.axis.latitude_aliases.append("latitude_prclim_mpd")
-    #fobs = cdms2.open(args.reference_data_path)
-    #dobs_orig = fobs(args.obsvar)
-    #fobs.close()
-    print("args.reference_data_path  =  ", args.reference_data_path)
     fobs = xr.open_dataset(args.reference_data_path, decode_times=False)
     dobs_orig = fobs[args.obsvar]
     fobs.close()
 
-    #obsgrid = dobs_orig.getGrid()
-    #print(" obsgrid =  ", obsgrid)
-
     ########################################
 
     # FCN TO COMPUTE GLOBAL ANNUAL RANGE AND MONSOON PRECIP INDEX
@@ -128,14 +118,11 @@ def monsoon_wang_runner(args):
     annrange_obs, mpi_obs = mpd(dobs_orig)
 
     # create monsoon domain mask based on observations: annual range > 2.5 mm/day
-
-    domain_mask_obs = xr.where(annrange_obs>thr, 1, 0)
-    domain_mask_obs.name = 'mask'
-    #domain_mask_obs.to_netcdf('/home/dong12/PMP_240131/pcmdi_metrics/pcmdi_metrics/monsoon_wang/domain_mask_obs.nc')
-
-    mpi_obs = mpi_obs.where(domain_mask_obs)
-    #sys.exit()
-
+    if args.obs_mask:
+        domain_mask_obs = xr.where(annrange_obs > thr, 1, 0)
+        domain_mask_obs.name = "mask"
+        mpi_obs = mpi_obs.where(domain_mask_obs)
+        nout_mpi_obs = os.path.join(outpathdata, "mpi_obs_masked.nc")
 
     #########################################
     # SETUP WHERE TO OUTPUT RESULTING DATA (netcdf)
@@ -171,12 +158,8 @@ def monsoon_wang_runner(args):
             gmods.append(mod)
 
     if len(gmods) == 0:
-        print("gmods = ", gmods)
         raise RuntimeError("No model file found!")
 
-    print("nout = ", nout)
-    #print("jout = ", jout)
-    print("gmods = ", gmods)
     #########################################
 
     egg_pth = resources.resource_path()
@@ -192,27 +175,8 @@ def monsoon_wang_runner(args):
         locals,
     )
 
-    #print("locals = ", locals)
-    #print("globals = ", globals)
-    print('os.path.join(egg_pth, "default_regions.py")  = ', os.path.join(egg_pth, "default_regions.py"))
-    # /home/dong12/miniconda3/envs/PMP_240423/share/pmp/default_regions.py
-
     regions_specs = locals["regions_specs"]
     doms = ["AllMW", "AllM", "NAMM", "SAMM", "NAFM", "SAFM", "ASM", "AUSM"]
-    #doms = ["AllMW"]#, "AllM", "NAMM", "SAMM", "NAFM", "SAFM", "ASM", "AUSM"]
-    #doms = ["AUSM"]
-    #doms = ["NAMM"]
-    #doms = ["NAFM"]
-    
-    #print('region_specs = ',regions_specs)
-    #print('region_specs["AllMW"] = ',regions_specs['AllMW'])
-    #print('region_specs["AllMW"]["domain"] = ',regions_specs['AllMW']['domain'])
-    #print('region_specs["AllM"]["domain"] = ',regions_specs['AllM']['domain'])
-    #print('region_specs["ASM"]["domain"] = ',regions_specs['ASM']['domain'])
-
-
-    #sys.exit()
-
 
     mpi_stats_dic = {}
     for i, mod in enumerate(gmods):
@@ -230,120 +194,50 @@ def monsoon_wang_runner(args):
 
         mpi_stats_dic[mod] = {}
 
-        print(
-            "******************************************************************************************"
-        )
         print("modelFile =  ", modelFile)
-        #f = cdms2.open(modelFile)
-        #d_orig = f(var)
         f = xr.open_dataset(modelFile)
         d_orig = f[var]
 
-
         annrange_mod, mpi_mod = mpd(d_orig)
 
-        domain_mask_mod = xr.where(annrange_mod>thr, 1, 0)
-        #domain_mask_mod.name = 'mask'
+        domain_mask_mod = xr.where(annrange_mod > thr, 1, 0)
         mpi_mod = mpi_mod.where(domain_mask_mod)
 
-        #print('mod_annrange.dims = ', annrange_mod.dims)
-        #print('obs_annrange_dims = ', annrange_obs.dims)
-        #print('mod_annrange.coords = ', annrange_mod.coords)
-        #print('obs_annrange_coords = ', annrange_obs.coords)
-
-
-#        sys.exit()
-
-#        annrange_mod = annrange_mod.regrid(
-#            obsgrid, regridTool="regrid2", regridMethod="conserve", mkCyclic=True
-#        )
-#        mpi_mod = mpi_mod.regrid(
-#            obsgrid, regridTool="regrid2", regridMethod="conserve", mkCyclic=True
-#        )
-
-        #print('annrange_obs =  ', annrange_obs)
         lats = annrange_obs.lat[0]
         latn = annrange_obs.lat[-1]
         lone = annrange_obs.lon[-1]
         lonw = annrange_obs.lon[0]
 
-        #annrange_obs.to_netcdf("annrange_obs.nc")
-
-        #annrange_obs = annrange_obs.interp(lat=mpi_mod.lat, lon=mpi_mod.lon, method='linear')
         annrange_obs = regrid(annrange_obs, annrange_mod)
-        
-        #annrange_obs.to_netcdf("annrange_obs_interp.nc")
-
-        #print('annrange_obs =  ', annrange_obs)
-        #print('annrange_obs =  ', annrange_obs.sel(lat=slice(lats, latn), lon=slice(lonw, lone)))
 
-        #mpi_obs = mpi_obs.interp(lat=mpi_mod.lat, lon=mpi_mod.lon, method='linear')
         mpi_obs = regrid(mpi_obs, mpi_mod)
 
-
-        #print("mpi_obs = ", mpi_obs)
-
-
         regions_specs = load_regions_specs()
-        #print("regions_specs - ", regions_specs)
-#        print("list(regions_specs.keys())",  list(regions_specs.keys()))
 
         for dom in doms:
             mpi_stats_dic[mod][dom] = {}
 
-            #reg_sel = regions_specs[dom]["domain"]
-
             print("dom =  ", dom)
 
-            #mpi_obs.to_netcdf("xd_mpi_obs.nc")
-
             mpi_obs_reg = region_subset(mpi_obs, dom)
-            #mpi_obs_reg = region_subset(mpi_obs, dom, data_var="pr", regions_specs=regions_specs)
-            #mpi_obs_reg = region_subset(mpi_obs, dom)#, var="pr", regions_specs=regions_specs)
-            #mpi_obs_reg = region_subset(mpi_obs, dom, var="pr")#, regions_specs=regions_specs)
-            #mpi_obs_reg = region_subset(mpi_obs, regions_specs, region=dom)#, var="pr", regions_specs=regions_specs)
-            #mpi_obs_reg = region_subset(mpi_obs, regions_specs=regions_specs, region=dom)
-
-            #print("mpi_obs_reg =  ", mpi_obs_reg)
-            #sys.exit()
-
-            #mpi_obs_reg = mpi_obs(reg_sel)
-
-            #mpi_obs_reg_sd = float(statistics.std(mpi_obs_reg, axis="xy"))
-            mpi_obs_reg_sd = mpi_obs_reg.std(dim=['lat', 'lon'])
-
-            #mpi_mod_reg = mpi_mod(reg_sel)
-            #mpi_mod_reg = region_subset(mpi_mod, regions_specs=regions_specs, region=dom)
+            mpi_obs_reg_sd = mpi_obs_reg.std(dim=["lat", "lon"])
             mpi_mod_reg = region_subset(mpi_mod, dom)
 
-            #cor = float(statistics.correlation(mpi_mod_reg, mpi_obs_reg, axis="xy"))
             da1_flat = mpi_mod_reg.values.ravel()
             da2_flat = mpi_obs_reg.values.ravel()
-            #print("da1_flat = ", da1_flat)
-            #print("da2_flat = ", da2_flat)
-
-            #cor = np.corrcoef(da1_flat, da2_flat)[0, 1]
-            cor = np.ma.corrcoef(np.ma.masked_invalid(da1_flat), np.ma.masked_invalid(da2_flat) )[0, 1]
 
-            print("cor = ", cor)
+            cor = np.ma.corrcoef(
+                np.ma.masked_invalid(da1_flat), np.ma.masked_invalid(da2_flat)
+            )[0, 1]
 
-            sys.exit()
-
-            #rms = float(statistics.rms(mpi_mod_reg, mpi_obs_reg, axis="xy"))
             squared_diff = (mpi_mod_reg - mpi_obs_reg) ** 2
             mean_squared_error = squared_diff.mean(skipna=True)
             rms = np.sqrt(mean_squared_error)
 
             rmsn = rms / mpi_obs_reg_sd
-            #del(mpi_mod_reg)
-            #del(mpi_obs_reg)
 
             #  DOMAIN SELECTED FROM GLOBAL ANNUAL RANGE FOR MODS AND OBS
-            #annrange_mod_dom = annrange_mod(reg_sel)
-            #annrange_obs_dom = annrange_obs(reg_sel)
 
-            #annrange_mod_dom = region_subset(annrange_mod, regions_specs=regions_specs, region=dom)
-            #annrange_obs_dom = region_subset(annrange_obs, regions_specs=regions_specs, region=dom)
             annrange_mod_dom = region_subset(annrange_mod, dom)
             annrange_obs_dom = region_subset(annrange_obs, dom)
 
@@ -361,26 +255,22 @@ def monsoon_wang_runner(args):
 
             # SAVE ANNRANGE AND HIT MISS AND FALSE ALARM FOR EACH MOD DOM
             fm = os.path.join(nout, "_".join([mod, dom, "wang-monsoon_xcdat.nc"]))
-            #g = cdms2.open(fm, "w")
-            #g.write(annrange_mod_dom)
-            #g.write(hitmap, dtype=numpy.int32)
-            #g.write(missmap, dtype=numpy.int32)
-            #g.write(falarmmap, dtype=numpy.int32)
-            #g.close()
-            ds_out = xr.Dataset({
+            ds_out = xr.Dataset(
+                {
                     "obsmap": annrange_obs_dom,
                     "modmap": annrange_mod_dom,
                     "hitmap": hitmap,
                     "missmap": missmap,
-                    "falarmmap": falarmmap
-                    })
+                    "falarmmap": falarmmap,
+                }
+            )
             ds_out.to_netcdf(fm)
         f.close()
 
         if np.isnan(cor):
+            print("invalid correlation values")
             sys.exit()
 
-
     #  OUTPUT METRICS TO JSON FILE
     OUT = pcmdi_metrics.io.base.Base(os.path.abspath(jout), json_filename)
 
diff --git a/pcmdi_metrics/monsoon_wang/output/monsoon_wang.json b/pcmdi_metrics/monsoon_wang/output/monsoon_wang.json
index 7c91674e6..b3eac7e26 100644
--- a/pcmdi_metrics/monsoon_wang/output/monsoon_wang.json
+++ b/pcmdi_metrics/monsoon_wang/output/monsoon_wang.json
@@ -1401,4 +1401,4 @@
         "domain",
         "statistic"
     ]
-}
\ No newline at end of file
+}
diff --git a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cdms.json b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cdms.json
index d5a29e28e..491c49035 100644
--- a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cdms.json
+++ b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cdms.json
@@ -99,4 +99,4 @@
         "domain",
         "statistic"
     ]
-}
\ No newline at end of file
+}
diff --git a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip5_xcdat_240906.json b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip5_xcdat_240906.json
index 6490729d5..1c232646c 100644
--- a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip5_xcdat_240906.json
+++ b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip5_xcdat_240906.json
@@ -1401,4 +1401,4 @@
         "domain",
         "statistic"
     ]
-}
\ No newline at end of file
+}
diff --git a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip6_xcdat_240906.json b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip6_xcdat_240906.json
index 8c101b784..2677ade21 100644
--- a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip6_xcdat_240906.json
+++ b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_cmip6_xcdat_240906.json
@@ -1317,4 +1317,4 @@
         "domain",
         "statistic"
     ]
-}
\ No newline at end of file
+}
diff --git a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_multi_model.json b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_multi_model.json
index 6440ca908..8f4b69da2 100644
--- a/pcmdi_metrics/monsoon_wang/output/monsoon_wang_multi_model.json
+++ b/pcmdi_metrics/monsoon_wang/output/monsoon_wang_multi_model.json
@@ -1317,4 +1317,4 @@
         "domain",
         "statistic"
     ]
-}
\ No newline at end of file
+}
diff --git a/pcmdi_metrics/monsoon_wang/test_param.py b/pcmdi_metrics/monsoon_wang/test_param.py
index 132ca4051..379ff2c4e 100644
--- a/pcmdi_metrics/monsoon_wang/test_param.py
+++ b/pcmdi_metrics/monsoon_wang/test_param.py
@@ -1,54 +1,164 @@
-import os
-
 #
 #  OPTIONS ARE SET BY USER IN THIS FILE AS INDICATED BELOW BY:
 #
 #
 
 # LIST OF MODEL VERSIONS TO BE TESTED
-modnames = ['CanCM4']
-modnames = ['CESM2']
-#modnames = ['IPSL-CM6A-LR']
-#modnames = ['ACCESS-CM2', 'ACCESS-ESM1-5', 'AWI-CM-1-1-MR', 'AWI-ESM-1-1-LR', 'BCC-CSM2-MR', 'BCC-ESM1', 'CAMS-CSM1-0', 'CanESM5-1', 'CanESM5', 'CAS-ESM2-0', 'CESM2-FV2', 'CESM2', 'CESM2-WACCM-FV2', 'CESM2-WACCM', 'CIESM', 'CMCC-CM2-HR4', 'CMCC-CM2-SR5', 'CMCC-ESM2', 'E3SM-1-0', 'E3SM-1-1-ECA', 'E3SM-1-1', 'E3SM-2-0', 'EC-Earth3-AerChem', 'EC-Earth3-CC', 'EC-Earth3', 'EC-Earth3-Veg-LR', 'EC-Earth3-Veg', 'FGOALS-f3-L', 'FIO-ESM-2-0', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G-CC', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'GISS-E2-2-H', 'ICON-ESM-LR', 'INM-CM4-8', 'INM-CM5-0', 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR-INCA', 'IPSL-CM6A-LR', 'KACE-1-0-G', 'KIOST-ESM', 'MCM-UA-1-0', 'MIROC6', 'MPI-ESM-1-2-HAM', 'MPI-ESM1-2-HR', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorCPM1', 'NorESM2-LM', 'NorESM2-MM', 'SAM0-UNICON', 'TaiESM1']
+modnames = ["CanCM4"]
+modnames = ["CESM2"]
+# modnames = ['IPSL-CM6A-LR']
+# modnames = ['ACCESS-CM2', 'ACCESS-ESM1-5', 'AWI-CM-1-1-MR', 'AWI-ESM-1-1-LR', 'BCC-CSM2-MR', 'BCC-ESM1', 'CAMS-CSM1-0', 'CanESM5-1', 'CanESM5', 'CAS-ESM2-0', 'CESM2-FV2', 'CESM2', 'CESM2-WACCM-FV2', 'CESM2-WACCM', 'CIESM', 'CMCC-CM2-HR4', 'CMCC-CM2-SR5', 'CMCC-ESM2', 'E3SM-1-0', 'E3SM-1-1-ECA', 'E3SM-1-1', 'E3SM-2-0', 'EC-Earth3-AerChem', 'EC-Earth3-CC', 'EC-Earth3', 'EC-Earth3-Veg-LR', 'EC-Earth3-Veg', 'FGOALS-f3-L', 'FIO-ESM-2-0', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G-CC', 'GISS-E2-1-G', 'GISS-E2-1-H', 'GISS-E2-2-G', 'GISS-E2-2-H', 'ICON-ESM-LR', 'INM-CM4-8', 'INM-CM5-0', 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR-INCA', 'IPSL-CM6A-LR', 'KACE-1-0-G', 'KIOST-ESM', 'MCM-UA-1-0', 'MIROC6', 'MPI-ESM-1-2-HAM', 'MPI-ESM1-2-HR', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorCPM1', 'NorESM2-LM', 'NorESM2-MM', 'SAM0-UNICON', 'TaiESM1']
 
 #  trouble model ICON-ESM-LR'
 
 
-#cmip6 models
-modnames = ['ACCESS-CM2', 'ACCESS-ESM1-5', 'AWI-ESM-1-1-LR', 'BCC-CSM2-MR', 'BCC-ESM1', 'CAMS-CSM1-0', 'CanESM5', 'CAS-ESM2-0', 'CESM2', 'CESM2-WACCM', 'CIESM', 'CMCC-ESM2', 'E3SM-1-0', 'E3SM-1-1', 'E3SM-2-0', 'EC-Earth3', 'FGOALS-f3-L', 'GFDL-CM4', 'GFDL-ESM4', 'GISS-E2-1-G', 'GISS-E2-1-H', 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR', 'KIOST-ESM', 'MIROC6', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorESM2-LM', 'TaiESM1']
-
-
-#modnames = [ 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR', 'KIOST-ESM', 'MIROC6', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorESM2-LM', 'TaiESM1']
-
-#cmip5 models
-#modnames = [ACCESS1-0, BCC-CSM1-1-M, BNU-ESM, CanCM4, CanESM2, CCSM4, CESM1-BGC, CESM1-CAM5, CESM1-FASTCHEM, CESM1-WACCM, CMCC-CESM, CMCC-CM, CMCC-CMS, CNRM-CM5-2, CNRM-CM5, CSIRO-Mk3-6-0, FGOALS-g2, FIO-ESM, GFDL-CM2p1, GFDL-CM3, GFDL-ESM2G, GFDL-ESM2M, GISS-E2-H-CC, GISS-E2-H, GISS-E2-R-CC, GISS-E2-R, HadCM3, HadGEM2-AO, INMCM4, IPSL-CM5A-LR, IPSL-CM5A-MR, IPSL-CM5B-LR, MIROC4h, MIROC5, MIROC-ESM-CHEM, MIROC-ESM, MPI-ESM-LR, MPI-ESM-MR, MPI-ESM-P, MRI-CGCM3, MRI-ESM1, NorESM1-ME, NorESM1-M]
-
-modnames = ['ACCESS1-0', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-CAM5', 'CESM1-WACCM', 'CMCC-CESM', 'CNRM-CM5', 'CSIRO-Mk3-6-0', 'FGOALS-g2', 'FIO-ESM', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadCM3', 'HadGEM2-AO', 'INMCM4', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC4h', 'MIROC5', 'MIROC-ESM', 'MPI-ESM-LR', 'MPI-ESM-MR', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M']
-
-modnames = ['CMCC-CESM', "CNRM-CM5"]
-modnames = ['ACCESS1-0', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-CAM5', 'CESM1-WACCM', 'CMCC-CESM', 'CNRM-CM5', 'CSIRO-Mk3-6-0']
-modnames = ['CESM1-WACCM', 'CMCC-CESM', 'CNRM-CM5', 'CSIRO-Mk3-6-0']
-modnames = ['CNRM-CM5', 'CSIRO-Mk3-6-0', 'FGOALS-g2', 'FIO-ESM', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadCM3', 'HadGEM2-AO', 'INMCM4', 'IPSL-CM5A-LR', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC4h', 'MIROC5', 'MIROC-ESM', 'MPI-ESM-LR', 'MPI-ESM-MR', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M','ACCESS1-0', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-CAM5', 'CESM1-WACCM', 'CMCC-CESM']
-#modnames = ['IPSL-CM5A-LR','CNRM-CM5', 'CSIRO-Mk3-6-0', 'FGOALS-g2', 'FIO-ESM', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadCM3', 'HadGEM2-AO', 'INMCM4', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC4h', 'MIROC5', 'MIROC-ESM', 'MPI-ESM-LR', 'MPI-ESM-MR', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M','ACCESS1-0', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-CAM5', 'CESM1-WACCM', 'CMCC-CESM']
+# cmip6 models
+modnames = [
+    "ACCESS-CM2",
+    "ACCESS-ESM1-5",
+    "AWI-ESM-1-1-LR",
+    "BCC-CSM2-MR",
+    "BCC-ESM1",
+    "CAMS-CSM1-0",
+    "CanESM5",
+    "CAS-ESM2-0",
+    "CESM2",
+    "CESM2-WACCM",
+    "CIESM",
+    "CMCC-ESM2",
+    "E3SM-1-0",
+    "E3SM-1-1",
+    "E3SM-2-0",
+    "EC-Earth3",
+    "FGOALS-f3-L",
+    "GFDL-CM4",
+    "GFDL-ESM4",
+    "GISS-E2-1-G",
+    "GISS-E2-1-H",
+    "IPSL-CM5A2-INCA",
+    "IPSL-CM6A-LR",
+    "KIOST-ESM",
+    "MIROC6",
+    "MPI-ESM1-2-LR",
+    "MRI-ESM2-0",
+    "NESM3",
+    "NorESM2-LM",
+    "TaiESM1",
+]
+
+
+# modnames = [ 'IPSL-CM5A2-INCA', 'IPSL-CM6A-LR', 'KIOST-ESM', 'MIROC6', 'MPI-ESM1-2-LR', 'MRI-ESM2-0', 'NESM3', 'NorESM2-LM', 'TaiESM1']
+
+# cmip5 models
+# modnames = [ACCESS1-0, BCC-CSM1-1-M, BNU-ESM, CanCM4, CanESM2, CCSM4, CESM1-BGC, CESM1-CAM5, CESM1-FASTCHEM, CESM1-WACCM, CMCC-CESM, CMCC-CM, CMCC-CMS, CNRM-CM5-2, CNRM-CM5, CSIRO-Mk3-6-0, FGOALS-g2, FIO-ESM, GFDL-CM2p1, GFDL-CM3, GFDL-ESM2G, GFDL-ESM2M, GISS-E2-H-CC, GISS-E2-H, GISS-E2-R-CC, GISS-E2-R, HadCM3, HadGEM2-AO, INMCM4, IPSL-CM5A-LR, IPSL-CM5A-MR, IPSL-CM5B-LR, MIROC4h, MIROC5, MIROC-ESM-CHEM, MIROC-ESM, MPI-ESM-LR, MPI-ESM-MR, MPI-ESM-P, MRI-CGCM3, MRI-ESM1, NorESM1-ME, NorESM1-M]
+
+modnames = [
+    "ACCESS1-0",
+    "BCC-CSM1-1-M",
+    "BNU-ESM",
+    "CanCM4",
+    "CanESM2",
+    "CCSM4",
+    "CESM1-CAM5",
+    "CESM1-WACCM",
+    "CMCC-CESM",
+    "CNRM-CM5",
+    "CSIRO-Mk3-6-0",
+    "FGOALS-g2",
+    "FIO-ESM",
+    "GFDL-CM3",
+    "GFDL-ESM2G",
+    "GFDL-ESM2M",
+    "GISS-E2-H",
+    "GISS-E2-R",
+    "HadCM3",
+    "HadGEM2-AO",
+    "INMCM4",
+    "IPSL-CM5A-LR",
+    "IPSL-CM5A-MR",
+    "IPSL-CM5B-LR",
+    "MIROC4h",
+    "MIROC5",
+    "MIROC-ESM",
+    "MPI-ESM-LR",
+    "MPI-ESM-MR",
+    "MRI-CGCM3",
+    "MRI-ESM1",
+    "NorESM1-M",
+]
+
+modnames = ["CMCC-CESM", "CNRM-CM5"]
+modnames = [
+    "ACCESS1-0",
+    "BCC-CSM1-1-M",
+    "BNU-ESM",
+    "CanCM4",
+    "CanESM2",
+    "CCSM4",
+    "CESM1-CAM5",
+    "CESM1-WACCM",
+    "CMCC-CESM",
+    "CNRM-CM5",
+    "CSIRO-Mk3-6-0",
+]
+modnames = ["CESM1-WACCM", "CMCC-CESM", "CNRM-CM5", "CSIRO-Mk3-6-0"]
+modnames = [
+    "CNRM-CM5",
+    "CSIRO-Mk3-6-0",
+    "FGOALS-g2",
+    "FIO-ESM",
+    "GFDL-CM3",
+    "GFDL-ESM2G",
+    "GFDL-ESM2M",
+    "GISS-E2-H",
+    "GISS-E2-R",
+    "HadCM3",
+    "HadGEM2-AO",
+    "INMCM4",
+    "IPSL-CM5A-LR",
+    "IPSL-CM5A-MR",
+    "IPSL-CM5B-LR",
+    "MIROC4h",
+    "MIROC5",
+    "MIROC-ESM",
+    "MPI-ESM-LR",
+    "MPI-ESM-MR",
+    "MRI-CGCM3",
+    "MRI-ESM1",
+    "NorESM1-M",
+    "ACCESS1-0",
+    "BCC-CSM1-1-M",
+    "BNU-ESM",
+    "CanCM4",
+    "CanESM2",
+    "CCSM4",
+    "CESM1-CAM5",
+    "CESM1-WACCM",
+    "CMCC-CESM",
+]
+# modnames = ['IPSL-CM5A-LR','CNRM-CM5', 'CSIRO-Mk3-6-0', 'FGOALS-g2', 'FIO-ESM', 'GFDL-CM3', 'GFDL-ESM2G', 'GFDL-ESM2M', 'GISS-E2-H', 'GISS-E2-R', 'HadCM3', 'HadGEM2-AO', 'INMCM4', 'IPSL-CM5A-MR', 'IPSL-CM5B-LR', 'MIROC4h', 'MIROC5', 'MIROC-ESM', 'MPI-ESM-LR', 'MPI-ESM-MR', 'MRI-CGCM3', 'MRI-ESM1', 'NorESM1-M','ACCESS1-0', 'BCC-CSM1-1-M', 'BNU-ESM', 'CanCM4', 'CanESM2', 'CCSM4', 'CESM1-CAM5', 'CESM1-WACCM', 'CMCC-CESM']
 
 # ROOT PATH FOR MODELS CLIMATOLOGIES
-#test_data_path = '$INPUT_DIR$/CMIP5_demo_clims/cmip5.historical.%(model).r1i1p1.mon.pr.198101-200512.AC.v20200426.nc'
-#test_data_path = '/p/css03/cmip5_css02/data/cmip5/output1/CCCma/CanCM4/historical/mon/atmos/Amon/r1i1p1/pr/1/pr_Amon_CanCM4_historical_r1i1p1_196101-200512.nc' 
+# test_data_path = '$INPUT_DIR$/CMIP5_demo_clims/cmip5.historical.%(model).r1i1p1.mon.pr.198101-200512.AC.v20200426.nc'
+# test_data_path = '/p/css03/cmip5_css02/data/cmip5/output1/CCCma/CanCM4/historical/mon/atmos/Amon/r1i1p1/pr/1/pr_Amon_CanCM4_historical_r1i1p1_196101-200512.nc'
 
 # cmip6
-#test_data_path = '/p/user_pub/pmp/pmp_results/pmp_v1.1.2/diagnostic_results/CMIP_CLIMS/cmip6/historical/v20230823/pr/cmip6.historical.%(model).r1i1p1f1.mon.pr.198101-200512.AC.v20230823.nc'
+# test_data_path = '/p/user_pub/pmp/pmp_results/pmp_v1.1.2/diagnostic_results/CMIP_CLIMS/cmip6/historical/v20230823/pr/cmip6.historical.%(model).r1i1p1f1.mon.pr.198101-200512.AC.v20230823.nc'
 
 # cmip5
-test_data_path = '/p/user_pub/pmp/pmp_results/pmp_v1.1.2/diagnostic_results/CMIP_CLIMS/cmip5/historical/v20230323/pr/cmip5.historical.%(model).r1i1p1.mon.pr.198101-200512.AC.v20230323.nc'
+test_data_path = "/p/user_pub/pmp/pmp_results/pmp_v1.1.2/diagnostic_results/CMIP_CLIMS/cmip5/historical/v20230323/pr/cmip5.historical.%(model).r1i1p1.mon.pr.198101-200512.AC.v20230323.nc"
 
 # ROOT PATH FOR OBSERVATIONS
-#reference_data_path = '$INPUT_DIR$/obs4MIPs_PCMDI_monthly/NOAA-NCEI/GPCP-2-3/mon/pr/gn/v20210727/pr_mon_GPCP-2-3_PCMDI_gn_197901-201907.nc'
-#reference_data_path = '/p/user_pub/PCMDIobs/obs4MIPs/NASA-GSFC/GPCP-Monthly-3-2/mon/pr/gn/v20240408/pr_mon_GPCP-Monthly-3-2_RSS_gn_198301-198312.nc'
-reference_data_path = '/p/user_pub/PCMDIobs/obs4MIPs_legacy/PCMDIobs2_clims/atmos/pr/TRMM-3B43v-7/pr_mon_TRMM-3B43v-7_BE_gn_199801-201712.v20200421.AC.nc'
+# reference_data_path = '$INPUT_DIR$/obs4MIPs_PCMDI_monthly/NOAA-NCEI/GPCP-2-3/mon/pr/gn/v20210727/pr_mon_GPCP-2-3_PCMDI_gn_197901-201907.nc'
+# reference_data_path = '/p/user_pub/PCMDIobs/obs4MIPs/NASA-GSFC/GPCP-Monthly-3-2/mon/pr/gn/v20240408/pr_mon_GPCP-Monthly-3-2_RSS_gn_198301-198312.nc'
+reference_data_path = "/p/user_pub/PCMDIobs/obs4MIPs_legacy/PCMDIobs2_clims/atmos/pr/TRMM-3B43v-7/pr_mon_TRMM-3B43v-7_BE_gn_199801-201712.v20200421.AC.nc"
 
 # DIRECTORY WHERE TO PUT RESULTS
-#results_dir = '$OUTPUT_DIR$/monsoon_wang'
-results_dir = '/home/dong12/PMP_240131/pcmdi_metrics/pcmdi_metrics/monsoon_wang/output'
+# results_dir = '$OUTPUT_DIR$/monsoon_wang'
+results_dir = "/home/dong12/PMP_240131/pcmdi_metrics/pcmdi_metrics/monsoon_wang/output"
 
 # Threshold
 threshold = 2.5 / 86400
+
+# monsoon domain mask based on observations
+obs_mask = True
diff --git a/share/DefArgsCIA.json b/share/DefArgsCIA.json
index cd33a055d..8507f33ba 100644
--- a/share/DefArgsCIA.json
+++ b/share/DefArgsCIA.json
@@ -163,4 +163,4 @@
         ],
         "help":"A list of variables to be processed"
     }
-}
\ No newline at end of file
+}
diff --git a/share/default_regions.py b/share/default_regions.py
index 4f21606a1..72eb7804a 100755
--- a/share/default_regions.py
+++ b/share/default_regions.py
@@ -50,7 +50,7 @@
     },
     # North African Monsoon
     "NAFM": {
-        #"domain": cdutil.region.domain(latitude=(0.0, 45.0), longitude=(310.0, 60.0))
+        # "domain": cdutil.region.domain(latitude=(0.0, 45.0), longitude=(310.0, 60.0))
         "domain": cdutil.region.domain(latitude=(0.0, 45.0), longitude=(-50, 60.0))
     },
     # South African Monsoon
diff --git a/test_landmask_generation.py b/test_landmask_generation.py
index 500f5962d..b4b84e37e 100644
--- a/test_landmask_generation.py
+++ b/test_landmask_generation.py
@@ -1,7 +1,9 @@
-from pcmdi_metrics.utils import create_land_sea_mask
+import os
+
 import xarray as xr
 import xcdat as xc
-import os
+
+from pcmdi_metrics.utils import create_land_sea_mask
 
 demo_dir = "/Users/lee1043/Documents/Research/git/pcmdi_metrics_20230620_pcmdi/pcmdi_metrics/doc/jupyter/Demo/demo_data/"
 demo_file = "CMIP5_demo_data/ts_Amon_ACCESS1-0_historical_r1i1p1_185001-200512.nc"
@@ -15,9 +17,10 @@
 mask2.plot()
 (mask2 - mask).plot()
 
-from pcmdi_metrics import resources
 import os
 
+from pcmdi_metrics import resources
+
 egg_pth = resources.resource_path()
 source_path = os.path.join(egg_pth, "navy_land.nc")
 
@@ -25,11 +28,12 @@
 
 ds_navy["sftlf"].plot()
 
-import cdms2, cdutil
+import cdms2
+import cdutil
 
 f = cdms2.open(dummy_input)
 
-dc = f('ts')
+dc = f("ts")
 
 mask_cdms = cdutil.generateLandSeaMask(dc)
 
@@ -38,10 +42,8 @@
 
 
 mask_cdms_da = xr.DataArray(
-    mask_cdms, 
-    coords=[list(lat), list(lon)], 
-    dims=["lat", "lon"], 
-    name=dc.id)
+    mask_cdms, coords=[list(lat), list(lon)], dims=["lat", "lon"], name=dc.id
+)
 
 
 mask_cdms_da.plot()
@@ -53,10 +55,9 @@
 (mask2 - mask_cdms_da).sum()
 
 #
-#%%time
-#a = mask.to_numpy()
+# %%time
+# a = mask.to_numpy()
 #
 #
-#%%time
-#b = mask.data
-
+# %%time
+# b = mask.data