diff --git a/da_functions.py b/da_functions.py
index 97299fd8c783994116ab2f88f2799a8cb12ebd2a..06e0e509638481527bf0098dd650306b3f80a156 100644
--- a/da_functions.py
+++ b/da_functions.py
@@ -262,8 +262,7 @@ def run_linear_advection_EnKF(m_const,da_const):
 
     #Construct localization matrix C if loc!=None
     if da_const["loc"]:
-        if da_const["loc_type"]=='gaussian':
-            C = gaussian_loc_matrix(da_const,m_const)
+        C = loc_matrix(da_const,m_const)
     else:
         C=np.ones([m_const["nx"],m_const["nx"]])
 
@@ -712,12 +711,35 @@ def var_reduction_estimate_iterative(states,m_const,da_const,j=0,dJ_update_flag=
 
 
 
-def gaussian_loc_matrix(da_const,m_const):
+def loc_matrix(da_const,m_const):
+    """
+    Creates the localization matrix, using either a gaussian or gaspari cohn function 
+    """
     C = np.zeros([m_const['nx'],m_const['nx']])
-    #The gaussian is slightly less trivial because of the periodic boundary domain, so what I do is I use the same function that i use to generate the initial conditions and then step it once for all other locations
-    #And of course i am not sure if i have my dimensions switched
-    C[:,0] = gaussian_initial_condition(m_const["x_grid"],da_const["loc_length"])
+    # I cheat a bit by mirroring the functions to accomodate for the repeating boundary conditions, but this should only lead to a maximum dx/2 error. 
+    if da_const['loc_type']=='gaussian': 
+        C[:,0] = gaussian_initial_condition(m_const["x_grid"],da_const["loc_length"])
+    if da_const['loc_type']=='gaspari_cohn': 
+        C[:,0] = gaspari_cohn(m_const["x_grid"],da_const["loc_length"])
+
     for i in range(1,m_const['nx']):
         C[:,i] = np.hstack([C[-1,i-1],C[:-1,i-1]])
     return C
 
+def gaspari_cohn(x,loc_length):
+    """Gaspari-Cohn function."""
+    
+    ra = x/loc_length
+    gp = np.zeros_like(ra)
+    i=np.where(ra<=1.)[0]
+    gp[i]=-0.25*ra[i]**5+0.5*ra[i]**4+0.625*ra[i]**3-5./3.*ra[i]**2+1.
+    i=np.where((ra>1.)*(ra<=2.))[0]
+    gp[i]=1./12.*ra[i]**5-0.5*ra[i]**4+0.625*ra[i]**3+5./3.*ra[i]**2-5.*ra[i]+4.-2./3./ra[i]
+        
+    #Now we mirror things to zero for the periodic boundary domain
+    half_idx = int(ra.shape[0]/2)
+    gp[-half_idx:] = gp[half_idx-1::-1]
+    return gp
+
+
+