diff --git a/components/isceobj/Alos2Proc/runFrameOffset.py b/components/isceobj/Alos2Proc/runFrameOffset.py
index 635e94c9..59356688 100644
--- a/components/isceobj/Alos2Proc/runFrameOffset.py
+++ b/components/isceobj/Alos2Proc/runFrameOffset.py
@@ -97,18 +97,23 @@ def frameOffset(track, image, outputfile, crossCorrelation=True, matchingMode=0)
         swath1 = track.frames[j-1].swaths[0]
         swath2 = track.frames[j].swaths[0]
 
+        #consider frame/swath azimuth sensing start differences caused by swath mosaicking
+        # tested with alos2App.py, alos2burstApp.py and alosStack
+        frame1 = track.frames[j-1]
+        frame2 = track.frames[j]
+        delta_az = -((frame2.sensingStart - frame1.sensingStart).total_seconds() - (swath2.sensingStart - swath1.sensingStart).total_seconds()) / swath1.azimuthLineInterval
 
         #offset from geometry
         offsetGeometrical = computeFrameOffset(swath1, swath2)
         rangeOffsetGeometrical.append(offsetGeometrical[0])
-        azimuthOffsetGeometrical.append(offsetGeometrical[1])
+        azimuthOffsetGeometrical.append(offsetGeometrical[1]+delta_az)
 
         #offset from cross-correlation
         if crossCorrelation:
             offsetMatching = estimateFrameOffset(swath1, swath2, image1, image2, matchingMode=matchingMode)
             if offsetMatching != None:
                 rangeOffsetMatching.append(offsetMatching[0])
-                azimuthOffsetMatching.append(offsetMatching[1])
+                azimuthOffsetMatching.append(offsetMatching[1]+delta_az)
             else:
                 print('******************************************************************')
                 print('WARNING: bad matching offset, we are forced to use')
diff --git a/contrib/stack/topsStack/README.md b/contrib/stack/topsStack/README.md
index 54b6ac7f..75f4454b 100644
--- a/contrib/stack/topsStack/README.md
+++ b/contrib/stack/topsStack/README.md
@@ -186,7 +186,7 @@ Ionospheric phase estimation has more requirements than regular InSAR processing
 In stack ionospheric phase estimation, acquistions with same swath starting ranges are put in a group. A network is formed within a group. Extra pairs are also processed to connect the different groups so that all acquistions are connected. But we need to estimate a phase offset for these extra pairs, which might not be accurate. Therefore, for a particualr swath starting ranges, if there are only a few acquistions, it's better to just discard them so that we don't have to estimate the phase offsets.
 
 ```
-s1_select_ion.py -dir data/slc -sn 33.550217/37.119545 -nr 10
+s1_select_ion.py -dir data/slc -sn 33.550217 37.119545 -nr 10
 ```
 
 Acquistions to be used need to fully cover the south/north bounds. After running this command, acquistion not to be used will be put in a folder named 'not_used'. It's OK to run this command multiple times.
@@ -204,6 +204,8 @@ An example --param_ion file 'ion_param.txt' is provided in the code directory. F
 stackSentinel.py -s ../data/slc -d ../data/dem/dem_1_arcsec/demLat_N32_N41_Lon_W113_W107.dem.wgs84 -b '33.550217 37.119545 -111.233932 -107.790451' -a ../data/s1_aux_cal -o ../data/orbit -C geometry -c 2 --param_ion ../code/ion_param.txt --num_connections_ion 3
 ```
 
+Note in 'ion_param.txt', if 'consider burst properties in ionosphere computation' is True, Coregistration options '-C', '--coregistration' in stackSentinel.py must be NESD.
+
 If ionospheric phase estimation is enabled in stackSentinel.py, it will generate the following run files. Here ***ns*** means number of steps in the original stack processing, which depends on the type of stack (slc, correlation, interferogram, and offset).
 
 -	run_ns+1_subband_and_resamp
@@ -251,6 +253,22 @@ Estimate ionospheric phase for each date. We highly recommend inspecting all pai
 
 Typical anomalies include dense fringes caused by phase unwrapping errors, and a range ramp as a result of errors in estimating phase offsets for pairs with different swath starting ranges (check pairs_diff_starting_ranges.txt).
 
+**run_ns+9_filtIonShift**
+
+Filter azimuth ionospheric shift.
+
+**run_ns+10_invertIonShift**
+
+Estimate azimuth ionospheric shift for each date. As in step **run_ns+8_invertIon**, check if there are anamolies.
+
+**run_ns+11_burstRampIon**
+
+Compute azimuth burst ramps as a result of ionosphere for each date.
+
+**run_ns+12_mergeBurstRampIon**
+
+Merge azimuth burst ramps.
+
 #### 3. run command files generated ####
 
 Run the commands sequentially.
@@ -261,7 +279,11 @@ Results from ionospheric phase estimation.
 
 -	reference and coreg_secondarys: now contains also subband burst SLCs
 -	ion: original ionospheric phase estimation results
+-       ion_azshift_dates: azimuth ionospheric shift for each acquistion
+-       ion_burst_ramp_dates: azimuth burst ramps caused by ionosphere for each acquistion
+-       ion_burst_ramp_merged_dates: merged azimuth burst ramps caused by ionosphere for each acquistion
 -	ion_dates: ionospheric phase for each acquistion
+-	ion/date1_date2/ion_cal/azshift.ion: azimuth ionospheric shift
 -	ion/date1_date2/ion_cal/filt.ion: filtered ionospheric phase
 -	ion/date1_date2/ion_cal/raw_no_projection.ion: original ionospheric phase
 -	ion/date1_date2/lower/merged/fine_look.unw: unwrapped lower band interferogram
@@ -273,6 +295,7 @@ If ionospheric phase estimation processing is swath by swath because of differen
 -	ion/date1_date2/lower/merged_IW*
 -	ion/date1_date2/upper/merged_IW*
 
+Unit of azimuth ionospheric shift is number of single look azimuth lines.
 After processing, we can plot ionospheric phase estimation results using plotIonPairs.py and plotIonDates.py. For example
 
 ```
