From 5a9430fe79e8f3ce8f7488589be4822ec0aeca14 Mon Sep 17 00:00:00 2001
From: "Boyd, Oliver S" <olboyd@usgs.gov>
Date: Mon, 27 Jan 2020 17:55:18 +0000
Subject: [PATCH] Correct near line 214 to ensure MinIndex and VolFrac are
 defined for water and ice.

---
 python/GeoPhys.py | 20 ++++++++++++++++++--
 1 file changed, 18 insertions(+), 2 deletions(-)

diff --git a/python/GeoPhys.py b/python/GeoPhys.py
index 4bc7124..66a2e02 100644
--- a/python/GeoPhys.py
+++ b/python/GeoPhys.py
@@ -99,6 +99,9 @@ def loadPar(x1,x2,y1,y2):
     mi = nc.variables['Mineral Index'][:]
     vf = nc.variables['Volume Fraction'][:]
 
+    MNF1[210:211] = 0
+    MNF2[:,210:211] = -5
+
     return MNF1,MNF2,WTD,mi,vf
 
 def CalcGeo(lat,lon,Depth,igt):
@@ -134,15 +137,23 @@ def CalcGeo(lat,lon,Depth,igt):
         US = np.zeros((212),dtype=np.float64)
         PS = np.zeros((212),dtype=np.float64)
         nm,nmi = np.shape(mi)
+        P = 1e6*np.ones((1),dtype=np.float64)
+        T = 350*np.ones((1),dtype=np.float64)
         for l in range(0,nmi):
             m = np.nonzero(mi[:,l] > 0)[0]
             if np.size(m) > 0:
-                P = 1e6*np.ones((1),dtype=np.float64)
-                T = 350*np.ones((1),dtype=np.float64)
                 Ks, Us, E, lmb, v, Vp, Vs, Ps, Vpv, Vpr, Vsv, Vsr, a = mv.CalcMV(mi[m,l],vf[m,l],T,P)
                 KS[l] = Ks
                 US[l] = Us
                 PS[l] = Ps
+            Ks, Us, E, lmb, v, Vp, Vs, Ps, Vpv, Vpr, Vsv, Vsr, a = mv.CalcMV(mi[0,210],vf[0,210],T,P)
+            KS[210] = Ks
+            US[210] = Us
+            PS[210] = Ps
+            Ks, Us, E, lmb, v, Vp, Vs, Ps, Vpv, Vpr, Vsv, Vsr, a = mv.CalcMV(mi[0,211],vf[0,211],T,P)
+            KS[211] = Ks
+            US[211] = Us
+            PS[211] = Ps
         for l in range(0,nLocs):
             Temp[l,:] = 300 + 0.025*Depth
     else:
@@ -200,6 +211,11 @@ def CalcGeo(lat,lon,Depth,igt):
                 m = np.nonzero(mi[:,PG[l]-1] > 0)[0]
                 MinIndex = mi[m,PG[l]-1]
                 VolFrac = vf[m,PG[l]-1]
+                if PG[l] > 210:
+                    MinIndex = np.ones((1,1))
+                    VolFrac = np.ones((1,1))
+                    MinIndex[0] = mi[0,PG[l]-1]
+                    VolFrac[0] = vf[0,PG[l]-1]
                 f = interp.interp1d(Depth,Temp[i,:],kind='linear',fill_value='extrapolate')
                 Vecs[2,:] = f(Dep)
 
-- 
GitLab