source: BOL/LMDZ_Setup/era2gcm.sh @ 5859

#!/bin/bash

. lmdz_env.sh

set -eu

#-----------------------------------------------------------------------------
#
#  Script for interpolating monthly ERA* files to the LMDZ grid;
#  it creates the folder structure required by the LMDZ_Setup scripts
#
#  NB: for "cleanest" guiding, the 1st step of the next month should be added at the end of each month.
#    If you need such precision, you should use the scripts here  - AND adjust "ERADIR" in script_SIMU !
#    wget http://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/INTERP_NUDGE
#      documented on LMDZpedia (search for "Guidage" )
#
#-----------------------------------------------------------------------------
# Requires : netcdf, ferret, nco, cdo
#-----------------------------------------------------------------------------

#------------------------------------------------
# User choices
#------------------------------------------------
# Periode :
mth_i=200001
mth_f=200212
#
vars="u v ta q"
vars="u v"
#
# Choix des fichiers de guidage ("rea"nalyses) : ERA5, ERAI, OPERA
rea="ERAI"
resol_rea=075

# Location of shared files
ERA_PATH="/gpfsstore/rech/psl/rpsl376/ergon/$rea/"

# New location.
# The idea is to get files from the LMDZ web site and store them locally
# once for all.
# http://lmdz.lmd.jussieu.fr/pub/3DInputData/$rea
ERA_PATH="$LMDZ_INIT/3DInputData/ReAnalyses/$rea"

GRID_DIR=./INIT

#-----------------------------------------------------------------------------
#Utilite du block suivant a re-examiner :
#-----------------------------------------------------------------------------
#L'utilisateur a maintenant des choix a faire en plus de mth* : type de guidage, et "rea"
# Alors c'est plus facile d'editer&modifier le script, puis lancer avec ./era2gcm.sh, que de donner tous les params
# TODO check w/ Adriana if we can remove those and put them as inline args (above)
while (($# > 0))
   do
   case $1 in
     "-h") cat <<........fin
           ./era2gcm.sh [-mthini initial_month] [-mthend final_month] [-grid_dir directory_containing_grille_gcm.nc]
........fin
     exit ;;
     "-grid_dir") GRID_DIR=$2 ; shift ; shift ;;
     "-mthini") mth_i=$2 ; shift ; shift ;;
     "-mthend") mth_f=$2 ; shift ; shift ;;
     *) $0 -h ; exit 1
   esac
done
#--Fin du bloc a examiner --------------------------------------------

tmin=1
resol=grilles_gcm.nc

GRID_FI=${GRID_DIR}/${resol}
if [ ! -f "$GRID_FI" ] ; then
   echo "Le fichier $GRID_DIR/$resol est nécessaire; créer le fichier $GRID_FI avec grilles_gcm_netcdf.e"
   exit 1
fi
mth=$mth_i