@@ -285,4 +308,5 @@ Relationships of the ionospheric phases:
 ```
 ion_dates/date1.ion - ion_dates/date2.ion = ion/date1_date2/ion_cal/filt.ion
 ion_dates/date1.ion - ion_dates/date2.ion = ionospheric phase in merged/interferograms/date1_date2/filt_fine.unw
+ion_burst_ramp_merged_dates/date1.float - ion_burst_ramp_merged_dates/date2.float = azimuth burst ramps caused by ionosphere in merged/interferograms/date1_date2/filt_fine.unw
 ```
diff --git a/contrib/stack/topsStack/Stack.py b/contrib/stack/topsStack/Stack.py
index 4fd860f5..580ac6c7 100644
--- a/contrib/stack/topsStack/Stack.py
+++ b/contrib/stack/topsStack/Stack.py
@@ -339,8 +339,33 @@ def filtIon(self, function):
         self.f.write('win_min : ' + '{}'.format(self.win_min) + '\n')
         self.f.write('win_max : ' + '{}'.format(self.win_max) + '\n')
 
+    def filtIonShift(self, function):
+        self.f.write('###################################'+'\n')
+        self.f.write(function + '\n')
+        self.f.write('filtIonShift : ' + '\n')
+        self.f.write('reference_stack : ' + self.reference_stack + '\n')
+        self.f.write('reference : ' + self.reference + '\n')
+        self.f.write('secondary : ' + self.secondary + '\n')
+        self.f.write('input : ' + self.input + '\n')
+        self.f.write('coherence : ' + self.coherence + '\n')
+        self.f.write('output : ' + self.output + '\n')
+        self.f.write('win_min : ' + '{}'.format(self.win_min) + '\n')
+        self.f.write('win_max : ' + '{}'.format(self.win_max) + '\n')
+        self.f.write('nrlks : ' + '{}'.format(self.nrlks) + '\n')
+        self.f.write('nalks : ' + '{}'.format(self.nalks) + '\n')
 
 
+    def burstRampIon(self, function):
+        self.f.write('###################################'+'\n')
+        self.f.write(function + '\n')
+        self.f.write('burstRampIon : ' + '\n')
+        self.f.write('reference_stack : ' + self.reference_stack + '\n')
+        self.f.write('input : ' + self.input + '\n')
+        self.f.write('output : ' + self.output + '\n')
+        self.f.write('nrlks : ' + '{}'.format(self.nrlks) + '\n')
+        self.f.write('nalks : ' + '{}'.format(self.nalks) + '\n')
+        self.f.write('ion_height : ' + '{}'.format(self.ion_height) + '\n')
+
 
     def write_wrapper_config2run_file(self, configName, line_cnt, numProcess = 1):
         # dispassionate list of commands for single process
