diff --git a/+gemini3d/+setup/+gridgen/makegrid_cart_3D.m b/+gemini3d/+setup/+gridgen/makegrid_cart_3D.m
index 6972c5e3..323b64c6 100644
--- a/+gemini3d/+setup/+gridgen/makegrid_cart_3D.m
+++ b/+gemini3d/+setup/+gridgen/makegrid_cart_3D.m
@@ -1,8 +1,15 @@
-function xgf = makegrid_cart_3D(p)
+function xgf = makegrid_cart_3D(p,year)
 %% make 3D grid and save to disk
+
+ if ~exist('year','var')
+     % If year is not specified, set it 1985 by default
+      year = 1985;
+ end
+ 
 arguments
   p (1,1) struct
 end
+
 %% create altitude grid
 %p.alt_min = 80e3;
 %p.alt_max = 1000e3;
@@ -63,7 +70,7 @@
 %DISTANCE EW AND NS (FROM ENU (or UEN in our case - cyclic permuted) COORD. SYSTEM) NEED TO BE CONVERTED TO DIPOLE SPHERICAL AND THEN
 %GLAT/GLONG - BASICALLY HERE WE ARE MAPPING THE CARTESIAN GRID ONTO THE
 %SURFACE OF A SPHERE THEN CONVERTING TO GEOGRAPHIC.
-[thetactr,phictr] = gemini3d.geog2geomag(p.glat, p.glon);    %get the magnetic coordinates of the grid center, based on user input
+[thetactr,phictr] = gemini3d.geog2geomag(p.glat, p.glon, year);    %get the magnetic coordinates of the grid center, based on user input
 
 %% Center of earth distance
 r=Re+z;
@@ -84,7 +91,7 @@
 phi=repmat(phi,[lx1,1,lx3]);
 
 %% COMPUTE THE GEOGRAPHIC COORDINATES OF EACH GRID POINT
-[glatgrid,glongrid] = gemini3d.geomag2geog(theta,phi);
+[glatgrid,glongrid] = gemini3d.geomag2geog(theta,phi,year);
 
 %% COMPUTE ECEF CARTESIAN IN CASE THEY ARE NEEDED
 xECEF=r.*sin(theta).*cos(phi);
diff --git a/+gemini3d/+setup/+gridgen/makegrid_tilteddipole_3D.m b/+gemini3d/+setup/+gridgen/makegrid_tilteddipole_3D.m
index fcaf18d4..0c4e916d 100644
--- a/+gemini3d/+setup/+gridgen/makegrid_tilteddipole_3D.m
+++ b/+gemini3d/+setup/+gridgen/makegrid_tilteddipole_3D.m
@@ -1,4 +1,4 @@
-function xgf = makegrid_tilteddipole_3D(dtheta,dphi,lpp,lqp,lphip,altmin,glat,glon,gridflag)
+function xgf = makegrid_tilteddipole_3D(dtheta,dphi,lpp,lqp,lphip,altmin,glat,glon,gridflag,year)
 
 %NOTE THAT INPUTS DTHETA AND DPHI ARE INTENDED TO REPRESENT THE FULL THETA
 %AND PHI EXTENTS OF
@@ -13,11 +13,18 @@
 %   properly.
 %(5) Terrestrial mag. moment and mass are hard-coded in.
 %(6) We assume target latitude is northern hemisphere
+%(7) 
 %
 %OTHER NOTES:
 % - Fortran code will interpret this grid as including ghost cells, so
 % if you want a dimension to be size "n" adjust requested grid size so that
 % it is "n+4"
+
+ if ~exist('year','var')
+     % If year is not specified, set it 1985 by default
+      year = 1985;
+ end
+ 
 arguments
     dtheta {mustBeNumeric}
     dphi {mustBeNumeric}
@@ -28,6 +35,7 @@
     glat {mustBeNumeric}
     glon {mustBeNumeric}
     gridflag {mustBeInteger}
+    year {mustBeInteger}
 end
 
 %% PAD GRID WITH GHOST CELLS
@@ -42,7 +50,7 @@
 
 
 %TD SPHERICAL LOCATION OF REQUESTED CENTER POINT
-[thetatd,phid]=gemini3d.geog2geomag(glat,glon);
+[thetatd,phid]=gemini3d.geog2geomag(glat,glon,year);
 
 
 %SETS THE EDGES OF THE GRID
@@ -368,7 +376,7 @@
 
 %xg.glat=(pi/2-theta)*180/pi; xg.glon=phi*180/pi*ones(lx(1),lx(2));
 for iphi=1:lphi
