source: BOL/Multi_atlas/make_axe2.sh @ 5304

#!/bin/bash -vx
# Par exemple:
# ./make_axe2.sh V6.5AIS8
# -----------------------------------------------------------------

local=`pwd`
#source ~/env_Multi_atlas.sh
#module purge
#module load gcc/9.4.0
#module load netcdf-fortran/4.5.3-parallel-netcdf
#module load cdo/2.0.6
#module load nco/5.0.1
#module load ncl/6.6.2
#module load anaconda3-py/2020.11
#module load ferret/7.6.0

source ${local}/env_Multi_atlas_axe2.sh
set -vx
pr_yr=1
pr_da=1

# -----------------------------------------------------------------
# Based on /home/fabric/LMDZ/UTILS/config.sh
# -----------------------------------------------------------------

runstxt=~/${local}/runs.txt

while test -n "${1}"; do
   case $1 in
      -runstxt) runstxt=$2 ; shift ;;
      -h) echo "Example: ./make_axe2.sh [-runstxt runs.description.file] V6.5AIS8" ;;
      *) comp=$1 ;;
   esac
   shift
done


WRK=$MULTIDIR/AXE2/WORK$$
if [ -d $WRK ] ; then WRK=$WRK$$ ; fi
mkdir -p $WRK ; cd $WRK

COMP_D=$MULTIDIR/$comp
if [ ! -d $COMP_D ] ; then echo "$COMP_D not found" && exit ; fi

# -----------------------------------------------------------------
# Liste des simulations de la comparasion $comp
# -----------------------------------------------------------------
liste_simy=""
DEF_FILE=$COMP_D/def.txt
for s in `awk ' {print $1} ' $DEF_FILE` ; do liste_sim="$liste_sim $s" ; done
for s in `awk ' {print $1"_"$2} ' $DEF_FILE` ; do liste_simy="$liste_simy $s" ; done

echo liste_sim: $liste_sim
echo liste_simy: $liste_simy
# -----------------------------------------------------------------
# -----------------------------------------------------------------
# Analyse de la distribution longitudinale des precipitation
# sur ocean et continent
# -----------------------------------------------------------------
# -----------------------------------------------------------------

echo OK0

outputdir=$COMP_D/AXE2/PR_YR
if [ -d $outputdir/PNG ] ; then pr_yr=0 ; fi
if [ $pr_yr = 1 -o 0 = 0 ] ; then

echo OK1

outputdir=$COMP_D/AXE2/PR_YR
for sub in [ PNG NC TXT ] ; do mkdir -p $outputdir/$sub ; done

echo OK2 $liste_sim

listey=""
for simy in $liste_simy ; do
  sim=`echo $simy | sed -e 's/_[0-9][0-9][0-9][0-9]_[0-9][0-9][0-9][0-9]//'`
  years=`echo $simy | sed -e 's/'$sim'_//'`
  simdir=`grep -w "^$sim "  $runstxt | awk ' { print $2 } '`
  ncks -v precip $simdir/ATM/Analyse/SE/${sim}_SE_${years}_1M_histmth.nc -O grid.nc
  cdo timavg -remapcon,grid.nc ${MAIN_CMOR}/OBS/pr.nc gpcp.nc
  cdo timavg -remapcon,grid.nc ${MAIN}/TRMM/prm2000-2009.nc trmm.nc

######################################################
# Sauvegarde d'un masque oceanique pour le diag MJO
######################################################
   if [ ! -f $WRK/pourc_oce.nc ] ; then ncks -d time_counter,1 -v pourc_oce  "$simdir/ATM/Analyse/SE/${sim}_SE_${years}_1M_histmth.nc" $WRK/pourc_oce.nc ; fi
   echo ncks -d time_counter,1 -v pourc_oce  "$simdir/ATM/Analyse/SE/${sim}_SE_${years}_1M_histmth.nc" $WRK/pourc_oce.nc
cat <<eod>| tmpyr.jnl
use "$simdir/ATM/Analyse/SE/${sim}_SE_${years}_1M_histmth.nc"
use gpcp.nc
use trmm.nc
let con=pourc_ter[d=1,l=1:12@ave]+pourc_lic[d=1,l=1:12@ave]
let oce=pourc_oce[d=1,l=1:12@ave]+pourc_sic[d=1,l=1:12@ave]
let ocem=if ( oce[l=1] ge 99. ) then 1
let conm=if ( con[l=1] ge 99. ) then 1
let mixm=if ( oce[l=1] gt 1. AND con[l=1] gt 1. )  then 1
let totm=tsol[d=1,l=1]*0.+1.