@@ -1330,7 +1355,7 @@ def unwrap_ion(self, pairs_same_starting_ranges_update, pairs_diff_starting_rang
                     configObj.defoMax = '2'
                     configObj.rangeLooks = '{}'.format(ionParamUsrObj.ION_numberRangeLooks0)
                     configObj.azimuthLooks = '{}'.format(ionParamUsrObj.ION_numberAzimuthLooks0)
-                    configObj.rmfilter = False
+                    configObj.rmFilter = False
                     configObj.unwMethod = 'snaphu'
                     configObj.unwrap(function)
 
@@ -1348,7 +1373,7 @@ def unwrap_ion(self, pairs_same_starting_ranges_update, pairs_diff_starting_rang
                     configObj.defoMax = '2'
                     configObj.rangeLooks = '{}'.format(ionParamUsrObj.ION_numberRangeLooks0)
                     configObj.azimuthLooks = '{}'.format(ionParamUsrObj.ION_numberAzimuthLooks0)
-                    configObj.rmfilter = False
+                    configObj.rmFilter = False
                     configObj.unwMethod = 'snaphu'
                     configObj.unwrap(function)
 
@@ -1542,6 +1567,105 @@ def invertIon(self):
         self.runf.write(self.text_cmd + cmd + '\n')
 
 
+    def filtIonShift(self, pairs):
+
+        ionParamUsrObj = ionParamUsr(self.param_ion)
+        ionParamUsrObj.configure()
+
+        line_cnt = 0
+        for p in pairs:
+            configName = os.path.join(self.config_path,'config_filtIonShift_{}_{}'.format(p[0], p[1]))
+            configObj = config(configName)
+            configObj.configure(self)
+            configObj.reference_stack = os.path.join(self.work_dir, 'reference')
+            configObj.reference = os.path.join(self.work_dir, 'coreg_secondarys', '{}'.format(p[0]))
+            configObj.secondary = os.path.join(self.work_dir, 'coreg_secondarys', '{}'.format(p[1]))
+            configObj.input = os.path.join(self.work_dir, 'ion', '{}_{}'.format(p[0], p[1]), 'ion_cal', 'filt.ion')
+            configObj.coherence = os.path.join(self.work_dir, 'ion', '{}_{}'.format(p[0], p[1]), 'ion_cal', 'raw_no_projection.cor')
+            configObj.output = os.path.join(self.work_dir, 'ion', '{}_{}'.format(p[0], p[1]), 'ion_cal', 'azshift.ion')
+            configObj.win_min = ionParamUsrObj.ION_ionshiftFilteringWinsizeMin
+            configObj.win_max = ionParamUsrObj.ION_ionshiftFilteringWinsizeMax
+            configObj.nrlks = ionParamUsrObj.ION_numberRangeLooks
+            configObj.nalks = ionParamUsrObj.ION_numberAzimuthLooks
+            configObj.filtIonShift('[Function-1]')
+            configObj.finalize()
+
+            line_cnt += 1
+            #line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt, self.numProcess)
+            line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt)
+            del configObj
+
+
+    def invertIonShift(self):
+
+        ionParamUsrObj = ionParamUsr(self.param_ion)
+        ionParamUsrObj.configure()
+        
+        ion_in = os.path.join(self.work_dir,'ion')
+        ion_out = os.path.join(self.work_dir,'ion_azshift_dates')
+        hgt = os.path.join(self.work_dir,'merged/geom_reference/hgt.rdr')
+
+        cmd = 'invertIon.py --idir {} --filename azshift.ion --odir {} --nrlks1 {} --nalks1 {} --nrlks2 {} --nalks2 {} --merged_geom {} --interp --msk_overlap'.format(ion_in,ion_out,ionParamUsrObj.ION_numberRangeLooks, ionParamUsrObj.ION_numberAzimuthLooks, self.rangeLooks, self.azimuthLooks,hgt)
+
+        self.runf.write(self.text_cmd + cmd + '\n')
+
+
+    def burstRampIon(self, dates):
+
+        ionParamUsrObj = ionParamUsr(self.param_ion)
+        ionParamUsrObj.configure()
+
+        line_cnt = 0
+        for p in dates:
+            configName = os.path.join(self.config_path,'config_burstRampIon_{}'.format(p))
+            configObj = config(configName)
+            configObj.configure(self)
+            configObj.reference_stack = os.path.join(self.work_dir, 'reference')
+            configObj.input = os.path.join(self.work_dir, 'ion_azshift_dates', '{}.ion'.format(p))
+            configObj.output = os.path.join(self.work_dir, 'ion_burst_ramp_dates', '{}'.format(p))
+            configObj.nrlks = self.rangeLooks
+            configObj.nalks = self.azimuthLooks
+            configObj.ion_height = ionParamUsrObj.ION_ionHeight
+            configObj.burstRampIon('[Function-1]')
+            configObj.finalize()
+
+            line_cnt += 1
+            #line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt, self.numProcess)
+            line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt)
+            del configObj
+
+
+    def mergeBurstRampIon(self, dates, stackReferenceDate):
+
+        line_cnt = 0
+        for d in dates:
+            configName = os.path.join(self.config_path,'config_mergeBurstRampIon_' + d)
+            configObj = config(configName)
+            configObj.configure(self)
+            configObj.stack = os.path.join(self.work_dir, 'stack')
+            if d == stackReferenceDate:
+                configObj.reference = os.path.join(self.work_dir, 'reference')
+            else:
+                configObj.reference = os.path.join(self.work_dir, 'coreg_secondarys/' + d)
+            configObj.dirName = os.path.join(self.work_dir, 'ion_burst_ramp_dates/' + d)
+            configObj.namePattern = 'burst*float'
+            configObj.mergedFile = os.path.join(self.work_dir, 'ion_burst_ramp_merged_dates/' + d + '.float')
+            configObj.mergeBurstsMethod = 'top'
+            if d == stackReferenceDate:
+                configObj.aligned = 'False'
+            else:
+                configObj.aligned = 'True'
+            configObj.validOnly = 'True'
+            configObj.useVirtualFiles = 'True'
+            configObj.multiLook = 'True'
+            configObj.mergeBurst('[Function-1]')
+            configObj.finalize()
+
+            line_cnt += 1
+            line_cnt = configObj.write_wrapper_config2run_file(configName, line_cnt)
+            del configObj
+
+
     def finalize(self):
         self.runf.close()
         #writeJobFile(self.run_outname)