-  [glats,glons]=gemini3d.geomag2geog(xg.theta(:,:,iphi),xg.phi(1,1,iphi)*ones(lq,lp));
+  [glats,glons]=gemini3d.geomag2geog(xg.theta(:,:,iphi),xg.phi(1,1,iphi)*ones(lq,lp),year);
   % only meant to work for one latitude at at time
   xg.glat(:,:,iphi)=glats;
   xg.glon(:,:,iphi)=glons;
diff --git a/+gemini3d/+setup/plot_mapgrid.m b/+gemini3d/+setup/plot_mapgrid.m
index db0a5538..4797b175 100644
--- a/+gemini3d/+setup/plot_mapgrid.m
+++ b/+gemini3d/+setup/plot_mapgrid.m
@@ -1,4 +1,9 @@
-function ha = plot_mapgrid(xg,flagsource,neuinfo)
+function ha = plot_mapgrid(xg,flagsource,neuinfo,year)
+
+ if ~exist('year','var')
+     % If year is not specified, set it 1985 by default
+      year = 1985;
+ end
 
 arguments
   xg (1,1) struct
@@ -12,7 +17,6 @@
 dmlon=max(mlon(:))-min(mlon(:));
 dmlat=max(mlat(:))-min(mlat(:));
 
-
 %% ORGANIZE INPUT STRUCTURE
 if flagsource ~= 0
     neugridtype=neuinfo.neugridtype;
@@ -29,7 +33,7 @@
         rhomax=neuinfo.rhomax;
     end %if
 
-    [sourcetheta,sourcephi]=geog2geomag(sourcelat,sourcelong);
+    [sourcetheta,sourcephi]=geog2geomag(sourcelat,sourcelong,year);
     sourcemlat=90-sourcetheta*180/pi;
     sourcemlon=sourcephi*180/pi;
 
@@ -421,7 +425,7 @@
         hold on;
         %ax=axis;
         %load coastlines;
-        [thetacoast,phicoast]=geog2geomag(coastlat,coastlon);
+        [thetacoast,phicoast]=geog2geomag(coastlat,coastlon,year);
         mlatcoast=90-thetacoast*180/pi;
         mloncoast=phicoast*180/pi;
 
diff --git a/+gemini3d/+tests/TestAUnit.m b/+gemini3d/+tests/TestAUnit.m
index 88a1444e..265320d2 100644
--- a/+gemini3d/+tests/TestAUnit.m
+++ b/+gemini3d/+tests/TestAUnit.m
@@ -119,16 +119,16 @@ function test_max_mpi(tc)
 
 function test_coord(tc)
 
-[lat, lon] = gemini3d.geomag2geog(pi/2, pi/2);
+[lat, lon] = gemini3d.geomag2geog(pi/2, pi/2, 2020);
 tc.verifyEqual([lat, lon], [0, 19], 'AbsTol', 1e-6, 'RelTol', 0.001)
 
-[theta, phi] = gemini3d.geog2geomag(0,0);
+[theta, phi] = gemini3d.geog2geomag(0,0, 2020);
 tc.verifyEqual([theta, phi], [1.50863496978059, 1.24485046147953], 'AbsTol', 1e-6, 'RelTol', 0.001)
 
-[alt, lon, lat] = gemini3d.UEN2geog(0,0,0, pi/2, pi/2);
+[alt, lon, lat] = gemini3d.UEN2geog(0,0,0, pi/2, pi/2, 2020);
 tc.verifyEqual([alt, lat, lon], [0, 0, 19], 'AbsTol', 1e-6, 'RelTol', 0.001)
 
-[z,x,y] = gemini3d.geog2UEN(0, 0, 0, pi/2, pi/2);
+[z,x,y] = gemini3d.geog2UEN(0, 0, 0, pi/2, pi/2, 2020);
 tc.verifyEqual([z,x,y], [0, -2076275.16205889, 395967.844181141], 'AbsTol', 1e-6, 'RelTol', 0.001)
 
 end
diff --git a/+gemini3d/UEN2geog.m b/+gemini3d/UEN2geog.m
index 4cf22c6d..aacc2d1a 100644
--- a/+gemini3d/UEN2geog.m
+++ b/+gemini3d/UEN2geog.m
@@ -1,10 +1,11 @@
-function [alt,glon,glat] = UEN2geog(z,x,y,thetactr,phictr)
+function [alt,glon,glat] = UEN2geog(z,x,y,thetactr,phictr,year)
 arguments
   z {mustBeNumeric}
   x {mustBeNumeric}
   y {mustBeNumeric}
   thetactr {mustBeNumeric}
   phictr {mustBeNumeric}
