EnKF/interpolate_observations/esa/ESA_to_ORCA1.bash

#!/bin/bash
#SBATCH -n 1
#SBATCH -t 96:00:00
#SBATCH -J interpolate_ESA
#SBATCH -o /esnas/obs/esa/scripts/prepare_obs_EnKF/interpolation/slurm-%j.out
#SBATCH -e /esnas/obs/esa/scripts/prepare_obs_EnKF/interpolation/slurm-%j.err

# Original: Pierre Mathiot, 2011
# Update  : Francois Massonnet, 2013
# Update at BSC: Francois Massonnet, 2016
#
# Interpolation of ESA sea ice concentration to ORCA grid
# Particular attention is paid to "flag" values. Concentration
# is only considered in "nominal" regions, that is, not the coastal
# areas, not the gap filling area around the north pole, not questionable 
# data etc.
#
#
# Questions: francois.massonnet@uclouvain.be
#            francois.massonnet@bsc.es

module load netCDF-Fortran/4.2-foss-2015a

#set -x
set -o nounset
set -o errexit

yearb=1992    # Years to process
yeare=2008           

grid=ORCA1    # Grid type

nbsmooth=2    # This is important. We need to interpolate the observations
              # on the model grid. However the model is on a coarser grid.
              # Therefore, before interpolation, we will smooth the input
              # observational data. This variable tells how many times we
              # have to apply the "smooth9" function of cdo. This function
              # makes a 2-D smoothing of each grid point by weighting with
              # itself and its neighbours. 
              # Roughly speaking, the variable nbsmooth should be the ratio
              # of model resolution to observational resolution.
              # At ORCA1, resolution is ~50 km at the poles and input data
              # from ESA is 25 km, hence a ratio of 2
              # Set nbsmooth to 0 if you don't want smoothing.


sosiedir=$HOME/sosie-2.6.2/bin

mask=/esnas/autosubmit/con_files/mesh_mask_nemo.Ec3.2_O1L75.nc
sourcedir=/esnas/obs/esa/original_files/seaice # Where original files are located
outdir=/esnas/obs/esa/interpolated/

scratchdir=/esnas/scratch/$USER/TMP

# Create a directory to work
# --------------------------
tmpdir=$scratchdir/TMP_${RANDOM}

echo "TMPDIR IS >>>>>> $tmpdir <<<<<<<"
mkdir -p $tmpdir
cp namelist* $tmpdir
cp sosie_mapping_ESA*.nc $tmpdir

cd $tmpdir

# 1. Get the model grid and the mask
if [ ! -f ${mask} ]
then
	echo "${mask} does not exist."
	exit
fi

ln -sf ${mask} mask_out.nc 

# Link sosie
ln -sf $sosiedir/sosie.x .