diff --git a/contrib/stack/topsStack/burstRampIon.py b/contrib/stack/topsStack/burstRampIon.py
new file mode 100755
index 00000000..d2c810fc
--- /dev/null
+++ b/contrib/stack/topsStack/burstRampIon.py
@@ -0,0 +1,169 @@
+#!/usr/bin/env python3
+
+# Author: Cunren Liang
+# Copyright 2022
+
+import os
+import copy
+import glob
+import shutil
+import argparse
+import numpy as np
+
+import isce
+import isceobj
+import s1a_isce_utils as ut
+from VRTManager import Swath
+from isceobj.TopsProc.runIon import computeDopplerOffset
+
+
+def createParser():
+    parser = argparse.ArgumentParser(description='compute burst phase ramps using azimuth ionospheric shift')
+    parser.add_argument('-k', '--reference_stack', type=str, dest='reference_stack', required=True,
+            help='Directory with the reference image of the stack')
+    parser.add_argument('-i', '--input', dest='input', type=str, required=True,
+            help='input azimuth ionospheric shift')
+    parser.add_argument('-o', '--output', dest='output', type=str, required=True,
+            help='Directory with output')
+    parser.add_argument('-r', '--nrlks', dest='nrlks', type=int, default=1, 
+            help='number of range looks of azimuth ionospheric shift. Default: 1')
+    parser.add_argument('-a', '--nalks', dest='nalks', type=int, default=1, 
+            help='number of azimuth looks of azimuth ionospheric shift. Default: 1')
+    parser.add_argument('-t', '--ion_height', dest='ion_height', type=float, default=200.0, 
+            help='height of ionospheric layer above the Earth surface in km. Default: 200.0')
+
+    return parser
+
+
+def cmdLineParse(iargs = None):
+    parser = createParser()
+    return parser.parse_args(args=iargs)
+
+
+def main(iargs=None):
+    '''
+    compute burst phase ramps using azimuth ionospheric shift
+    both regular and subband bursts are merged in a box defined by reference of stack
+
+    '''
+    inps = cmdLineParse(iargs)
+
+    #ionospheric layer height (m)
+    ionHeight = inps.ion_height * 1000.0
+    #earth's radius (m)
+    earthRadius = 6371 * 1000.0
+
+    img = isceobj.createImage()
+    img.load(inps.input + '.xml')
+    width = img.width
+    length = img.length
+    ionShift = np.fromfile(inps.input, dtype=np.float32).reshape(length, width)
+
+    swathList = ut.getSwathList(inps.reference_stack)
+    frameReferenceAll = [ut.loadProduct(os.path.join(inps.reference_stack, 'IW{0}.xml'.format(swath))) for swath in swathList]
+
+    refSwaths = [Swath(x) for x in frameReferenceAll]
+    topSwath = min(refSwaths, key = lambda x: x.sensingStart)
+    botSwath = max(refSwaths, key = lambda x: x.sensingStop)
+    leftSwath = min(refSwaths, key = lambda x: x.nearRange)
+    rightSwath = max(refSwaths, key = lambda x: x.farRange)
+
+    totalWidth  = int(np.round((rightSwath.farRange - leftSwath.nearRange)/leftSwath.dr + 1))
+    totalLength = int(np.round((botSwath.sensingStop - topSwath.sensingStart).total_seconds()/topSwath.dt + 1 ))
+
+    #!!!should CHECK if input ionospheric shift is really in this geometry
+    nearRange = leftSwath.nearRange
+    sensingStart = topSwath.sensingStart
+    dr = leftSwath.dr
+    dt = topSwath.dt
+
+    for swath in swathList:
+        frameReference = ut.loadProduct(os.path.join(inps.reference_stack, 'IW{0}.xml'.format(swath)))
+        swathObj = Swath(frameReference)
+
+        minBurst = frameReference.bursts[0].burstNumber
+        maxBurst = frameReference.bursts[-1].burstNumber
+        #!!!should CHECK if input ionospheric shift is really in this geometry
+        if minBurst==maxBurst:
+            print('Skipping processing of swath {0}'.format(swath))
+            continue
+
+        midBurstIndex = round((minBurst + maxBurst) / 2.0) - minBurst
+        midBurst = frameReference.bursts[midBurstIndex]
+
+        outputDir = os.path.join(inps.output, 'IW{0}'.format(swath))
+        os.makedirs(outputDir, exist_ok=True)
+
+        #compute mean offset: indices start from 0
+        firstLine   = int(np.round((swathObj.sensingStart - sensingStart).total_seconds()/(dt*inps.nalks))) + 1
+        lastLine    = int(np.round((swathObj.sensingStop - sensingStart).total_seconds()/(dt*inps.nalks))) - 1
+        firstColumn = int(np.round((swathObj.nearRange - nearRange)/(dr*inps.nrlks)))     + 1
+        lastColumn  = int(np.round((swathObj.farRange - nearRange)/(dr*inps.nrlks)))     - 1
+        ionShiftSwath = ionShift[firstLine:lastLine+1, firstColumn:lastColumn+1]
+        ionShiftSwathValid = ionShiftSwath[np.nonzero(ionShiftSwath!=0)]
+
+        if ionShiftSwathValid.size == 0:
+            ionShiftSwathMean = 0.0
+            print('mean azimuth ionospheric shift of swath {}: {} single look azimuth lines'.format(swath, ionShiftSwathMean))
+        else:
+            ionShiftSwathMean = np.mean(ionShiftSwathValid, dtype=np.float64)
+            print('mean azimuth ionospheric shift of swath {}: {} single look azimuth lines'.format(swath, ionShiftSwathMean))
+
+        for burst in frameReference.bursts:
+            #compute mean offset: indices start from 0
+            firstLine   = int(np.round((burst.sensingStart - sensingStart).total_seconds()/(dt*inps.nalks))) + 1
+            lastLine    = int(np.round((burst.sensingStop - sensingStart).total_seconds()/(dt*inps.nalks))) - 1
+            firstColumn = int(np.round((burst.startingRange - nearRange)/(dr*inps.nrlks))) + 1
+            lastColumn  = int(np.round((burst.startingRange + (burst.numberOfSamples - 1) * dr - nearRange)/(dr*inps.nrlks))) - 1
+            ionShiftBurst = ionShift[firstLine:lastLine+1, firstColumn:lastColumn+1]
+            
+            ionShiftBurstValid = ionShiftBurst[np.nonzero(ionShiftBurst!=0)]
+            if ionShiftBurstValid.size < (lastLine - firstLine + 1) * (lastColumn - firstColumn + 1) / 2.0:
+                ionShiftBurstMean = 0.0
+                print('mean azimuth ionospheric shift of burst {}: 0.0 single look azimuth lines'.format(burst.burstNumber))
+            else:
+                #ionShiftBurstMean should use both phaseRamp1 and phaseRamp2, while
+                #ionShiftSwathMean should use phaseRamp2 only as in (to be consistent with) previous ESD
+                #The above is tested against runIon.py in topsApp.py
+                #ionShiftBurstMean = np.mean(ionShiftBurstValid, dtype=np.float64) - ionShiftSwathMean
+                ionShiftBurstMean = np.mean(ionShiftBurstValid, dtype=np.float64)
+                print('mean azimuth ionospheric shift of burst {}: {} single look azimuth lines'.format(burst.burstNumber, ionShiftBurstMean))
+
+            #compute burst phase ramps
+            (dopplerOffset, Ka) = computeDopplerOffset(burst, 1, burst.numberOfLines, 1, burst.numberOfSamples, nrlks=1, nalks=1)
+
+            #satellite height
+            satHeight = np.linalg.norm(burst.orbit.interpolateOrbit(burst.sensingMid, method='hermite').getPosition())
+            #orgininal doppler offset should be multiplied by this ratio
+            ratio = ionHeight/(satHeight-earthRadius)
+
+            #phaseRamp1 and phaseRamp2 have same sign according to Liang et. al., 2019.
+            #phase ramp due to imaging geometry
+            phaseRamp1 = dopplerOffset * burst.azimuthTimeInterval * ratio * (burst.azimuthTimeInterval * Ka[None,:] * 4.0 * np.pi)
+            #phase ramp due to non-zero doppler centroid frequency, ESD
+            phaseRamp2 = dopplerOffset * burst.azimuthTimeInterval * Ka[None,:] * 2.0 * np.pi * burst.azimuthTimeInterval
+            phaseRamp = (phaseRamp1 + phaseRamp2) * ionShiftBurstMean - phaseRamp2 * ionShiftSwathMean
+
+            outfile = os.path.join(outputDir, '%s_%02d.float'%('burst', burst.burstNumber))
+            
+            phaseRamp.astype(np.float32).tofile(outfile)
+
+            image = isceobj.createImage()
+            image.setDataType('FLOAT')
+            image.setFilename(outfile)
+            image.extraFilename = outfile + '.vrt'
+            image.setWidth(burst.numberOfSamples)
+            image.setLength(burst.numberOfLines)
+            image.renderHdr()
+
+
+
+if __name__ == '__main__':
+    '''
+    Main driver.
+    '''
+    # Main Driver
+    main()
+
+
+
diff --git a/contrib/stack/topsStack/filtIonShift.py b/contrib/stack/topsStack/filtIonShift.py
new file mode 100755
index 00000000..3bec81c9
--- /dev/null
+++ b/contrib/stack/topsStack/filtIonShift.py
@@ -0,0 +1,242 @@
+#!/usr/bin/env python3
+
+# Author: Cunren Liang
+# Copyright 2022
+
+import os
+import copy
+import glob
+import shutil
+import argparse
+import numpy as np
+
+import isce
+import isceobj
+from isceobj.TopsProc.runIon import adaptive_gaussian
+from isceobj.TopsProc.runIon import weight_fitting
+from isceobj.TopsProc.runIon import fit_surface
+from isceobj.TopsProc.runIon import cal_surface
+import s1a_isce_utils as ut
+
+
+def createParser():
+    parser = argparse.ArgumentParser(description='compute and filter azimuth ionospheric shift [unit: masBurst.azimuthTimeInterval]')
+    parser.add_argument('-k', '--reference_stack', type=str, dest='reference_stack', required=True,
+                        help='Directory with the reference image of the stack')
+    parser.add_argument('-f', '--reference', type=str, dest='reference', required=True,
+                        help='Directory with the reference image (coregistered or reference of the stack)')
+    parser.add_argument('-s', '--secondary', type=str, dest='secondary', required=True,
+                        help='Directory with the secondary image (coregistered or reference of the stack)')
+    parser.add_argument('-i', '--input', dest='input', type=str, required=True,
+            help='input ionosphere')
+    parser.add_argument('-c', '--coherence', dest='coherence', type=str, required=True,
+            help='coherence, e.g. raw_no_projection.cor')
+    parser.add_argument('-o', '--output', dest='output', type=str, required=True,
+            help='output azimuth ionospheric shift')
+    parser.add_argument('-b', '--win_min', dest='win_min', type=int, default=100, 
+            help='minimum filtering window size')
+    parser.add_argument('-d', '--win_max', dest='win_max', type=int, default=200, 
+            help='maximum filtering window size')
+    parser.add_argument('-r', '--nrlks', dest='nrlks', type=int, default=1, 
+            help='number of range looks. Default: 1')
+    parser.add_argument('-a', '--nalks', dest='nalks', type=int, default=1, 
+            help='number of azimuth looks. Default: 1')
+    #parser.add_argument('-m', '--masked_areas', dest='masked_areas', type=int, nargs='+', action='append', default=None,
+    #        help='This is a 2-d list. Each element in the 2-D list is a four-element list: [firstLine, lastLine, firstColumn, lastColumn], with line/column numbers starting with 1. If one of the four elements is specified with -1, the program will use firstLine/lastLine/firstColumn/lastColumn instead. e.g. two areas masked out: --masked_areas 10 20 10 20 --masked_areas 110 120 110 120')
+
+    return parser
+
+
+def cmdLineParse(iargs = None):
+    parser = createParser()
+    return parser.parse_args(args=iargs)
+
+
+def main(iargs=None):
+    '''
+    check overlap among all acquistions, only keep the bursts that in the common overlap,
+    and then renumber the bursts.
+    '''
+    inps = cmdLineParse(iargs)
+
+    '''
+    calculate azimuth shift caused by ionosphere using ionospheric phase
+    '''
+
+    #################################################
+    #SET PARAMETERS HERE
+    #gaussian filtering window size
+    #size = np.int(np.around(width / 12.0))
+    #size = ionParam.ionshiftFilteringWinsize
+    size_max = inps.win_max
+    size_min = inps.win_min
+
+    #THESE SHOULD BE GOOD ENOUGH, NO NEED TO SET IN setup(self)
+    #if applying polynomial fitting
+    #0: no fitting, 1: with fitting
+    fit = 0
+    corThresholdIonshift = 0.85
+    #################################################
+
+
+####################################################################
+    #STEP 1. GET DERIVATIVE OF IONOSPHERE
+####################################################################
+
+    #get files
+    ionfile = inps.input
+    #we are using filtered ionosphere, so we should use coherence file that is not projected.
+    #corfile = os.path.join(ionParam.ionDirname, ionParam.ioncalDirname, ionParam.ionCor)
+    corfile = inps.coherence
+    img = isceobj.createImage()
+    img.load(ionfile + '.xml')
+    width = img.width
+    length = img.length
+    amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
+    ion = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
+    cor = (np.fromfile(corfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
+
+    ########################################################################################
+    #AFTER COHERENCE IS RESAMPLED AT grd2ion, THERE ARE SOME WIRED VALUES
+    cor[np.nonzero(cor<0)] = 0.0
+    cor[np.nonzero(cor>1)] = 0.0
+    ########################################################################################
+
+    #get the azimuth derivative of ionosphere
+    dion = np.diff(ion, axis=0)
+    dion = np.concatenate((dion, np.zeros((1,width))), axis=0)
+
+    #remove the samples affected by zeros
+    flag_ion0 = (ion!=0)
+    #moving down by one line
+    flag_ion1 = np.roll(flag_ion0, 1, axis=0)
+    flag_ion1[0,:] = 0
+    #moving up by one line
+    flag_ion2 = np.roll(flag_ion0, -1, axis=0)
+    flag_ion2[-1,:] = 0
+    #now remove the samples affected by zeros
+    flag_ion = flag_ion0 * flag_ion1 * flag_ion2
+    dion *= flag_ion
+
+    flag = flag_ion * (cor>corThresholdIonshift)
+    index = np.nonzero(flag)
+
+
+####################################################################
+    #STEP 2. FIT A POLYNOMIAL TO THE DERIVATIVE OF IONOSPHERE
+####################################################################
+
+    order = 3
+
+    #look for data to use
+    point_index = np.nonzero(flag)
+    m = point_index[0].shape[0]
+
+    #calculate input index matrix
+    x0=np.matlib.repmat(np.arange(width), length, 1)
+    y0=np.matlib.repmat(np.arange(length).reshape(length, 1), 1, width)
+
+    x = x0[point_index].reshape(m, 1)
+    y = y0[point_index].reshape(m, 1)
+    z = dion[point_index].reshape(m, 1)
+    w = cor[point_index].reshape(m, 1)
+
+    #convert to higher precision type before use
+    x=np.asfarray(x,np.float64)
+    y=np.asfarray(y,np.float64)
+    z=np.asfarray(z,np.float64)
+    w=np.asfarray(w,np.float64)
+    coeff = fit_surface(x, y, z, w, order)
+
+    rgindex = np.arange(width)
+    azindex = np.arange(length).reshape(length, 1)
+    #convert to higher precision type before use
+    rgindex=np.asfarray(rgindex,np.float64)
+    azindex=np.asfarray(azindex,np.float64)
+    dion_fit = cal_surface(rgindex, azindex, coeff, order)
+
+    #no fitting
+    if fit == 0:
+        dion_fit *= 0
+    dion_res = (dion - dion_fit)*(dion!=0)
+
+
+####################################################################
+    #STEP 3. FILTER THE RESIDUAL OF THE DERIVATIVE OF IONOSPHERE
+####################################################################
+
+    #this will be affected by low coherence areas like water, so not use this.
+    #filter the derivation of ionosphere
+    #if size % 2 == 0:
+    #    size += 1
+    #sigma = size / 2.0
+
+    #g2d = gaussian(size, sigma, scale=1.0)
+    #scale = ss.fftconvolve((dion_res!=0), g2d, mode='same')
+    #dion_filt = ss.fftconvolve(dion_res, g2d, mode='same') / (scale + (scale==0))
+
+    #minimize the effect of low coherence pixels
+    cor[np.nonzero( (cor<0.85)*(cor!=0) )] = 0.00001
+    dion_filt = adaptive_gaussian(dion_res, cor, size_max, size_min)
+
+    dion = (dion_fit + dion_filt)*(dion!=0)
+
+    #return dion
+
+
+####################################################################
+    #STEP 4. CONVERT TO AZIMUTH SHIFT
+####################################################################
+
+    #!!! use number of swaths in reference for now
+
+    #these are coregistered secondaries, so there should be <= swaths in reference of entire stack
+    referenceSwathList = ut.getSwathList(inps.reference)
+    secondarySwathList = ut.getSwathList(inps.secondary)
+    swathList = list(sorted(set(referenceSwathList+secondarySwathList)))
+
+    #swathList = ut.getSwathList(inps.reference)
+    firstSwath = None
+    for swath in swathList:
+        frameReference = ut.loadProduct(os.path.join(inps.reference_stack, 'IW{0}.xml'.format(swath)))
+
+        minBurst = frameReference.bursts[0].burstNumber
+        maxBurst = frameReference.bursts[-1].burstNumber
+        if minBurst==maxBurst:
+            print('Skipping processing of swath {0}'.format(swath))
+            continue
+        else:
+            firstSwath = swath
+            break
+    if firstSwath is None:
+        raise Exception('no valid swaths!')
+
+    midBurstIndex = round((minBurst + maxBurst) / 2.0) - minBurst
+    masBurst = frameReference.bursts[midBurstIndex]
+
+    #shift casued by ionosphere [unit: masBurst.azimuthTimeInterval]
+    rng = masBurst.rangePixelSize * ((np.arange(width))*inps.nrlks + (inps.nrlks - 1.0) / 2.0) + masBurst.startingRange
+    Ka = masBurst.azimuthFMRate(rng)
+    ionShift = dion / (masBurst.azimuthTimeInterval * inps.nalks) / (4.0 * np.pi) / Ka[None, :] / masBurst.azimuthTimeInterval
+
+    #output
+    outfile = inps.output
+    os.makedirs(os.path.dirname(outfile), exist_ok=True)
+    ion = np.zeros((length*2, width), dtype=np.float32)
+    ion[0:length*2:2, :] = amp
+    ion[1:length*2:2, :] = ionShift
+    ion.astype(np.float32).tofile(outfile)
+    img.filename = outfile
+    img.extraFilename = outfile + '.vrt'
+    img.renderHdr()
+
+
+if __name__ == '__main__':
+    '''
+    Main driver.
+    '''
+    # Main Driver
+    main()
+
+
+
diff --git a/contrib/stack/topsStack/invertIon.py b/contrib/stack/topsStack/invertIon.py
index b3e5f1e9..43abd5f1 100755
--- a/contrib/stack/topsStack/invertIon.py
+++ b/contrib/stack/topsStack/invertIon.py
@@ -96,6 +96,8 @@ def cmdLineParse():
     parser = argparse.ArgumentParser(description='least squares estimation')
     parser.add_argument('--idir', dest='idir', type=str, required=True,
             help = 'input directory where each pair (YYYYMMDD_YYYYMMDD) is located. only folders are recognized')
+    parser.add_argument('--filename', dest='filename', type=str, default='filt.ion',
+            help = 'input file name. default: filt.ion')
     parser.add_argument('--odir', dest='odir', type=str, required=True,
             help = 'output directory for estimated result of each date')
     parser.add_argument('--zro_date', dest='zro_date', type=str, default=None,
@@ -205,7 +207,7 @@ def cmdLineParse():
     ndate = len(dates)
     npair = len(pairs)
 
-    ionfile = os.path.join(idir, pairs[0], 'ion_cal', 'filt.ion')
+    ionfile = os.path.join(idir, pairs[0], 'ion_cal', inps.filename)
 
     img = isceobj.createImage()
     img.load(ionfile+'.xml')
@@ -219,7 +221,7 @@ def cmdLineParse():
     wls = False
     stdPairs = np.ones((npair, length, width), dtype=np.float32)
     for i in range(npair):
-        ionfile = os.path.join(idir, pairs[i], 'ion_cal', 'filt.ion')
+        ionfile = os.path.join(idir, pairs[i], 'ion_cal', inps.filename)
         ionPairs[i, :, :] = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[1:length*2:2, :]
         #flag of valid/invalid is defined by amplitde image
         amp = (np.fromfile(ionfile, dtype=np.float32).reshape(length*2, width))[0:length*2:2, :]
diff --git a/contrib/stack/topsStack/ion_param.txt b/contrib/stack/topsStack/ion_param.txt
index 3896804d..a01bbbbf 100644
--- a/contrib/stack/topsStack/ion_param.txt
+++ b/contrib/stack/topsStack/ion_param.txt
@@ -2,6 +2,9 @@
 ###the values below are the default values used by the module
 ###remove # to set the parameters
 
+#consider burst properties in ionosphere computation: False
+#height of ionosphere layer in km: 200.0
+
 #maximum window size for filtering ionosphere phase:         200
 #minimum window size for filtering ionosphere phase:         100
 #maximum window size for filtering ionosphere azimuth shift: 150
diff --git a/contrib/stack/topsStack/s1_select_ion.py b/contrib/stack/topsStack/s1_select_ion.py
index ec517d1a..75d70551 100755
--- a/contrib/stack/topsStack/s1_select_ion.py
+++ b/contrib/stack/topsStack/s1_select_ion.py
@@ -422,8 +422,8 @@ def cmdLineParse():
     parser = argparse.ArgumentParser( description='select Sentinel-1A/B acquistions good for ionosphere correction. not used slices are moved to folder: not_used')
     parser.add_argument('-dir', dest='dir', type=str, required=True,
             help = 'directory containing the "S1*_IW_SLC_*.zip" files')
-    parser.add_argument('-sn', dest='sn', type=str, required=True,
-            help='south/north bound of area of interest, format: south/north')
+    parser.add_argument('-sn', dest='sn', type=float, nargs=2, required=True,
+            help='south and north bound of area of interest, format: -sn south north')
     parser.add_argument('-nr', dest='nr', type=int, default=10,
             help = 'minimum number of acquisitions for same starting ranges. default: 10')
 
@@ -438,7 +438,8 @@ def cmdLineParse():
 if __name__ == '__main__':
 
     inps = cmdLineParse()
-    s,n=[float(x) for x in inps.sn.split('/')]
+    #s,n=[float(x) for x in inps.sn.split('/')]
+    s,n=inps.sn
 
     #group the slices
     group = get_group(inps.dir)
diff --git a/contrib/stack/topsStack/stackSentinel.py b/contrib/stack/topsStack/stackSentinel.py
index ef92360e..24639d18 100755
--- a/contrib/stack/topsStack/stackSentinel.py
+++ b/contrib/stack/topsStack/stackSentinel.py
@@ -14,7 +14,7 @@
 import isce
 import isceobj
 from isceobj.Sensor.TOPS.Sentinel1 import Sentinel1
-from topsStack.Stack import config, run, sentinelSLC
+from topsStack.Stack import config, run, sentinelSLC, ionParamUsr
 
 
 helpstr = """
@@ -607,8 +607,9 @@ def checkCurrentStatusIonosphere(inps):
     acquisitionDates, stackReferenceDate, secondaryDates, safe_dict = get_dates(inps)
 
     pairs_same_starting_ranges, pairs_diff_starting_ranges = selectNeighborPairsIonosphere(safe_dict, inps.num_connections_ion)
+    dateListIon = getDatesIonosphere(pairs_same_starting_ranges, pairs_diff_starting_ranges)
     pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update = excludeExistingPairsIonosphere(pairs_same_starting_ranges, pairs_diff_starting_ranges, inps.work_dir)
-    dateListIon = getDatesIonosphere(pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update)
+    dateListIonUpdate = getDatesIonosphere(pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update)
 
     #report pairs of different swath starting ranges.
     pdiff = 'ionosphere phase estimation pairs with different swath starting ranges\n'
@@ -623,7 +624,7 @@ def checkCurrentStatusIonosphere(inps):
     with open('pairs_diff_starting_ranges.txt', 'w') as f:
         f.write(pdiff)
 
-    return dateListIon, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict
+    return dateListIonUpdate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, dateListIon
 
 
 ########################################
@@ -821,12 +822,12 @@ def offsetStack(inps, acquisitionDates, stackReferenceDate, secondaryDates, safe
     return i
 
 
-def ionosphereStack(inps, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i):
+def ionosphereStack(inps, dateListIonUpdate, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i):
 
     i+=1
     runObj = run()
     runObj.configure(inps, 'run_{:02d}_subband_and_resamp'.format(i))
-    runObj.subband_and_resamp(dateListIon, stackReferenceDate)
+    runObj.subband_and_resamp(dateListIonUpdate, stackReferenceDate)
     runObj.finalize()
 
     i+=1
@@ -871,6 +872,38 @@ def ionosphereStack(inps, dateListIon, stackReferenceDate, pairs_same_starting_r
     runObj.invertIon()
     runObj.finalize()
 
+
+    ionParamUsrObj = ionParamUsr(inps.param_ion)
+    ionParamUsrObj.configure()
+    if ionParamUsrObj.ION_considerBurstProperties:
+        i+=1
+        runObj = run()
+        runObj.configure(inps, 'run_{:02d}_filtIonShift'.format(i))
+        runObj.filtIonShift(pairs_same_starting_ranges_update + pairs_diff_starting_ranges_update)
+        runObj.finalize()
+
+        #From now on, need to reprocess all dates/pairs, not only newly added
+        i+=1
+        runObj = run()
+        runObj.configure(inps, 'run_{:02d}_invertIonShift'.format(i))
+        runObj.invertIonShift()
+        runObj.finalize()
+
+        i+=1
+        runObj = run()
+        runObj.configure(inps, 'run_{:02d}_burstRampIon'.format(i))
+        runObj.burstRampIon(dateListIon)
+        runObj.finalize()
+
+        i+=1
+        runObj = run()
+        runObj.configure(inps, 'run_{:02d}_mergeBurstRampIon'.format(i))
+        runObj.mergeBurstRampIon(dateListIon, stackReferenceDate)
+        runObj.finalize()
+
+    del ionParamUsrObj
+        
+
     return i
 
 
@@ -1015,8 +1048,19 @@ def main(iargs=None):
     elif not os.path.isfile(inps.param_ion):
         print("Ion parameter file is missing. Ionospheric estimation will not be done.")
     else:
-        dateListIon, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict = checkCurrentStatusIonosphere(inps)
-        i = ionosphereStack(inps, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i)
+        #############################################
+        #check if consider burst properties
+        ionParamUsrObj = ionParamUsr(inps.param_ion)
+        ionParamUsrObj.configure()
+        if ionParamUsrObj.ION_considerBurstProperties:
+            if inps.coregistration != 'NESD':
+                raise Exception("Coregistration options '-C', '--coregistration' must be NESD when 'consider burst properties in ionosphere computation' is True")
+        del ionParamUsrObj
+        #############################################
+
+        dateListIonUpdate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, dateListIon = checkCurrentStatusIonosphere(inps)
+        i = ionosphereStack(inps, dateListIonUpdate, dateListIon, stackReferenceDate, pairs_same_starting_ranges_update, pairs_diff_starting_ranges_update, safe_dict, i)
+
 
 if __name__ == "__main__":
     # Main engine