+  year {mustBeNumeric}
 end
 
 % converts magnetic up, north, east coordinates into geographic
@@ -43,6 +44,6 @@
 
 
 %Now convert the magnetic to geographic using our simple transformation
-[glat,glon] = gemini3d.geomag2geog(theta,phi);
+[glat,glon] = gemini3d.geomag2geog(theta,phi,year);
 
 end
diff --git a/+gemini3d/geog2UEN.m b/+gemini3d/geog2UEN.m
index 132e1a89..98079623 100644
--- a/+gemini3d/geog2UEN.m
+++ b/+gemini3d/geog2UEN.m
@@ -1,4 +1,4 @@
-function [z,x,y] = geog2UEN(alt,glon,glat,thetactr,phictr)
+function [z,x,y] = geog2UEN(alt,glon,glat,thetactr,phictr,year)
 arguments
   alt {mustBeNumeric}
   glon {mustBeNumeric}
@@ -16,7 +16,7 @@
 z=alt;
 
 %Convert to geomganetic coordinates
-[theta,phi] = gemini3d.geog2geomag(glat,glon);
+[theta,phi] = gemini3d.geog2geomag(glat,glon,year);
 
 
 %Convert to northward distance in meters
diff --git a/+gemini3d/geog2geomag.m b/+gemini3d/geog2geomag.m
index 136abd2c..9e8434fd 100644
--- a/+gemini3d/geog2geomag.m
+++ b/+gemini3d/geog2geomag.m
@@ -1,11 +1,45 @@
-function [thetat,phit] = geog2geomag(lat,lon)
+function [thetat,phit] = geog2geomag(lat,lon,year)
+
+ if ~exist('year','var')
+     % If year is not specified, set it 1985 by default
+      year = 1985;
+ end
+ 
+ 
 arguments
   lat {mustBeNumeric}
   lon {mustBeNumeric}
+  year {mustBeNumeric}
 end
 
-thetan = deg2rad(11);
-phin = deg2rad(289);
+% Schmidt semi-normalised spherical harmonic coefficients are from http://wdc.kugi.kyoto-u.ac.jp/igrf/coef/igrf13coeffs.html
+% Calculation of thetan and phin are based on Millward et al., 1996
+% For year < 1985, use the same coefficients as for 1985
+case ismember(year, 1920:1987)
+thetan=deg2rad(1.1000e+01);
+phin=deg2rad(2.8900e+02);
+case ismember(year, 1988:1992)
+thetan=deg2rad(1.0862e+01);
+phin=deg2rad(2.8887e+02);
+case ismember(year, 1993:1997)
+thetan=deg2rad(1.0677e+01);
+phin=deg2rad(2.8858e+02);
+case ismember(year, 1998:2002)
+thetan=deg2rad(1.0457e+01);
+phin=deg2rad(2.8843e+02);
+case ismember(year, 2003:2007)
+thetan=deg2rad(1.0252e+01);
+phin=deg2rad(2.8819e+02);
+case ismember(year, 2008:2012)
+thetan=deg2rad(9.9840e+00);
+phin=deg2rad(2.8779e+02);
+case ismember(year, 2013:2017)
+thetan=deg2rad(9.6869e+00);
+phin=deg2rad(2.8739e+02);
+case ismember(year, 2018:2025)
+thetan=deg2rad(9.4105e+00);
+phin=deg2rad(2.8732e+02);
+end
 
 %enforce [0,360] longitude
 loncorrected = mod(lon, 360);
diff --git a/+gemini3d/geomag2geog.m b/+gemini3d/geomag2geog.m
index aeef9a63..c01e67be 100644
--- a/+gemini3d/geomag2geog.m
+++ b/+gemini3d/geomag2geog.m
@@ -1,11 +1,44 @@
-function [lat,lon] = geomag2geog(thetat,phit)
+function [lat,lon] = geomag2geog(thetat,phit,year)
+
+ if ~exist('year','var')
+     % If year is not specified, set it 1985 by default
+      year = 1985;
+ end
+ 
 arguments
   thetat {mustBeNumeric}
   phit {mustBeNumeric}
+  year {mustBeNumeric}
 end
 
