Fix PSM for Mg/Ca

cfr/da/enkf.py

@@ -240,6 +240,17 @@ def run(self, recon_yrs=np.arange(1, 2001), recon_loc_rad=25000, recon_timescale
                     recon_timescale=recon_timescale,
                     debug=debug,
                 )
+        elif recon_sampling_mode == 'presampled':
+            if verbose: p_header(f'>>> Prior sampling done ahead ...')
+            recon_period = (np.min(recon_yrs), np.max(recon_yrs))
+            vn = list(self.prior.keys())[0]
+            prior_time = self.prior[vn].da.time.values
+            self.prior_sample_years = prior_time
+            self.prior_sample_idx = list(range(len(prior_time)))
+            self.gen_Ye()
+            self.gen_Xb()
+            for yr_idx, target_yr in enumerate(tqdm(recon_yrs, desc='KF updating')):
+                self.Xa[yr_idx] = self.update_yr(target_yr, recon_loc_rad, recon_timescale, debug=debug, verbose=verbose, allownan=allownan)
 
         else:
             raise ValueError('Wrong recon_sampling_mode.')

cfr/psm.py

@@ -1252,7 +1252,7 @@ def __init__(self, pobj=None, climate_required=['sst', 'sss']):
         self.pobj = pobj
         self.climate_required = climate_required
 
-    def forward(self, seed=2333, species='all_sea', clean=1, pH=1, omega=1, sw=2, H=1, age=15):
+    def forward(self, seed=2333, species='all', clean=1, pH=1, omega=1, sw=2, H=1, age=15):
         vn1 = f'model.{self.climate_required[0]}'
         vn2 = f'model.{self.climate_required[0]}'
         res, sigma = pb.MgCa_forward(
@@ -1264,7 +1264,7 @@ def forward(self, seed=2333, species='all_sea', clean=1, pH=1, omega=1, sw=2, H=
         pp = ProxyRecord(
             pid=self.pobj.pid,
             time=self.pobj.clim[vn1].da.time,
-            value=pct[1],  # median
+            value=np.exp(pct[1]),  # median; the raw output is ln(Mg/Ca), so taking exp to recover Mg/Ca
             lat=self.pobj.lat,
             lon=self.pobj.lon,
             ptype=self.pobj.ptype,

cfr/reconjob.py

@@ -775,8 +775,10 @@ def save_recon(self, save_path, compress_params=None, verbose=False, output_full
                     iobw =  utils.geo_mean(da, lat_min=-20, lat_max=20, lon_min=40 ,lon_max=100)
                     ds['iobw'] = iobw
 
-            if not output_full_ens: da = da.mean(dim='ens')
-            ds[vn] = da
+            if output_full_ens:
+                ds[vn] = da
+            else:
+                ds[vn] = da.mean(dim='ens')
 
         encoding_dict = {}
         for k in self.recon_fields.keys():

setup.py

@@ -5,7 +5,7 @@
 
 setup(
     name='cfr',  # required
-    version='2024.6.26',
+    version='2024.7.1',
     description='cfr: a Python package for Climate Field Reconstruction',
     long_description=long_description,
     long_description_content_type='text/x-rst',