for year in `seq ${yearb} ${yeare}`
do
  for month in 01 02 03 04 05 06 07 08 09 10 11 12
  do
    echo $year-${month}
    
    rm -f listfiles.txt
    ls $sourcedir/${year}/$month/ESACCI-SEAICE-L4-SICONC-SSMI-NH25kmEASE2-????????-fv01.11.nc > list_files.txt
  
    for file in `cat list_files.txt`
    do
      tag=`basename $file | sed -e "s/ESACCI-SEAICE-L4-SICONC-SSMI-NH25kmEASE2-//" -e "s/-fv01.11.nc//"`
      echo $tag

      # Check that the file is also available for the other hemisphere
      nfile=`ls $sourcedir/${year}/${month}/ESACCI-SEAICE-L4-SICONC-SSMI-?H25kmEASE2-${tag}-fv01.11.nc | wc -l`
  
      if [[ $nfile == 2 ]]
      then
        echo "Two files were found for the time stamp ${tag}!"
        
        for hemi in nh sh
        do
          HEMI=`echo $hemi | tr '[:lower:]' '[:upper:]'`

          filein=$sourcedir/${year}/${month}/ESACCI-SEAICE-L4-SICONC-SSMI-${HEMI}25kmEASE2-${tag}-fv01.11.nc
          cp -f $filein data_in_${hemi}-${tag}.nc
          
          # We only take points that have nominal (status_flag = 0) quality
          cdo chname,status_flag,mask -eqc,0  \
              -selvar,status_flag data_in_${hemi}-${tag}.nc mask_in_${hemi}-${tag}.nc
  
          # Smoothing input data
          if [[ $nbsmooth == 0 ]]
          then
            echo "WARNING - YOU DID NOT SMOOTH INPUT DATA"
            echo "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
            mv data_in_${hemi}-${tag}.nc data_in_smooth_${hemi}-${tag}.nc 
          else
            mv data_in_${hemi}-${tag}.nc data_in_smooth0_${hemi}-${tag}.nc
            for jsmooth in `seq 1 $nbsmooth`
            do
              cdo smooth9 data_in_smooth$(( $jsmooth - 1 ))_${hemi}-${tag}.nc \
                          data_in_smooth${jsmooth}_${hemi}-${tag}.nc
            done
            mv data_in_smooth${jsmooth}_${hemi}-${tag}.nc data_in_smooth_${hemi}-${tag}.nc
          fi # if smoothing is necessary
  
  
          # Interpolation of land-sea mask
          echo "Interpolation of land-sea mask"
          sed -e "s/TTAARRGGEETT/${grid}/" \
              -e "s/HHEEMMII/${hemi}/" \
              -e "s/TTAAGG/${tag}/" \
                  namelist_ESA-mask > namelist
  
          ./sosie.x
     
          # Since we just did an interpolation, the mask will not be sharp:
          # There will be values not equal to 1. Let's remask the result
          # The threshold is pretty arbitrary here, but to be conservative
          # we just ignore points different from one.
  
          cdo eqc,1 tmask_ESA-${hemi}-${grid}_${tag}.nc ESA-${hemi}_mask_on_${grid}_${tag}.nc
  
  
          # 2. Interpolation of concentration
          sed -e "s/TTAARRGGEETT/${grid}/" \
              -e "s/HHEEMMII/${hemi}/" \
              -e "s/TTAAGG/${tag}/" \
                  namelist_ESA-conc > namelist
          
          ./sosie.x
  
          #    Mask the interpolated ice concentration using the interpolated mask.
          #    Also, put missing values returned by sosie (<0) to zero, so that 
          #    we'll be able to add northern and southern hemispheres later
  
          cdo setmisstoc,0 \
             -mul -selvar,at_i  at_i_ESA-${hemi}-${grid}_${tag}.nc \
                  -selvar,tmask ESA-${hemi}_mask_on_${grid}_${tag}.nc \
                                at_i_ESA-${hemi}-${grid}_${tag}_masked.nc
  
          # 3. Interpolation of standard deviation
          #
          # This is a point to worry about, although there is not enough
          # at this stage to do better. Interpolation of second-order moments
          # should account for the fact that errors are correlated, however
          # we don't have this information from the OSI-SAF product. 
          # The best we can do is to assume that these errors are not correlated
          # and can be interpolated as is.
  
          sed -e "s/TTAARRGGEETT/${grid}/" \
              -e "s/HHEEMMII/${hemi}/" \
              -e "s/TTAAGG/${tag}/" \
                  namelist_ESA-error > namelist
  
          ./sosie.x
  
          #    Mask the interpolated error using the interpolated mask
          cdo setmisstoc,0 \
              -mul -selvar,error_at_i error_at_i_ESA-${hemi}-${grid}_${tag}.nc \
                  -selvar,tmask ESA-${hemi}_mask_on_${grid}_${tag}.nc \
                             error_at_i_ESA-${hemi}-${grid}_${tag}_masked.nc
        done # for each hemisphere
        
        # Merge North and South data in one file
        for var in at_i error_at_i
        do
          ncbo -O -F -v ${var} --op_typ=add ${var}_ESA-nh-${grid}_${tag}_masked.nc \
                                            ${var}_ESA-sh-${grid}_${tag}_masked.nc \
                                            ${var}_ESA-${grid}_${tag}.nc
        done
  
        # Dump errors into concentration files
        ncks -F -A -v error_at_i error_at_i_ESA-${grid}_${tag}.nc at_i_ESA-${grid}_${tag}.nc
        mv at_i_ESA-${grid}_${tag}.nc $outdir
        # Do a bit of cleaning here!
        rm -f at_i_ESA-?h-${grid}_${tag}_masked.nc                    \
              error_at_i_ESA-?h-${grid}_${tag}_masked.nc              \
              error_at_i_ESA-?h-${grid}_${tag}_masked.nc              \
              error_at_i_ESA-${grid}_${tag}.nc                        \
              ESA-?h_mask_on_${grid}_${tag}.nc                        \
              error_at_i_ESA-?h-${grid}_${tag}.nc                     \
              at_i_ESA-?h-${grid}_${tag}.nc                           \
              tmask_ESA-?h-${grid}_${tag}.nc                          \
              data_in_smooth*${tag}.nc                                   \
              mask_in_?h-${tag}.nc                                       
      else
        echo "Only one file was found for the time stamp $tag"
        echo "This is not coded yet"
      fi # if there are two files
    done # for each file of the year
  done # for each month
done # for each year