diff --git a/smoothing/README b/smoothing/README
new file mode 100644
index 0000000000000000000000000000000000000000..2b629f9e2e2077d1b5fb7f57fa69f0f3fbe3a857
--- /dev/null
+++ b/smoothing/README
@@ -0,0 +1,17 @@
+Smoothing package:
+------------------
+This is a suite of programs to obtain smoothed restarts from a given GLORYS restart.
+The structure of the code is as follows:
+run_smoothing.bash
+    |
+    ---> smooth_fields.R
+    |
+    ---> replace_variables.bash
+
+The bash-script "run_smoothing.bash" sets the parameters to call the respective R and bash scripts and runs them.
+
+Instructions for running:
+-------------------------
+For an ORCA1L46 setup, and two start-dates, the script has run on moore
+ in an interactive queue specifying:  qrsh -l h_vmem=20G,s_rt=02:00:00,h_rt=02:00:00
+
diff --git a/smoothing/replace_variables.bash b/smoothing/replace_variables.bash
new file mode 100755
index 0000000000000000000000000000000000000000..c5bc1fbaee1d6cbb39be2c7565f31678dfd64751
--- /dev/null
+++ b/smoothing/replace_variables.bash
@@ -0,0 +1,23 @@
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
+#
+# This is a bash script that needs to run after the smoothing program "smooth_fields.R"
+# It replaces the fields in a file by the corresponding smoothed fields, obtained from "smooth_fields.R"
+#
+# Author: I. Andreu-Burillo (IC3/CFU, December 2012)
+#
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
+
+
+# Work one start-date in the list at a time
+for sdate in ${list_of_dates}; do
+  cd ${outpath}/$sdate
+  cp ${inpath}${sdate}_${setup}.nc ${outpath}/$sdate/restart_glorys2v1_oce_${sdate}_${setup}.nc
+  ncwa -O -a t restart_glorys2v1_oce_${sdate}_${setup}.nc restart_glorys2v1_oce_${sdate}_${setup}.nc || echo "This file does not have a degenerate t dimension"
+  for cvar in ${varlist}; do
+    ncks -A -v ${cvar} restart_glorys2v1_oce_${sdate}_${setup}_${cvar}.nc restart_glorys2v1_oce_${sdate}_${setup}.nc
+  done
+  mv ${outpath}/$sdate/restart_glorys2v1_oce_${sdate}_${setup}.nc ${outpath}/restart_glorys2v1_oce_${sdate}_${setup}.nc && rm -R ${outpath}/$sdate
+done
+
+echo "end of creation of smoothed restart-files"
+echo "Please, fetch the files under " ${outpath}
diff --git a/smoothing/run_smoothing.bash b/smoothing/run_smoothing.bash
new file mode 100755
index 0000000000000000000000000000000000000000..0c02e318ef00153a19a11215f7f895e048fc4031
--- /dev/null
+++ b/smoothing/run_smoothing.bash
@@ -0,0 +1,107 @@
+#!/bin/bash
+set -evx
+
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
+#
+# This script runs the suite that produces smoother restarts
+# from a given GLORYS restart. Up to present, it has been used when the
+# model crashes due to too low salinities close to coasts.
+#
+# To run, adapt the parameters in sections "To be set by the user"
+# For an ORCA1L46 setup, and two start-dates, the script has run on moore
+# in an interactive queue specifying:
+# qrsh -l h_vmem=20G,s_rt=02:00:00,h_rt=02:00:00
+#
+# Author: isabel.andreu-burillo@ic3.cat (IC3/CFU)
+# Created: December 2012
+# Modified: July 2012
+#
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
+
+## To be set by the user >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+## Switches to turn on/off (1/0) each step
+# Smoothing step, it smoothes the variables in varlist and creates new NetCDF files
+first_step=1
+# Replacing step, it replaces the raw variables in the original file,
+# with the smoothed variables in a new restart file
+second_step=1
+
+## NEMO grid indices for which the model crashes and around which there needs to be smoothing
+ip=193
+jp=635
+## Number of points around the problem point (ip,jp) defining the region of smoothing
+delt=5
+### List of variables to be smoothed
+varlist=(tn tb sn sb un ub vn vb)
+## Year(s), month(s), day(s) of the restart
+yearlist=(1992 1993)
+monthlist=(11 11)
+daylist=(03 02) 
+list_of_dates=()
+setup=ORCA025L46
+inpath=/cfu/scratch/iandreu/Smooth_restarts/restart_glorys2v1_oce_
+outpath=/scratch/$USER/SmoothRestarts/
+## <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< To be set by the user
+
+mkdir -p /scratch/$USER/SmoothRestarts/
+# Create the list of dates to be treated
+nsdatesm1=$((${#yearlist[@]}-1))
+echo $nsdatesm1
+for (( jloop=0;jloop<=${nsdatesm1};jloop=$(($jloop+1)) )) ; do
+  echo ${yearlist[$jloop]}
+  yyyy=${yearlist[$jloop]}
+  month=${monthlist[$jloop]}
+  mm=`printf "%2.2i" ${month}`
+  day=${daylist[$jloop]}
+  dd=`printf "%2.2i" ${day}`
+  list_of_dates=${list_of_dates}" "$yyyy$mm$dd
+  echo ${list_of_dates}
+  mkdir -p /scratch/$USER/SmoothRestarts/$yyyy$mm$dd
+
+  # Producing the smooth fields from the original ones
+  cvlist=`echo ${varlist[@]} | sed -e s/\ /\',\'/g`
+  if [ ${first_step} -eq 1 ]; then
+    cat > header_1.R << EOF
+    ## Clean up before starting
+    rm(list=ls())
+    ip=$ip
+    jp=$jp
+    delt=$delt
+    yyyy=$yyyy
+    month=$month
+    day=$day
+    cvarlist=c('${cvlist[@]}')
+    list_of_dates="${listofdates}"
+    setup="${setup}"
+    inpath="${inpath}"
+    outpath="${outpath}"
+EOF
+    cat header_1.R smooth_fields.R > first_step.R
+    R CMD BATCH ./first_step.R && {
+    rm header_1.R*
+    rm first_step.R*
+    }
+  fi
+done
+
+# Replacing the original fields with the smoothed fields in a new restart file
+if [ ${second_step} -eq 1 ]; then
+  cat > header_2.bash << EOF
+  #!/bin/bash
+  set -evx
+  ip=$ip
+  jp=$jp
+  delt=$delt
+  setup=${setup}
+  inpath=$inpath
+  outpath=$outpath
+  list_of_dates="`echo ${list_of_dates[@]}`"
+  varlist="${varlist[@]}"
+EOF
+  cat header_2.bash replace_variables.bash > second_step.bash
+  chmod u+x second_step.bash
+  ./second_step.bash && {
+  rm header_2.bash ;
+  rm second_step.bash;
+  }
+fi
diff --git a/smoothing/smooth_fields.R b/smoothing/smooth_fields.R
new file mode 100644
index 0000000000000000000000000000000000000000..acfcca5a475376f1289945a3b0989189d910323c
--- /dev/null
+++ b/smoothing/smooth_fields.R
@@ -0,0 +1,116 @@
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
+#
+# R program to smooth given ocean fields around a point (ip,jp)
+# throughout the water column (only wet points)
+#
+# Once the program has run, it must be followed by the bash script "replace_variables.bash"
+# under directory /home/iandreu/Scripts/Shell
+#
+# Author: I. Andreu-Burillo (IC3/CFU, December 2012)
+#
+# NOTE:
+# =====
+# At this moment, the smoothing is implemented in Loop A,
+# which corresponds to a coastal region with the continent on the West.
+# The smoothing is performed from the open ocean towards the coast,
+# to bring smoother information into the smaller scales.
+# If smoothing needs to be applied in a coastal region with a different orientation,
+# a different loop will need to be written and a switch implemented.
+# 
+# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
+
+## Load necessary libraries
+library(ncdf)
+mm= sprintf('%2.2i',month)
+dd= sprintf('%2.2i',day)
+
+## Treat the different variables that need smoothing one after the other
+for (cvar in cvarlist) {
+
+  ##
+  ## Read in the variable from the original restart file
+  ##
+  indata.nc=open.ncdf(paste(inpath,yyyy,mm,dd,"_",setup,".nc",sep=''))
+  x=get.var.ncdf(indata.nc,"x")
+  y=get.var.ncdf(indata.nc,"y")
+  z=get.var.ncdf(indata.nc,"z")
+  var_field=get.var.ncdf(indata.nc,cvar)
+  close.ncdf(indata.nc)
+  jpk=length(z)
+  print(jpk)
+
+  ## 3D variables come as var_field[x,y,z]
+
+  ## Recalculate the variable value at each point as the weighted average of 
+  ## it's value with the surrounding points.
+
+  ## NOTE:
+  ## =====
+  ## Loop A is a particular application to the Vietnam Coast.
+  ## We want to start the smoothing off-shore and sweep towards the shore, 
+  ## to bring the information of higher salinities towards the coast.
+  ## This is why both i and j indices decrease in the loop.
+  ## Loop A may have to be adapted to the specific configuration in case
+  ## of neighbouring coasts.
+
+  for (jk in seq(1,jpk,1)) {
+    ## Loop A --> From the NW towards SE
+    for (jj in seq((jp+delt),(jp-delt),-1)) {
+      for (ji in seq((ip+delt),(ip-delt),-1)) {
+        ## Treating only wet points
+        if (var_field[ji,jj,jk] > 0) {
+          icount = 1
+          zvar = 0
+          ## We average around each point that requires smoothing
+          for (jjj in seq((jj +1),(jj - 1), -1)) {
+            for (jii in seq((ji + 1),(ji - 1), -1)) {
+              if (var_field[jii,jjj,jk] > 0) {
+                zvar=zvar + var_field[jii,jjj,jk]
+                icount = icount + 1
+              }
+            }
+          }
+          max_var=max(var_field[seq((ji - 1),(ji + 1),1),seq((jj -1),(jj + 1),1),jk])
+          var_field[ji,jj,jk]=(var_field[ji,jj,jk]+zvar)/icount
+          if ( cvar == 'sn' | cvar == 'sb' ) {
+            # If the field is salinity and we are under the threshold "1",
+            # overwrite this with a higher value
+            if  (var_field[ji,jj,jk] < 1 ) {
+              var_field[ji,jj,jk] = max_var
+              icount = 0
+              zvar=0
+              for (jjj in seq((jj +1),(jj - 1), -1)) {
+                for (jii in seq((ji + 1),(ji - 1), -1)) {
+                  if (var_field[jii,jjj,jk] > 0) {
+                    zvar=zvar + var_field[jii,jjj,jk]
+                    icount = icount + 1
+                  }
+                }
+              }
+              var_field[ji,jj,jk]=(var_field[ji,jj,jk]+zvar)/icount
+            }
+          }
+          icount = 0
+          zvar=0
+        }
+      }
+    }
+    ## <-- end of Loop A
+  }
+
+  ##
+  ## Define variables for the new NetCDF file containing the smooth field
+  ##
+  x_new=dim.def.ncdf("x",'',x)
+  y_new=dim.def.ncdf("y",'',y)
+  z_new=dim.def.ncdf("z",'',z)
+
+  var_new=var.def.ncdf(cvar,'',list(x_new,y_new,z_new),1.e+20)
+
+  ## A different NetCDF file for each (cvar) variable in the (cvarlist) list
+  nc=create.ncdf(paste(outpath,yyyy,mm,dd,"/restart_glorys2v1_oce_",yyyy,mm,dd,"_",setup,"_",cvar,".nc",sep=''),var_new)
+  put.var.ncdf(nc,var_new,var_field)
+  close.ncdf(nc)
+  ## Clear up some space and key variables
+  rm(x, y, z,var_field, x_new, y_new, z_new, var_new)
+}