source: BOL/LMDZ_Setup/era2gcm_tuto.sh @ 5441

#!/bin/bash

. lmdz_env.sh

#-----------------------------------------------------------------------------
#
#  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
#-----------------------------------------------------------------------------
#@JZ module load netcdf/4.7.2-mpi-cuda ferret nco cdo
#
# Check&modify section "User choices", then run the script with ./era2gcm.sh
#-----------------------------------------------------------------------------
# ATTENTION : VERIFIEZ 
#  --> l'accès (repositoire ou telechargement) a la réanalyse ERA souhaitée (variable ANA_DIR)
#    Contact pour demander acces aux fichiers sur Jean-Zay : 
#        Sophie Bouffiès-Cloché, IPSL : Sophie.Bouffies-Cloche@ipsl.fr

#  --> la disponibilite des fichiers pour la periode souhaitee (variables file*) 
#-----------------------------------------------------------------------------

#------------------------------------------------
# User choices
#------------------------------------------------
# Periode :
mth_i=202201
mth_f=202201
#
# Guidage en vent(u&v) et/ou temperature, humidite
guide_uv="y"
guide_t="n"
guide_q="n"
#
# Choix des fichiers de guidage ("rea"nalyses) : ERA5, ERAI, OPERA
rea="ERA5"
#------------------------------------------------

#@JZ ERA5_PATH="/gpfsstore/rech/psl/rpsl376/ergon/ERA5/"
#@JZ ERAI_PATH="/gpfsstore/rech/psl/rpsl376/ergon/ERAI/"
#@ADS ERA5_PATH="/lus/work/CT1/cad15221/abarral/ERA_TMP/"  # /!\ temp for test only !
#@ADS ERAI_PATH=""  # /!\ temp for tests only!
#@ADS echo "/!\ Nudging files are not available on Adastra for now !"; exit 1  # comment this line for tests
#@ADS SCRATCH=$SCRATCHDIR

set -u  # raise error if path if unset on machine
echo "Paths: $ERA5_PATH $ERAI_PATH $SCRATCH"
set +u

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 ;;
     "-grille_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


#####################################################################
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="$ERA5_PATH/NETCDF/GLOBAL_025/hourly"
      suf="ap1e5.GLOBAL_025"
     else
      ANA_DIR="$ERA5_PATH/NETCDF/GLOBAL_025/4xdaily"
      suf="aphe5.GLOBAL_025"
     fi
   elif [ "$rea" = "ERAI" ] ; then
      ANA_DIR="$ERAI_PATH/NETCDF/GLOBAL_1125/4xdaily"
      suf="aphei.GLOBAL_1125"
   fi

varu=u
varv=v
vart=ta # peut etre parfois juste "t"
varq=q
if [ "$rea" = "ERA5" ] ; then varq=r ; fi
#varp=msl

if [ "$rea" = "ERAI" ] ; then
  # variables en format "short" doivent etre transformees en "float" via NCO 
  # This is done here with ncap2 ; also possible: "ncpdq --overwrite --unpack fin.nc fout.nc"
  fushort="$ANA_DIR/AN_PL/$an/u.$an$mois.$suf.nc"
  fvshort="$ANA_DIR/AN_PL/$an/v.$an$mois.$suf.nc"
  ftshort="$ANA_DIR/AN_PL/$an/ta.$an$mois.$suf.nc"
  fqshort="$ANA_DIR/AN_PL/$an/r.$an$mois.$suf.nc"
  fileu="$SCRATCH/u.$an$mois.$suf.nc"
  filev="$SCRATCH/v.$an$mois.$suf.nc"
  filet="$SCRATCH/ta.$an$mois.$suf.nc"
  fileq="$SCRATCH/r.$an$mois.$suf.nc"
  if [ "$guide_uv" = "y" ] ; then
    ncap2 -s 'u=float(u)' $fushort $fileu
    ncap2 -s 'v=float(v)' $fvshort $filev
  fi
  if [ "$guide_t" = "y" ] ; then ncap2 -s 'ta=float(ta)' $ftshort $filet ; fi
  if [ "$guide_q" = "y" ] ; then ncap2 -s 'q=float(q)' $fqshort $fileq ; fi
elif [ "$rea" = "ERA5" -a $an -ge 2022 ] ; then
  echo Extract 0,6,12,18 hours from ERA5 hourly files
  fu1h="$ANA_DIR/AN_PL/$an/u.$an$mois.$suf.nc"
  fv1h="$ANA_DIR/AN_PL/$an/v.$an$mois.$suf.nc"
  ft1h="$ANA_DIR/AN_PL/$an/ta.$an$mois.$suf.nc"
  fq1h="$ANA_DIR/AN_PL/$an/r.$an$mois.$suf.nc"
  fileu="$SCRATCH/u.$an$mois.$suf.nc"
  filev="$SCRATCH/v.$an$mois.$suf.nc"
  filet="$SCRATCH/ta.$an$mois.$suf.nc"
  fileq="$SCRATCH/r.$an$mois.$suf.nc"
  if [ "$guide_uv" = "y" ] ; then  
    cdo selhour,0,6,12,18 $fu1h $fileu
    cdo selhour,0,6,12,18 $fv1h $filev
  fi
  if [ "$guide_t" = "y" ] ; then cdo selhour,0,6,12,18 $ft1h $filet ; fi
  if [ "$guide_q" = "y" ] ; then cdo selhour,0,6,12,18 $fq1h $fileq ; fi
