Changeset 27354 for issm/trunk-jpl/src/m/classes/stochasticforcing.py

Timestamp:

11/04/22 04:51:01 (2 years ago)

Author:

vverjans

Message:

NEW added autoregression and stochastic capabilities for subglacial discharge

File:

• issm/trunk-jpl/src/m/classes/stochasticforcing.py (modified) (6 diffs)

issm/trunk-jpl/src/m/classes/stochasticforcing.py

@@ -36 +36 @@
         s += '{}\n'.format(fielddisplay(self, 'randomflag', 'whether to apply real randomness (true) or pseudo-randomness with fixed seed (false)'))
         s += 'Available fields:\n'
+        s += '   BasalforcingsDeepwaterMeltingRatearma\n'
         s += '   BasalforcingsSpatialDeepwaterMeltingRate\n'
         s += '   DefaultCalving\n'
         s += '   FloatingMeltRate\n'
         s += '   FrictionWaterPressure\n'
+        s += '   FrictionCoulombWaterPressure\n'
+        s += '   FrictionSchoofWaterPressure\n'
         s += '   FrontalForcingsRignotarma (thermal forcing)\n'
+        s += '   FrontalForcingsSubglacialDischargearma\n'
         s += '   SMBarma\n'
         s += '   SMBforcing\n'
@@ -122 +126 @@
         indSMBarma   = -1  # About to check for index of SMBarma
         indTFarma    = -1  # About to check for index of FrontalForcingsRignotarma
+        indSdarma    = -1  # About to check for index of FrontalForcingsSubglacialDischargearma
         indBDWarma   = -1  # About to check for index of BasalforcingsDeepwaterMeltingRatearma
         if ('SMBarma' in self.fields):
@@ -133 +138 @@
             if(md.frontalforcings.arma_timestep<self.stochastictimestep):
                 raise TypeError('FrontalForcingsRignotarma cannot have a timestep shorter than stochastictimestep')
+        if ('FrontalForcingsSubglacialDischargearma' in self.fields):
+            indSdarma = self.fields.index('FrontalForcingsSubglacialDischargearma')  # Index of Sdarma, now check for consistency with other timesteps
+            dimensions[indSdarma] = md.frontalforcings.num_basins
+            if(md.frontalforcings.sd_arma_timestep<self.stochastictimestep):
+                raise TypeError('FrontalForcingsSubglacialDischargearma cannot have a timestep shorter than stochastictimestep')
         if ('BasalforcingsDeepwaterMeltingRatearma' in self.fields):
             indBDWarma = self.fields.index('BasalforcingsDeepwaterMeltingRatearma')  # Index of BDWarma, now check for consistency with other timesteps
@@ -152 +162 @@
             if((md.frontalforcings.arma_timestep != md.basalforcings.arma_timestep) and np.any(covsum != 0)):
                 raise IOError('FrontalForcingsRignotarma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
+        if (indSMBarma != -1 and indSdarma != -1):  # Both ARMA models are used: check ARMA time step consistency
+            covsum = self.covariance[np.sum(dimensions[0:indSMBarma]).astype(int):np.sum(dimensions[0:indSMBarma + 1]).astype(int), np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int)]
+            if((md.smb.arma_timestep != md.frontalforcings.sd_arma_timestep) and np.any(covsum != 0)):
+                raise IOError('SMBarma and FrontalForcingsSubglacialDischargearma have different arma_timestep and non-zero covariance')
+        if (indSdarma != -1 and indBDWarma != -1):  # Both ARMA models are used: check ARMA time step consistency
+            covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indBDWarma]).astype(int):np.sum(dimensions[0:indBDWarma + 1]).astype(int)]
+            if((md.frontalforcings.sd_arma_timestep != md.basalforcings.arma_timestep) and np.any(covsum != 0)):
+                raise IOError('FrontalForcingsSubglacialDischargearma and BasalforcingsDeepwaterMeltingRatearma have different arma_timestep and non-zero covariance')
+        if (indSdarma != -1 and indTFarma != -1):  # Both ARMA models are used: check ARMA time step consistency
+            covsum = self.covariance[np.sum(dimensions[0:indSdarma]).astype(int):np.sum(dimensions[0:indSdarma + 1]).astype(int), np.sum(dimensions[0:indTFarma]).astype(int):np.sum(dimensions[0:indTFarma + 1]).astype(int)]
+            if((md.frontalforcings.sd_arma_timestep != md.frontalforcings.arma_timestep) and np.any(covsum != 0)):
+                raise IOError('FrontalForcingsSubglacialDischargearma and FrontalForcingsRignotarma have different arma_timestep and non-zero covariance')
+
 
         md = checkfield(md, 'fieldname', 'stochasticforcing.isstochasticforcing', 'values', [0, 1])
@@ -187 +210 @@
                     dimensions[ind] = md.smb.num_basins
                 if (field == 'FrontalForcingsRignotarma'):
+                    dimensions[ind] = md.frontalforcings.num_basins
+                if (field == 'FrontalForcingsSubglacialDischargearma'):
                     dimensions[ind] = md.frontalforcings.num_basins
                 if (field == 'BasalforcingsDeepwaterMeltingRatearma'):
@@ -237 +262 @@
                      'FrictionSchoofWaterPressure': 'frictionschoof',
                      'FrontalForcingsRignotarma': 'frontalforcingsrignotarma',
+                     'FrontalForcingsSubglacialDischargearma': 'frontalforcingsrignotarma',
                      'SMBarma': 'SMBarma',
                      'SMBforcing': 'SMBforcing'}