#####################################################################
#  LOOP ON MONTHS
#####################################################################
while (( mth <= mth_f )) ; do

   echo "mth $mth"
   mois=$(echo "$mth" | cut -c 5-)
   an=$(echo "$mth" | cut -c -4)
   ndays=( 31 28 31 30 31 30 31 31 30 31 30 31 )
   months=( jan feb mar apr may jun jul aug sep oct nov dec )
   if [ $(( an % 4 )) = 0 ] ; then ndays[1]=29 ; fi
   imois=$(( ${mois#0} - 1 ))
   month=${months[$imois]}
   nday=${ndays[$imois]}
   tmax=$(( nday * 4 ))
   echo "PARAMETRES CALENDAIRES $imois $month $nday $tmax"


   iip1=$(ncdump -h "$GRID_FI" | grep lonu | head -1 | awk ' { print $3 } ')
   jjm=$(ncdump -h "$GRID_FI" | grep latv | head -1 | awk ' { print $3 } ')
   (( jjp1 = jjm + 1 ))
#   \rm t2.nc ua.nc va.nc sst.nc u.nc v.nc T.nc ts.nc

#####################################################################
# Choix de la periode temporelle
#####################################################################

   t0="l=$tmin"
   t1tn="l=${tmin}:${tmax}"

#####################################################################
# Lien avec les fichiers netcdf contenant les d0 ECMWF 
#####################################################################
   echo "-------- liens de telechargement a actualiser ----"
   if [ "$rea" = "ERA5" ] ; then
     if [[ $an -ge 2022 ]] ; then
      ANA_DIR="$ERA_PATH/NETCDF/GLOBAL_025/hourly"
      suf="ap1e5.GLOBAL_025"
     else
      ANA_DIR="$ERA_PATH/NETCDF/GLOBAL_025/4xdaily"
      suf="aphe5.GLOBAL_025"
     fi
   elif [ "$rea" = "ERAI" ] ; then
      ANA_DIR="$ERA_PATH/NETCDF/GLOBAL_$resol_rea/4xdaily"
      suf="aphei.GLOBAL_$resol_rea"
   fi

   if [ ! -d $ANA_DIR ] ; then echo Directory $ANA_DIR does not exist ; exit 1 ; fi

   #################################################################################
   # Loop on variables among u, v, ta, q
   #################################################################################

   for var in $vars ; do

     case $var in
         ta) filename=T.nc ;;
         q) filename=hur.nc ;;
         *) filename=$var.nc
     esac
     outd=GUIDE_${rea}/$an/$mois
     out_file=$outd/$filename

     if [ ! -f $out_file ] ; then
        mkdir -p $outd
        file_rea="$ANA_DIR/AN_PL/$an/$var.$an$mois.$suf.nc"

        ############################################################################
        # Automatically downloading reanalysis if needed
        ############################################################################
        if [ "$( echo $ERA_PATH | grep 3DInputData )" != "" ] ; then
             if [ ! -f $file_rea ] ; then
                mkdir -p $( dirname $file_rea )
                cd $( dirname $file_rea )
                wget http://lmdz.lmd.jussieu.fr/pub/3DInputData/$( echo $file_rea | sed -e 's/^.*3DInputData//' )
                cd -
             fi
         fi
         ############################################################################
         # Checking the avaibility of reanalysis native file
         ############################################################################
         if [ ! -f $file_rea ] ; then echo Reanalysis file $file_rea missing ; exit 1 ; fi

         ############################################################################
         # Special treatments
         ############################################################################
         if [ "$rea" = "ERAI" ] ; then
             # original data in "short" should be changed to "float" for ferret
             ncap2 -s $var'=float('$var')' $file_rea -O tmp_in.nc
         elif [ "$rea" = "ERA5" -a $an -ge 2022 ] ; then
             # Recent analysis available on a daily basis
             cdo selhour,0,6,12,18 $file_rea -O tmp_in.nc
         else
             ln -sf $file_rea tmp_in.nc
         fi

         ############################################################################
         # Dimensions of horizontal grid
         ############################################################################
         case $var in
            v) imjm="i=1:$iip1,j=1:$jjm" ;;
            *) imjm="i=1:$iip1,j=1:$jjp1"
         esac

         ############################################################################
         # Grid to be used from grilles_gcm.nc
         ############################################################################
         case $var in
            u|v) grid=grille_$var ;;
            *) grid=grille_s
         esac

         ############################################################################
         # Interpolating with ferret
         ############################################################################
         cat <<.........eod......... >| tmp.jnl
         ! NB : Augmenter la memoire (plus de 512) peut faire planter
         set memory/size=512
         use "$GRID_FI"
         use tmp_in.nc
         define axis/t="1-${month}-${an}:00:00":"${nday}-${month}-${an}:18:00":6/units=hours thour
         ! Pour regrid horizontal on utilise la variable grille_s(lonv,latu) du fichier grilles_gcm.nc
         !! Alors il faut renommer la variable à la fin (ncrename), car le code cherche "UWND", "VWND".
         save/clobber/file="tmp_out.nc" $var[d=2,gxy=${grid}[d=1],$imjm,$t0,gt=thour@asn]
         repeat/$t1tn save/file="tmp_out.nc"/append $var[d=2,gxy=${grid}[d=1],$imjm,gt=thour@asn]
.........eod.........

         set +e
         ferret -nojnl <<.........eod.........
         go tmp.jnl
         quit
         fi
.........eod.........
         if [ $? != 0 ] ; then
            echo Something went wrong when running ferret with script tmp.jnl :
            cat tmp.jnl
            echo On $PWD
            exit 1
         else
            set -e
            mv tmp_out.nc $out_file
         fi

        ############################################################################
        # Changing variable names
        ############################################################################
         case $var in
            u) out_var=UWND ;;
            v) out_var=VWND ;;
            r) out_var=RH ;;
            ta) out_var=AIR
         esac
         #Ferret a ecrit $varu en majuscules, l'equivalent en bash est ${varu^^} 
         #   (Note : inversement, ${varu,,} passe $varu de majuscules en minuscules)
         ncrename -v ${var^^},$out_var $out_file
     fi
   done


   ############################################################################
   # Updating month counter
   ############################################################################
   echo AN MTH $an $mois
   (( mth = $mth + 1 ))
   if [ $mois = 12 ] ; then
      (( an = $an + 1 ))
      mth=${an}01
   fi

done
ln -sf GUIDE_$rea GUIDE