let filtre=var
go tmp0.jnl TOT,8,totm
frame/file=tot.gif
let filtre=if ( oce[l=1] ge 99.9 ) then var*oce/100.
go tmp0.jnl OCE,9,ocem
frame/file=oce.gif
let filtre=if ( con[l=1] ge 99.9 ) then var*con/100.
go tmp0.jnl CONT,10,conm
frame/file=con.gif
let filtre=if ( con[l=1] gt 0.1 AND oce[l=1] gt 0.1 ) then var
go tmp0.jnl MIXT,11,mixm
frame/file=mix.gif
list/x=-10:10/y=12:18 "PRlmdzAMMA",86400*precip[d=1,i=@ave,j=@ave,l=1:12@ave]
list/x=-10:10/y=12:18 "PRtrmmAMMA",r[d=3,i=@ave,j=@ave,l=1:12@ave]
list/x=-10:10/y=12:18 "PRgpcpAMMA",86400*pr[d=2,i=@ave,j=@ave,l=1:12@ave]
eod

cat <<eod>| tmp0.jnl
let mask=\$3
let var=86400*precip[d=1]*mask[i=@ave,j=@ave]
plot/line=\$2/vlim=0:10/title="$sim $years pr \$1 (mm/d)" filtre[i=@ave,l=1:12@ave]
list/y=-50:50 "PRlmdz\$1",filtre[i=@ave,j=@ave,l=1:12@ave]
let var=r[d=3]*mask
plot/line=1/o/title="TRMM" filtre[i=@ave,l=1]
list/y=-50:50 "PRtrmm\$1",filtre[i=@ave,j=@ave,l=1:12@ave]
let var=86400*pr[d=2]*mask
plot/line=7/o/title="GPCP" filtre[i=@ave,l=1:12@ave]
list/y=-50:50 "PRgpcp\$1",filtre[i=@ave,j=@ave,l=1:12@ave]
quit
eod

\rm gotmpyr.jnl
cat<<eof>gotmpyr.jnl
go tmpyr.jnl
quit
eof
echo ferret<gotmpyr.jnl   | grep "^I.*.PR" | awk ' { print $4 , $5 } ' > $outputdir/TXT/${sim}_${years}
ferret<gotmpyr.jnl | grep "^I.*.PR" | awk ' { print $4 , $5 } ' > $outputdir/TXT/${sim}_${years}
for type in con mix tot oce ; do
convert -density 144 $type.gif $outputdir/PNG/${sim}_${years}_$type.png
listey="$listey ${sim}_${years}_$type"
done

done

if [ -f $COMP_D/entete.html ] ; then cat $COMP_D/entete.html > $outputdir/index.html ; fi
~/Multi_atlas/concat_html.sh $COMP_D/AXE2/PR_YR/PNG "OK" "$listey" 4 >> $outputdir/index.html


fi # pr_yr=1

# -----------------------------------------------------------------
# -----------------------------------------------------------------
# Analyse de la variabilité jour à jour des pluies
# -----------------------------------------------------------------
# -----------------------------------------------------------------
outputdir=$COMP_D/AXE2/PR_DAY
# if [ -d $outputdir ] ; then pr_da = 0 ; fi

if [ $pr_da = 1 ] ; then

listey=""

outputdir=$COMP_D/AXE2/PR_DAY
for sub in [ PNG NC TXT ] ; do mkdir -p $outputdir/$sub ; done

# PLOT SIMULATION RESULTS
first=1
for simy in $liste_simy ; do
  sim=`echo $simy | sed -e 's/_[0-9][0-9][0-9][0-9]_[0-9][0-9][0-9][0-9]//'`
  years=`echo $simy | sed -e 's/'$sim'_//'`

  yi=`echo $years | cut -d_ -f1`
  yf=`echo $years | cut -d_ -f2`
  echo FICHIERS POUR LA SIMULATION $sim annees $years
  simdir=`grep -w "^$sim "  $runstxt | awk ' { print $2 } '`
  run=`basename $simdir `
  prd=`ls $simdir/ATM/Analyse/TS_DA/${run}_${yi}*_${yf}*_1D_precip.nc | tail -1`
  prc=`ls $simdir/ATM/Analyse/TS_DA/${run}_${yi}*_${yf}*_1D_pluc.nc | tail -1`
  prl=`ls $simdir/ATM/Analyse/TS_DA/${run}_${yi}*_${yf}*_1D_plul.nc | tail -1`
  ln -s $prd $outputdir/NC/${run}_${yi}_${yf}.nc
  prm=${MAIN_SE}/ORIG/${run}_SE_${years}_1M_histmth.nc
  if [ ! -f ${prm} ]; then
     prm=$simdir/ATM/Analyse/SE/${run}_SE_${years}_1M_histmth.nc
  fi
  if [ $first = 1 ] ; then
     remap=$prm ; first=0
     trmm=$outputdir/NC/TRMM.nc
     if [ ! -f $trmm ] ; then
        ncks -v precip $prm -O $outputdir/NC/grid.nc
        ###cdo remapcon,$outputdir/NC/grid.nc /thredds/ipsl/fabric/lmdz/TRMM/pr2000-2009.nc $trmm
        cdo remapcon,$outputdir/NC/grid.nc ${MAIN}/TRMM/pr2000-2009.nc $trmm
     fi
     ts_da="TRMM"
  fi
  echo PRD $prd
  if [ -f "$prd" ] ; then  ts_da="$ts_da ${run}_${yi}_${yf}" ; fi
