diff --git a/Spot_polynyas.py b/Spot_polynyas.py
new file mode 100644
index 0000000000000000000000000000000000000000..c4ec8e4c7f866f0c4f45cb7fa9c66a8d0873e599
--- /dev/null
+++ b/Spot_polynyas.py
@@ -0,0 +1,393 @@
+# This script takes Sea Ice data, and spot polynyas with a fixed threshold
+# Authors: C. Pelletier edited by N. Pirlet
+# Cite this paper for credits:
+
+import numpy as np
+import netCDF4 as nc
+import os
+import sys
+import time
+
+# Adapt------------
+year_start = 2001
+year_stop = 2017
+
+siconc_thres = 1.00
+sithic_thres = 0.25
+
+lfi = False # If you want polynyas close to lfi to be coastal polynyas (True) (NO_LFI => False)
+
+is_cyclic = True
+if(is_cyclic):
+ n_ovlp = 1 # West-East bondary overlap
+else:
+ n_ovlp = 0
+
+
+# Contient and area files
+input_geometry = "/cofast/npirlet/eANT_025/eANT025_mesh_mask_v2.nc"
+input_cell_area = "/cofast/npirlet/eANT_025/eANT025_cell_area.nc"
+
+#--------------------
+
+fillv = -10
+open_oce_value = 0
+iced_value = 1
+antarctica_pol_value = 2
+island_pol_value = 3
+ocean_pol_value = 4
+
+
+
+fillv_float = -1.e20
+
+
+
+geom = nc.Dataset(input_geometry, mode='r')
+if(is_cyclic):
+ # Creation of "land" mask
+ tmaskutil = geom.variables['tmaskutil'][0,0:370,:-n_ovlp] == 0
+ isf_draft = geom.variables['isf_draft'][0,0:370,:-n_ovlp] > 0
+ msk_cont = np.full(np.shape(tmaskutil),False)
+ msk_cont[tmaskutil == True] = True
+ msk_cont[isf_draft == True] = True
+ is_cont = msk_cont
+ is_ant = msk_cont
+ print(np.shape(msk_cont))
+
+ if lfi == True:
+ mask_LFI = np.load("/cofast/npirlet/eANT_025/mask_LFI_MAMJJASO_2001-2017.npy")
+ mask_LFI = mask_LFI[:,:-n_ovlp]
+ msk_cont[mask_LFI >= 50] = True
+else:
+ # Creation of "land" mask
+ tmaskutil = geom.variables['tmaskutil'][0,0:370,:] == 0
+ isf_draft = geom.variables['isf_draft'][0,0:370,:] > 0
+ msk_cont = np.full(np.shape(tmaskutil),False)
+ msk_cont[tmaskutil == True] = True
+ msk_cont[isf_draft == True] = True
+ is_cont = msk_cont
+ is_ant = msk_cont
+ print(np.shape(msk_cont))
+
+ if lfi == True:
+ mask_LFI = np.load("/cofast/npirlet/eANT_025/mask_LFI_MAMJJASO_2001-2017.npy")
+ msk_cont[mask_LFI >= 50] = True
+
+lat = geom.variables['nav_lat'][0:370,:]
+lon = geom.variables['nav_lon'][0:370,:]
+geom.close()
+
+areac = nc.Dataset(input_cell_area, mode='r')
+if(is_cyclic):
+ cell_area = areac.variables['areacello'][0,0:370,:-n_ovlp]
+ print(np.shape(cell_area))
+else:
+ cell_area = areac.variables['areacello'][0,0:370,:]
+ print(np.shape(cell_area))
+areac.close()
+
+ny, nx = np.shape(cell_area)
+
+
+is_border = np.zeros(shape=[ny, nx], dtype=bool)
+
+is_border[:1,:] = True
+is_border[-1:,:] = True
+if(not is_cyclic):
+ is_border[:,:1] = True
+ is_border[:,-1:] = True
+
+
+
+is_cont_neigh = np.ndarray(shape=[ny, nx, 4], dtype=bool)
+is_cont_neigh[:,:,0] = np.roll(is_cont, shift=+1, axis=-2)
+is_cont_neigh[:,:,1] = np.roll(is_cont, shift=-1, axis=-2)
+is_cont_neigh[:,:,2] = np.roll(is_cont, shift=+1, axis=-1)
+is_cont_neigh[:,:,3] = np.roll(is_cont, shift=-1, axis=-1)
+
+is_ant_neigh = np.ndarray(shape=[ny, nx, 4], dtype=bool)
+is_ant_neigh[:,:,0] = np.roll(is_ant, shift=+1, axis=-2)
+is_ant_neigh[:,:,1] = np.roll(is_ant, shift=-1, axis=-2)
+is_ant_neigh[:,:,2] = np.roll(is_ant, shift=+1, axis=-1)
+is_ant_neigh[:,:,3] = np.roll(is_ant, shift=-1, axis=-1)
+
+
+
+old_age_polynya = np.zeros(shape=[ny, nx], dtype=int)
+
+
+for year in range(year_start,year_stop+1):
+
+ # Data path to adapt
+ input_seaice_conc = "/cofast/npirlet/LUCIA/NO_LFI_40/NO_LFI_40_1d_siconc_y{}.nc".format(year)
+ input_seaice_thick = "/cofast/npirlet/LUCIA/NO_LFI_40/NO_LFI_40_1d_sivolu_y{}.nc".format(year)
+ outfile = "/cofast/npirlet/LUCIA/NO_LFI_40/polynyas/polynyas_NO_LFI_40_y{}_30cm.nc".format(year)
+ #-------------------
+
+ # Load sea ice concentration data
+ tmpdat = nc.Dataset(input_seaice_conc, mode='r')
+
+ if(is_cyclic):
+ siconc = tmpdat.variables['siconc'][:,0:370,:-n_ovlp]
+ else:
+ siconc = tmpdat.variables['siconc'][:,0:370,:]
+
+ tmpdat.close()
+
+
+ # Load sea ice thickness data
+ tmpdat = nc.Dataset(input_seaice_thick, mode='r')
+
+ if(is_cyclic):
+ sithic = tmpdat.variables['sivolu'][:,0:370,:-n_ovlp]
+ else:
+ sithic = tmpdat.variables['sivolu'][:,0:370,:]
+
+ time_var = tmpdat.variables['time_counter'][:]
+ long_name_time = tmpdat.variables['time_counter'].standard_name
+ units_time = tmpdat.variables['time_counter'].units
+ calendar_time = tmpdat.variables['time_counter'].calendar
+ tmpdat.close()
+
+
+ # Begin process
+
+ nt = (np.shape(sithic))[0]
+
+ idx_open_patch = fillv * np.ones(shape=[nt, ny, nx], dtype=int)
+
+ area_polynya = np.zeros(shape=[nt, ny, nx], dtype=float)
+
+ age_polynya = np.zeros(shape=[nt, ny, nx], dtype=int)
+
+ idx_polynya = np.zeros(shape=[nt,ny,nx], dtype=int)
+
+
+
+
+ for t in range(0,nt):
+
+ cnt_pol = 0
+ cnt_oce_pol = 0
+ cnt_island_pol = 0
+ cnt_ant_pol = 0
+
+ unsorted_idx_polynya = np.zeros(shape=[ny, nx], dtype=int)
+ all_poly_size = []
+
+ t_start = time.time()
+
+ is_open = np.logical_and(np.logical_not(is_cont), siconc[t,:,:] < siconc_thres, sithic[t,:,:] < sithic_thres)
+ idx_open_patch[t,:,:] = np.where(np.logical_and(np.logical_not(is_open), np.logical_not(is_cont)), iced_value, idx_open_patch[t,:,:])
+
+
+
+ is_open_neigh = np.ndarray(shape=[ny,nx,4], dtype=bool)
+
+ is_open_neigh[:,:,0] = np.roll(is_open, shift=+1, axis=-2)
+ is_open_neigh[:,:,1] = np.roll(is_open, shift=-1, axis=-2)
+ is_open_neigh[:,:,2] = np.roll(is_open, shift=+1, axis=-1)
+ is_open_neigh[:,:,3] = np.roll(is_open, shift=-1, axis=-1)
+
+
+ explored = np.zeros(shape=[ny, nx], dtype=bool)
+ candidates = np.transpose(np.nonzero(np.logical_and(is_open, np.logical_not(explored))))
+
+
+ curr_patch_idx = 0
+
+ sizes_patch = []
+ touched_cont = []
+
+
+ while(np.size(candidates) > 0):
+
+ curr_patch_idx+=1
+
+ is_in_currpatch = np.zeros(shape=[ny, nx], dtype=bool)
+ is_in_currpatch[candidates[0,0], candidates[0,1]] = True
+
+ keep_going = True
+ cnt_explore = 0
+
+ while(keep_going):
+
+ n_cp = np.count_nonzero(is_in_currpatch)
+
+ # add x+ neighbors
+ is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,0], is_in_currpatch), shift=-1, axis=0))
+
+ # add x- neighbors
+ is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,1], is_in_currpatch), shift=1, axis=0))
+
+ # add y+ neighbors
+ is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,2], is_in_currpatch), shift=-1, axis=1))
+
+ # add y- neighbors
+ is_in_currpatch = np.logical_or(is_in_currpatch, np.roll(np.logical_and(is_open_neigh[:,:,3], is_in_currpatch), shift=1, axis=1))
+
+
+ n_old = n_cp
+ n_cp = np.count_nonzero(is_in_currpatch)
+
+ n_new = n_cp - n_old
+ keep_going = (n_new > 0)
+ cnt_explore+=1
+
+
+ if( np.any( np.logical_and(is_in_currpatch, is_border))):
+ # patch touching border: open ocean
+
+ idx_open_patch[t,:,:] = np.where( is_in_currpatch, open_oce_value, idx_open_patch[t,:,:] )
+
+ else:
+ # patch doesn't touch border: polynya
+
+ cnt_pol+=1
+
+ curr_area_pol = np.sum(np.ma.array(data = cell_area, mask = np.logical_not(is_in_currpatch)))
+
+ all_poly_size.append(curr_area_pol)
+
+
+ area_polynya[t,:,:] = np.where(is_in_currpatch, curr_area_pol, area_polynya[t,:,:])
+
+ unsorted_idx_polynya = np.where(is_in_currpatch, cnt_pol, unsorted_idx_polynya)
+
+
+ # polynya age
+
+ if(np.any( np.logical_and( old_age_polynya > 0, is_in_currpatch ) ) ):
+ # polynya previously existed: its age is the maximum age of any of the current polynya node taken from the old time step
+ prev_age_tmp = np.ma.array(data = old_age_polynya, mask = np.logical_not(is_in_currpatch))
+
+ age_polynya[t,:,:] = np.where(is_in_currpatch, np.ma.amax(prev_age_tmp) + 1, age_polynya[t,:,:])
+
+ else:
+ # new polynya: aged 1
+ age_polynya[t,:,:] = np.where(is_in_currpatch, 1, age_polynya[t,:,:])
+
+
+
+ tmp_mask = is_in_currpatch[:,:,None] * (np.ones(shape=[4], dtype=bool))[None,None,:]
+ check_cont = np.logical_and(is_cont_neigh, tmp_mask)
+
+ if(np.any(check_cont)):
+ # coastal polynya: Antarctica or island
+
+
+ check_ant = np.logical_and(is_ant_neigh, tmp_mask)
+ if(np.any(check_ant)):
+ # polynya touches Antarctica
+ idx_open_patch[t,:,:] = np.where(is_in_currpatch, antarctica_pol_value, idx_open_patch[t,:,:])
+
+
+ cnt_ant_pol+=1
+
+ else:
+ # insular (non-Antarctica) polynya
+ idx_open_patch[t,:,:] = np.where(is_in_currpatch, island_pol_value, idx_open_patch[t,:,:])
+
+
+ cnt_island_pol+=1
+
+
+ else: # np.any(check_cont))
+ # the polynya doesn't touch masked cell: it's an ocean polynya
+
+ idx_open_patch[t,:,:] = np.where(is_in_currpatch, ocean_pol_value, idx_open_patch[t,:,:])
+
+ cnt_oce_pol+=1
+
+
+ sizes_patch.append(n_cp)
+
+ # we have explored the current patch
+ explored = np.where(is_in_currpatch, True, explored)
+
+ # the remaining polynya candidates are the ice-free cells which have not been explored
+ candidates = np.transpose(np.nonzero(np.logical_and(is_open, np.logical_not(explored))))
+
+ # age: fillvalue over continent, 0 over free-ocean or sea-ice cover
+ age_polynya[t,:,:] = np.where(is_cont, fillv, np.where( np.logical_or(idx_open_patch[t,:,:] == open_oce_value, idx_open_patch[t,:,:] == iced_value),0, age_polynya[t,:,:]))
+
+ # area polynya: fillvalue over continent
+ area_polynya[t,:,:] = np.where(is_cont, fillv_float, area_polynya[t,:,:])
+
+ # sort the polynya index by decreasing polynya extent
+ idx_polynya[t,:,:] = unsorted_idx_polynya
+ argsort = (np.argsort(all_poly_size))[::-1]
+
+ for i in range(0,cnt_pol):
+ idx_polynya[t,:,:] = np.where(unsorted_idx_polynya == argsort[i]+1, i + 1, idx_polynya[t,:,:])
+
+ idx_polynya[t,:,:] = np.where(is_cont, fillv, idx_polynya[t,:,:])
+
+ walltime_ts = time.time() - t_start
+ print('Time step %d'%t+': spotted %d'%cnt_pol+' polynya: %d'%cnt_ant_pol+' attached to Antarctica, %d'%cnt_island_pol+' attached to an island and %d'%cnt_oce_pol+' oceanic (walltime %.1f'%walltime_ts+' sec.)')
+
+
+
+ out_nc = nc.Dataset(outfile, mode='w')
+
+ out_nc.createDimension('time_counter', 0)
+ out_nc.createDimension('axis_nbounds', 2)
+ out_nc.createDimension('y', ny)
+ out_nc.createDimension('x', nx+n_ovlp)
+
+
+ tmp = out_nc.createVariable(varname="time_counter", datatype='f8', dimensions=['time_counter'])
+ tmp.standard_name = long_name_time
+ tmp.long_name = long_name_time
+ tmp.units = units_time
+ tmp.calendar = calendar_time
+ tmp[:] = time_var[:]
+
+ tmp = out_nc.createVariable(varname="latitude", datatype='f8', dimensions=['y', 'x'])
+ tmp.standard_name = "latitude"
+ tmp.units = "degrees North"
+ tmp[:] = lat
+
+ tmp = out_nc.createVariable(varname="longitude", datatype='f8', dimensions=['y', 'x'])
+ tmp.standard_name = "longitude"
+ tmp.units = "degrees East"
+ tmp[:] = lon
+
+
+ tmp = out_nc.createVariable(varname='polynya_status', datatype='f8', dimensions=['time_counter', 'y', 'x'], fill_value=fillv)
+ tmp.coordinates = 'time_counter latitude longitude'
+ tmp.description = "polnynya status. masked = continent; %d"%open_oce_value+' = open ocean; %d'%iced_value+" = sea-ice cover; %d"%antarctica_pol_value+" = polynya attached to Antarctica; %d"%island_pol_value+" = polynya attached to an island; %d"%ocean_pol_value+" = oceanic polynya"
+ if(is_cyclic):
+ tmp[:] = np.dstack((idx_open_patch, idx_open_patch[:,:,:n_ovlp]))
+ else:
+ tmp[:] = idx_open_patch
+
+
+ tmp = out_nc.createVariable(varname='area_polynya', datatype='f4', dimensions=['time_counter', 'y', 'x'], fill_value = fillv_float)
+ if(is_cyclic):
+ tmp[:] = np.dstack((area_polynya, area_polynya[:,:,:n_ovlp]))
+ else:
+ tmp[:] = area_polynya
+ tmp.coordinates = 'time_counter latitude longitude'
+ tmp.units = 'm2'
+
+
+ tmp = out_nc.createVariable(varname='age_polynya', datatype='i2', dimensions=['time_counter', 'y', 'x'], fill_value = fillv)
+ tmp.coordinates = 'time_counter latitude longitude'
+ tmp.units = 'time step (days)'
+ if(is_cyclic):
+ tmp[:] = np.dstack((age_polynya, age_polynya[:,:,:n_ovlp]))
+ else:
+ tmp[:] = age_polynya
+
+
+ tmp = out_nc.createVariable(varname='idx_polynya', datatype='i2', dimensions=['time_counter', 'y', 'x'], fill_value = fillv)
+ tmp.coordinates = 'time_counter latitude longitude'
+ tmp.description = 'index counting polynyas (in decreasing order of extent)'
+ if(is_cyclic):
+ tmp[:] = np.dstack((idx_polynya, idx_polynya[:,:,:n_ovlp]))
+ else:
+ tmp[:] = idx_polynya
+
+ out_nc.close()