else 
 fileu="$ANA_DIR/AN_PL/$an/u.$an$mois.$suf.nc"
 filev="$ANA_DIR/AN_PL/$an/v.$an$mois.$suf.nc"
 filet="$ANA_DIR/AN_PL/$an/ta.$an$mois.$suf.nc"
 fileq="$ANA_DIR/AN_PL/$an/r.$an$mois.$suf.nc"
fi

# verifier disponibilite des fichiers
if [ "$guide_uv" = "y" ] ; then  ls $fileu ; ls $filev ; fi
if [ "$guide_t" = "y" ] ; then  ls $filet ; fi
if [ "$guide_q" = "y" ] ; then  ls $fileq ; fi

outd=GUIDE_${rea}/$an/$mois
mkdir -p $outd

# effacer lien vers repertoire des vents ERA interpoles (utilise dans script_SIMU) s'il existe deja
# recreer lien vers le repertoire "rea" choisi
# NB : si GUIDE existe en tant que repertoire, non pas comme lien, il n'est pas affecte
rm -f GUIDE
ln -s GUIDE_${rea} GUIDE

ij="i=1:$iip1,j=1:$jjp1"
ijm="i=1:$iip1,j=1:$jjm"

###################################################################3
# scripts ferret pour interpolation
###################################################################3

# --- guide uv ---

if [ "$guide_uv" = "y" ] ; then


cat << eod >| tmp_uv.jnl
! NB : Augmenter la memoire (plus de 512) peut faire planter
set memory/size=512

use "$GRID_FI"
use "$fileu"
use "$filev"

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


!! ATTENTION : pour ecrire le fichier en simple, non pas double precision :
!! on utilise directement $varu = "float" dans le fichier d'origine (non pas une variable definie avec "let" dans ferret)
!! Alors il faut renommer la variable à la fin (ncrename), car le code cherche "UWND", "VWND".
save/clobber/file="$outd/u.nc" $varu[d=2,gxy=grille_u[d=1],$ij,$t0,gt=thour@asn]
repeat/$t1tn save/file="$outd/u.nc"/append $varu[d=2,gxy=grille_u[d=1],$ij,gt=thour@asn]
save/clobber/file="$outd/v.nc" $varv[d=3,gxy=grille_v[d=1],$ijm,$t0,gt=thour@asn]
repeat/$t1tn save/file="$outd/v.nc"/append $varv[d=3,gxy=grille_v[d=1],$ijm,gt=thour@asn]

eod

ferret -nojnl <<eod
go tmp_uv.jnl
quit
eod
#Ferret a ecrit $varu en majuscules, l'equivalent en bash est ${varu^^} 
#   (Note : inversement, ${varu,,} passe $varu de majuscules en minuscules)
ncrename -v ${varu^^},UWND $outd/u.nc 
ncrename -v ${varv^^},VWND $outd/v.nc

fi  # ------- fin guide_uv --------

# --- guide_t  ---

if [ "$guide_t" = "y" ] ; then

cat << eod >| tmp_t.jnl
set memory/size=512

use "$GRID_FI"
use "$filet"

define axis/t="1-${month}-${an}:00:00":"${nday}-${month}-${an}:18:00":6/units=hours thour

save/clobber/file="$outd/T.nc" $vart[d=2,gxy=grille_s[d=1],$ij,$t0,gt=thour@asn]
repeat/$t1tn save/file="$outd/T.nc"/append $vart[d=2,gxy=grille_s[d=1],$ij,gt=thour@asn]
eod

ferret -nojnl <<eod
go tmp_t.jnl
quit
eod
ncrename -v ${vart^^},AIR $outd/T.nc

fi  # ------- fin guide_t --------

# --- guide_q  ---

if [ "$guide_q" = "y" ] ; then

cat << eod >| tmp_q.jnl
set memory/size=512

use "$GRID_FI"
use "$fileq"

define axis/t="1-${month}-${an}:00:00":"${nday}-${month}-${an}:18:00":6/units=hours thour

save/clobber/file="$outd/hur.nc" $varq[d=2,gxy=grille_s[d=1],$ij,$t0,gt=thour@asn]
repeat/$t1tn save/file="$outd/hur.nc"/append $varq[d=2,gxy=grille_s[d=1],$ij,gt=thour@asn]
eod

ferret -nojnl << eod
go tmp_q.jnl
quit
eod
ncrename -v ${varq^^},RH $outd/hur.nc 

fi  # ------- fin guide_q --------


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

done