done



for sim in $ts_da ; do
echo Tracer de $sim
if [ $sim = TRMM ] ; then
    var=r
else
    var='86400*precip'
fi

nom=`echo $sim | sed -e 's/_/ /g'`
cat <<eod>| tmp.jnl
use "$WRK/pourc_oce.nc"
use "$outputdir/NC/$sim.nc"
DEFINE VIEWPORT/XLIM=0.,1./YLIM= .5,1. V1
DEFINE VIEWPORT/XLIM=0.,1./YLIM= 0.2,0.7 V2
DEFINE VIEWPORT/XLIM=0.,1./YLIM= 0.,0.35 V3
set memory/size=1000
let var=$var
let dv=var-var[l=@sbx:30]
let dv2=dv*dv
let dv2oce=if ( pourc_oce[d=1,l=1] ge 99. ) then dv2
reg/y=-50:50
fill/pal=rain_cmyk/nolab/set/lev=(3,15,3)(Inf) (dv2[l=@ave])^0.5
list "PRd1-30",(dv2[i=@ave,j=@ave,l=@ave])^0.5
! ppl axlsze,0.14,0.14
ppl fill
label/nouser 3.,6.7,0,0,0.24 "$nom"
label/nouser 3.,6.2,0,0,0.24 "Std deviation, pr(1d)-pr(30d)"
go land
frame/file=hf.gif
let dv=var[l=@sbx:20]-var[l=@sbx:120]
fill/pal=rain_cmyk/nolab/set/lev=(1,5,1)(Inf) (dv2[l=@ave])^0.5
list "PRd20-120",(dv2[i=@ave,j=@ave,l=@ave])^0.5
list/x=50:180/y=-20:20 "PRd20-120MJO",(dv2oce[i=@ave,j=@ave,l=@ave])^0.5
ppl fill
go land
label/nouser 3.,6.7,0,0,0.24 "$nom"
label/nouser 3,6.2,0,0,0.24 "Std deviation, pr(20d)-pr(120d)"
frame/file=is.gif
let pr1= if ( $var gt 1 ) then 1 else 0
fill/pal=rain_cmyk/nolab/set/lev=(0.1,1,0.15) pr1[l=@ave]
list "PRnd1",pr1[i=@ave,j=@ave,l=@ave]
ppl fill
label/nouser 3.,6.7,0,0,0.24 "$nom"
label/nouser 3.,6.2,0,0,0.24 "time frac. wth daily pr. above 1mm/day"
go land
frame/file=nd.gif
plot/line=7/vlim=0:9/title="$nom, PR (mm/day)" ${var}[i=@ave,l=@ave]
let pr10= if ( $var gt 10 ) then $var else 0
let pr20= if ( $var gt 20 ) then $var else 0
let pr50= if ( $var gt 50 ) then $var else 0
let pr100= if ( $var gt 100 ) then $var else 0
let pr200= if ( $var gt 200 ) then $var else 0
plot/o/line=8/title="PR>10mm/d" pr10[i=@ave,l=@ave]
plot/o/line=9/title="PR>20mm/d" pr20[i=@ave,l=@ave]
plot/o/line=10/title="PR>50mm/d" pr50[i=@ave,l=@ave]
plot/o/line=11/title="PR>100mm/d" pr100[i=@ave,l=@ave]
plot/o/line=12/title="PR>200mm/d" pr200[i=@ave,l=@ave]
frame/file=zon.gif
let dpr10=$var-pr10
let dpr20=pr10-pr20
let dpr50=pr20-pr50
! list "PRTOT",${var}[i=@ave,j=@ave,l=@ave]
list "PR0-10",dpr10[i=@ave,j=@ave,l=@ave]
list "PR10-20",dpr20[i=@ave,j=@ave,l=@ave]
list "PR20-50",dpr50[i=@ave,j=@ave,l=@ave]
list "PR50",pr50[i=@ave,j=@ave,l=@ave]
quit
eod