-thetan = deg2rad(11);
-phin = deg2rad(289);
+% Schmidt semi-normalised spherical harmonic coefficients are from http://wdc.kugi.kyoto-u.ac.jp/igrf/coef/igrf13coeffs.html
+% Calculation of thetan and phin are based on Millward et al., 1996
+% For year < 1985, use the same coefficients as for 1985
+case ismember(year, 1920:1987)
+thetan=deg2rad(1.1000e+01);
+phin=deg2rad(2.8900e+02);
+case ismember(year, 1988:1992)
+thetan=deg2rad(1.0862e+01);
+phin=deg2rad(2.8887e+02);
+case ismember(year, 1993:1997)
+thetan=deg2rad(1.0677e+01);
+phin=deg2rad(2.8858e+02);
+case ismember(year, 1998:2002)
+thetan=deg2rad(1.0457e+01);
+phin=deg2rad(2.8843e+02);
+case ismember(year, 2003:2007)
+thetan=deg2rad(1.0252e+01);
+phin=deg2rad(2.8819e+02);
+case ismember(year, 2008:2012)
+thetan=deg2rad(9.9840e+00);
+phin=deg2rad(2.8779e+02);
+case ismember(year, 2013:2017)
+thetan=deg2rad(9.6869e+00);
+phin=deg2rad(2.8739e+02);
+case ismember(year, 2018:2025)
+thetan=deg2rad(9.4105e+00);
+phin=deg2rad(2.8732e+02);
+end
 
 %enforce phit = [0,2pi]
 phit = mod(phit, 2*pi);
diff --git a/Examples/TECplot_map.m b/Examples/TECplot_map.m
index 609ec88a..178fd5f9 100644
--- a/Examples/TECplot_map.m
+++ b/Examples/TECplot_map.m
@@ -68,7 +68,7 @@ function TECplot_map(direc)
 
     %ADD A MAP OF COASTLINES
     load("coastlines", "coastlat", "coastlon")
-    [thetacoast,phicoast] = gemini3d.geog2geomag(coastlat,coastlon);
+    [thetacoast,phicoast] = gemini3d.geog2geomag(coastlat,coastlon,year);
     mlatcoast=90-thetacoast*180/pi;
     mloncoast=phicoast*180/pi;
 
diff --git a/Examples/grid_examples.m b/Examples/grid_examples.m
index 5639ff3c..49ea1c0b 100644
--- a/Examples/grid_examples.m
+++ b/Examples/grid_examples.m
@@ -55,7 +55,7 @@ function grid_examples()
 gridflag=1;
 flagsource=0;
 iscurv=true;
-
+year = 2020; % Adjust magnetic coordinates in accordance to the year required
 
 %% GEOGRAPHIC COORDINATES OF NEUTRAL SOURCE (OR GRID CENTER)
 neuinfo.sourcelat=[];
@@ -84,7 +84,7 @@ function grid_examples()
 
 
 %% FOR USERS INFO CONVERT SOURCE LOCATION TO GEOMAG
-[sourcetheta,sourcephi]=gemini3d.geog2geomag(neuinfo.sourcelat,neuinfo.sourcelong);
+[sourcetheta,sourcephi]=gemini3d.geog2geomag(neuinfo.sourcelat,neuinfo.sourcelong,year);
 sourcemlat=90-sourcetheta*180/pi;
 sourcemlon=sourcephi*180/pi;
 
@@ -92,7 +92,7 @@ function grid_examples()
 %% RUN THE GRID GENERATION CODE
 if ~exist('xg', 'var')
   if iscurv
-    xg = gemini3d.setup.gridgen.makegrid_tilteddipole_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag);
+    xg = gemini3d.setup.gridgen.makegrid_tilteddipole_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag,year);
 %    xg=makegrid_tilteddipole_nonuniform_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag);
 %    xg=makegrid_tilteddipole_nonuniform_oneside_3D(dtheta,dphi,lp,lq,lphi,altmin,glat,glon,gridflag);
   else
diff --git a/Examples/magplot_fort_map.m b/Examples/magplot_fort_map.m
index 824b5fcb..f4c21a4a 100644
--- a/Examples/magplot_fort_map.m
+++ b/Examples/magplot_fort_map.m
@@ -205,7 +205,7 @@
     %ADD A MAP OF COASTLINES
 %    if (license('test','Map_Toolbox'))
         load coastlines;
-        [thetacoast,phicoast]= gemini3d.geog2geomag(coastlat,coastlon);
+        [thetacoast,phicoast]= gemini3d.geog2geomag(coastlat,coastlon,year);
         mlatcoast=90-thetacoast*180/pi;
         mloncoast=phicoast*180/pi;