source: BOL/LMDZ_Setup/era2gcm_tuto.sh @ 4948

#!/bin/bash

#-----------------------------------------------------------------------------
#
#  Script for interpolating monthly ERA* files to the LMDZ grid;
#  it creates the folder structure required by the tutorial_prod 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
#-----------------------------------------------------------------------------
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'acces (repositoire ou telechargement) a la reanalyse ERA souhaitee (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
#------------------------------------------------

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
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
   esac
done
#--Fin du block 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 necessaire
   echo creer le fichier $GRID_FI avec grilles_gcm_netcdf.e
   exit
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 ))
   (( iim  = $iip1 - 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='/gpfsstore/rech/psl/rpsl376/ergon/ERA5/NETCDF/GLOBAL_025/hourly'
      suf="ap1e5.GLOBAL_025"
     else
      ANA_DIR='/gpfsstore/rech/psl/rpsl376/ergon/ERA5/NETCDF/GLOBAL_025/4xdaily'
      suf="aphe5.GLOBAL_025"
     fi
   elif [ "$rea" = "ERAI" ] ; then
      #ANA_DIR='http://prodn.idris.fr/thredds/dodsC/ipsl_private/rpsl376/ERAI/NETCDF/GLOBAL_1125/4xdaily'
      ANA_DIR='/gpfsstore/rech/psl/rpsl376/ergon/ERAI/NETCDF/GLOBAL_1125/4xdaily'
      #/gpfsstore/rech/psl/rpsl376/ergon/ERAI/NETCDF/GLOBAL_075/4xdaily/AN_PL/2006/u.200601.aphei.GLOBAL_075.nc 
      suf="aphei.GLOBAL_1125"
   else  #OPERA
      #ANA_DIR='http://prodn.idris.fr/thredds/dodsC/ipsl_private/rpsl376/OPERA/NETCDF/4xdaily'
      ANA_DIR='/gpfsstore/rech/psl/rpsl376/ergon/OPERA/NETCDF/4xdaily'
      suf="aph.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