\rm gotmp.jnl
cat<<eof>gotmp.jnl
go tmp.jnl
quit
eof
ferret<gotmp.jnl  | grep "^I.*.PR" | awk ' { print $4 , $5 } ' > $outputdir/TXT/$sim

for type in hf is nd zon ; do
convert $type.gif $outputdir/PNG/${sim}_$type.png
listey="$listey ${sim}_$type"
done
done

echo $ts_da
if [ -f $COMP_D/entete.html ] ; then cat $COMP_D/entete.html > $outputdir/index.html ; fi

~/Multi_atlas/concat_html.sh $COMP_D/AXE2/PR_DAY/PNG "OK" "$listey" 4 >> $outputdir/index.html

fi # pr_da=1


#######################################################################
# Bar chart bilan
#######################################################################

cd $COMP_D/AXE2
liste_simy=""
DEF_FILE=$COMP_D/def.txt
for s in `awk ' {print $1"_"$2 } ' $DEF_FILE` ; do liste_simy="$liste_simy $s" ; done
echo $liste_simy


ref=1
mkdir -p XMGR
for sim in $liste_simy ; do
   if [  $ref = 1 ] ; then
      clims="lmdz gpcp trmm"
   else
      clims=lmdz
   fi
   ref=0
   for clim in $clims ; do \rm -f $clim ; for type in TOT OCE CONT MIXT AMMA  ; do grep $type PR_YR/TXT/$sim | grep $clim | awk ' { print $2 } '  >> $clim ; done ; done
   for clim in $clims ; do for type in TOT OCE CONT MIXT AMMA  ; do grep $type PR_YR/TXT/$sim | grep $clim | awk ' { print $2 } '   ; done ; done
    \mv -f lmdz XMGR/$sim
   awk ' { print $2 } ' PR_DAY/TXT/$sim >> XMGR/$sim
done

mv trmm XMGR/TRMM
mv gpcp XMGR/GPCP
awk ' { print $2     } ' PR_DAY/TXT/TRMM >> XMGR/TRMM
awk ' { print $2 * 0 } ' PR_DAY/TXT/TRMM >> XMGR/GPCP
cd XMGR



nsims=`echo $liste_simy | wc -w | awk ' { print  $1 + 2  } '`
size=`echo $nsims | awk ' { print 3.5 / $1 } '`

cat <<eod>| precip.param
g0 type Chart
    world -0.5, 0, 12.5, 8


    yaxis  label "PR (mm/day) "
    yaxis  tick major 1


    xaxis  ticklabel font 0
    xaxis  ticklabel color 1
    xaxis  tick place both
    xaxis  tick spec type both
    xaxis  tick spec 13
    xaxis  tick major 0, 0
    xaxis  ticklabel 0, "TOT"
    xaxis  tick major 1, 1
    xaxis  ticklabel 1, "OCE"
    xaxis  tick major 2, 2
    xaxis  ticklabel 2, "CONT"
    xaxis  tick major 3, 3
    xaxis  ticklabel 3, "MIXT"
    xaxis  tick major 4, 4
    xaxis  ticklabel 4, "AMMA"
    xaxis  tick major 5, 5
    xaxis  ticklabel 5, "SIG1-30"
    xaxis  tick major 6, 6
    xaxis  ticklabel 6, "SIG20-120"
    xaxis  tick major 7, 7
    xaxis  ticklabel 7, "MJO"
    xaxis  tick major 8, 8
    xaxis  ticklabel 8, "FD(PR>1)"
    xaxis  tick major 9, 9
    xaxis  ticklabel 9, "<10"
    xaxis  tick major 10, 10
    xaxis  ticklabel 10, "10-20"
    xaxis  tick major 11, 11
    xaxis  ticklabel 11, "20-50"
    xaxis  tick major 12, 12
    xaxis  ticklabel 12, ">50"
    xaxis  tick major 13, 13
    xaxis  ticklabel char size 1.3
    xaxis  ticklabel angle 90

    legend 0.8, 0.85
eod

color=( 1 1 1 2 3 4 5 6 8 9 10 11 12 13 14 15 )
pattern=( 10 6 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 )

is=0
while [ $is != $nsims ] ; do
cat <<eod>> precip.param
s$is type bar
s$is symbol fill pattern ${pattern[$is]}
s$is symbol color ${color[$is]}
s$is symbol fill color ${color[$is]}
s$is symbol size $size
s$is line linestyle 0
eod
is=$(( $is + 1 ))
done

xmgrace TRMM GPCP $liste_simy -param precip.param -legend load -hardcopy -hdevice EPS -printfile tmp.eps
epstopdf tmp.eps
convert -density 144 tmp.pdf tmp.png