#!/bin/bash

#set -vx

source ~/env_Multi_atlas.sh

echo OBSDIR is defined in env.sh : $OBSDIR

years=$1
SEAS=$2
run=$3
zone=$4

vars='pr rstt rlut rtt crelt crest crett hfns hfls eva hfss tas sst hurs tauu tauv psl zg500 rsts rsut rsutcs rlutcs rsds rsdscs rsus rsuscs rlds rldscs rlus rsdt albs albt cress crels crets rts rah rahcs rahcre rlah rlahcs rlahcre rsah rsahcs rsahcre prw rttcs rsttcs cllcalipso clmcalipso clhcalipso cll clm clh'

if [ "$run" = "OBS" ]; then
  SIM="$run"
else
  SIM="$run"_"$years"
fi

echo ferret is $ferret

force_create=0
GR=VLR
MAINDIR=$DODSDIR/$login/lmdz
WRK=$MAINDIR/WORK/biais$$
orig=$MAINDIR/$GR/$SEAS

if [ ! -d $orig ] ; then mkdir $orig ; fi
if [ ! -f $orig/OBS/NC ] ; then
   ln -s $OBSDIR/$GR/$SEAS/OBS/NC $orig/OBS/.
fi
   

if [ -d $WRK ] ; then WRK=$WRK$$ ; fi
if [ -d $WRK ] ; then echo $WRK existe deja ; exit ; fi

mkdir -p $WRK

cd $orig
for sim in $SIM ; do

   echo $sim
   bad=0

   for var in $vars ; do

     #echo $orig/$sim/$d
     for d in $zone ; do
      echo $orig/$sim/$d
        mkdir -p $orig/$sim/$d
        if [ $force_create = 1 ] ; then
           \rm -f $orig/$sim/$d/$var.nc $orig/$sim/$d/$var
        fi
     done

     cd $WRK
     \rm -rf tmp* ferret*
     varobs=$var
     varmod=$var
     obs=$var
     echo $pal

# contour levels for biases
     case $var in
           pr) lev='(-INF)(0.5)(1,4,1)(6,14,2)(INF)' ; levd='(-Inf)(-5)(-2,2,1)(-0.5)(-0.2)(0.2)(0.5)(5)(Inf)'; pal=rain_cmyk ;;
           prl) lev='(-INF)(0.5)(1,4,1)(6,14,2)(INF)' ; levd='(-Inf)(-5)(-2,2,1)(-0.5)(-0.2)(0.2)(0.5)(5)(Inf)'; pal=rain_cmyk ; varmod=prl ; varobs=86400.*pr ;;
           prc) lev='(-INF)(0.5)(1,4,1)(6,14,2)(INF)' ; levd='(-Inf)(-5)(-2,2,1)(-0.5)(-0.2)(0.2)(0.5)(5)(Inf)'; pal=rain_cmyk ; varmod=prc ; varobs=86400.*pr ;;
           hfns) lev='(-INF)(-200,200,50)(-75)(-25)(25)(75)(INF)'; levd='(-Inf)(-80,80,20)(Inf)' ;;
           tas) lev='(-INF)(-60,20,10)(6,26,4)(24)(27)(28)(30)(INF)' ; levd='(-Inf)(-8,8)(-4,4)(-2,2,1)(Inf)' ; varmod=-273.15+tas ; varobs=-273.15+tas ;;
           sst) lev='(-INF)(-60,20,10)(6,26,4)(24)(27)(28)(30)(INF)' ; levd='(-Inf)(-8,8)(-4,4)(-2,2,1)(Inf)' ; varmod=-273.15+tas ; varobs=-273.15+sst ;;
           tasc) lev='(-INF)(-60,20,10)(6,26,4)(24)(27)(28)(30)(INF)' ; levd='(-Inf)(-8,8)(-4,4)(-2,2,1)(Inf)' ;;
           tauu) lev='(-INF)(-.16,.16,0.02)(INF)' ; levd='(-1.)(-0.16,0.16,0.02)(1.)' ;;
           tauv) lev='(-INF)(-0.16,0.16,0.02)(INF)' ; levd='(-1.)(-0.16,0.16,0.02)(1.)' ;;
           psl) lev='(-INF)(975,1030,5)(1040,1100,10)(INF)' ; levd='(-Inf)(-100,-60,10)(-20,20,2.5)(20,100,10)(Inf)' ;;
           zg500) lev='(-INF)(4900,5800,50)(INF)' ; levd='(-Inf)(-260,100,20)(Inf)' ;;
           sfcWind) levd='(-Inf)(-1.5,1.5,0.2)(Inf)' ;;
           hurs) levd='(-Inf)(-10,10,1)(Inf)' ;;
           ts) levd='(-Inf)(-3,3,0.5)(Inf)' ; varobs=tsk ;;
           deltat) levd='(-Inf)(-1,1,0.1)(Inf)' ; varmod=tsmtas ;;
           huss) levd='(-Inf)(-2.4,2.4,0.3)(Inf)' ; varmod=1000.*huss ;;
           prw) lev='(-Inf)(0.,65.,5)(Inf)' ; levd='(-Inf)(-30.,30.,5)(Inf)' ;;
           hfls) lev='(-Inf)(0,200,20)(Inf)' ; levd='(-Inf)(-50,-10,10)(-5,5,5)(5,50,10)(Inf)' ;;
           eva) lev='(-Inf)(0.5)(1,4,1)(6,14,2)(Inf)' ; levd='(-Inf)(-5)(-2,2,1)(-0.5)(-0.2)(0.2)(0.5)(5)(Inf)' ;;
           hfss) lev='(-Inf)(0,120,10)(Inf)' ; levd='(-Inf)(-50,-10,10)(-5,5,5)(5,50,10)(Inf)' ;;
           hfns) levd='(-Inf)(-50,-10,10)(-5,5,5)(5,50,10)(Inf)' ;;
           bils) levd='(-Inf)(-50,-10,10)(-5,5,5)(5,50,10)(Inf)' ; varobs=bils ;;
           cress) lev='(-Inf)(-120,-10,10)(Inf)' ; levd="$radd" ;;
           crels) lev='(-Inf)(0,100,10)(Inf)' ; levd="$radd" ;;
           crets) lev='(-Inf)(-60,60,10)(Inf)' ; levd="$radd" ;;
           crest) lev='(-Inf)(-120,-10,10)(Inf)' ; levd="$radd" ;;
           crelt) lev='(-Inf)(0,70,5)(Inf)' ; levd="$radd" ;;
           crett) lev='(-Inf)(-60,60,10)(Inf)' ; levd="$radd" ;;
           rts) lev='(-Inf)(-20,200,10)(Inf)' ; levd="$radd" ;;
           rtt) lev='(-Inf)(-100,100,10)(Inf)' ; levd="$radd" ;;
           rttcs) lev='(-Inf)(-100,100,10)(Inf)' ; levd="$radd" ;;
           rlut) lev='(-Inf)(150,320,10)(Inf)' ; levd="$radd" ;;
           rlutcs) lev='(-Inf)(150,350,10)(Inf)' ; levd="$radd" ;;
           rlus) lev='(-Inf)(120,400,20)(440)(480)(Inf)' ; levd="$radd" ;;
           rlds) lev='(-Inf)(100,420,20)(Inf)' ; levd="$radd" ;;
           rldscs) lev='(-Inf)(80,320,20)(Inf)' ; levd="$radd" ;;
           rsds) lev='(-Inf)(80,320,20)(Inf)' ; levd="$radd" ;;
           rsdscs) lev='(-Inf)(80,320,20)(Inf)' ; levd="$radd" ;;
           rsuscs) lev='(-Inf)(20,150,10)(Inf)' ; levd="$radd" ;;
           rsus) lev='(-Inf)(10,150,10)(Inf)' ; levd="$radd" ;;
           rsutcs) lev='(-Inf)(10)(30)(50,160,10)(180)(Inf)' ; levd="$radd" ;;
           rstt) lev='(-Inf)(0,330,10)(Inf)' ; levd="$radd" ;;
           rsttcs) lev='(-Inf)(0,330,10)(Inf)' ; levd="$radd" ;;
           rsts) lev='(-Inf)(0,330,10)(Inf)' ; levd="$radd" ;;
           rstscs) lev='(-Inf)(0,400,10)(Inf)' ; levd="$radd" ;;
           rsut) lev='(-Inf)(50,160,10)(Inf)' ; levd="$radd" ;;
           rsutcs) lev='(-Inf)(50,160,10)(Inf)' ; levd="$radd" ;;
           rsdt) lev='(-Inf)(0,700,10)(Inf)' ; levd="$radd" ;;
           albt|albs) lev='(-Inf)(0,100,5)(Inf)' ; levd='(-Inf)(-50,-10,10)(-5,5,5)(10,50,10)(Inf)' ;;
           cll) lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ; varobs=cllcalipso ;;
           clm) lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ; varobs=clmcalipso ;;
           clh) lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ; varobs=clhcalipso ;;
           clt) lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ; varobs=cltcalipso ;;
           cltcalipso|clhcalipso|clmcalipso|cllcalipso) lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           FO1cllcalipso|FO2cllcalipso|FO3cllcalipso) varmod=cllcalipso ; lev='(1)(2,24,2)(30,90,10)(INF)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           FO1clmcalipso|FO2clmcalipso|FO3clmcalipso) varmod=clmcalipso ; lev='(1)(2,24,2)(30,90,10)(INF)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           FO1clhcalipso|FO2clhcalipso|FO3clhcalipso) varmod=clhcalipso ; lev='(1)(2,24,2)(30,90,10)(INF)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           clcalipso) lev='(1,4,1)(4,32,2)(40,70,10)(Inf)' ; levd='(-Inf)(-30,30,5)(Inf)' ;;
           rah) lev='(-INF)(-180,0,10)(INF)'; levd="$rahd" ;;
           rahcs) lev='(-INF)(-180,0,10)(INF)'; levd="$rahd" ;;
           rahcre) lev='(-INF)(-40,70,5)(INF)' ; levd="$rahd" ;;
           rlah) lev='(-INF)(-280,-50,10)(INF)'; levd="$rahd" ;;
           rlahcs) lev='(-INF)(-280,-50,5)(INF)' ; levd="$rahd" ;;
           rlahcre) lev='(-INF)(-60,60,5)(INF)' ; levd="$rahd" ;;
           rsah) lev='(-INF)(0,150,10)(INF)'; levd="$rahd" ;;
           rsahcs) lev='(-INF)(0,150,10)(INF)' ; levd="$rahd" ;;
           rsahcre) lev='(-INF)(-20,20,5)(INF)' ; levd="$rahd" ;;
           cltoce|clhoce|clmoce|clloce) vlim=0:100 ; lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           cltter|clhter|clmter|cllter) vlim=0:100 ; lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           cltsic|clhsic|clmsic|cllsic) vlim=0:100 ; lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           cltlic|clhlic|clmlic|clllic) vlim=0:100 ; lev='(-Inf)(0,10,2)(0,100,10)(Inf)' ; levd='(-Inf)(-70)(-40,40,5)(70)(Inf)' ;;
           *) echo variable $var non prevue ; exit
     esac
     echo FIGURE FIGURE FIGURE $sim $title $var
     cat <<eod>| open.jnl
        use "$orig/$sim/NC/all.nc"
eod

        cat <<eod>| vardef.jnl
        let cress=rsds-rsus-rsdscs+rsuscs
        let crels=rlds-rldscs
        let rstt=rsdt-rsut
        let rsttcs=rsdt-rsutcs
        let rsts=rsds-rsus
        let rtt=rstt-rlut
        let rttcs=rsttcs-rlutcs
        let albt=100*(rsut/rsdt)
        let albs=100*(rsus/rsds)
        let crets=cress+crels
        let crest=rsutcs-rsut
        let crelt=rlutcs-rlut
        let crett=crest+crelt
        let rts=rsds-rsus+rlds-rlus
        let bil=rts-hfls-hfss
        !let hfns=hfls+hfss+rts
        !let hfns=bils
        let eva=hfls*0.03456
        let tsmtas=ts-tas
        let tsk=ts+273.15
!       let pslhPa=psl/100.
        let tasc=tas-273.15
     
        let rlah=rlus-rlds-rlut
        let rlahcs=rlus-rldscs-rlutcs
        let rlahcre=rlah-rlahcs

        let rsah=rsdt-rsut+rsus-rsds
        let rsahcs=rsdt-rsutcs+rsuscs-rsdscs
        let rsahcre=rsah-rsahcs

        let rah=rsah+rlah
        let rahcs=rsahcs+rlahcs
        let rahcre=rah-rahcs

        let clloce = if oce[d=1] ge 99 then cllcalipso
        let cllter = if ter[d=1] ge 99 then cllcalipso
        let cllsic = if sic[d=1] ge 99 then cllcalipso
        let clllic = if lic[d=1] ge 99 then cllcalipso

        let clmoce = if oce[d=1] ge 99 then clmcalipso
        let clmter = if ter[d=1] ge 99 then clmcalipso
        let clmsic = if sic[d=1] ge 99 then clmcalipso
        let clmlic = if lic[d=1] ge 99 then clmcalipso

        let clhoce = if oce[d=1] ge 99 then clhcalipso
        let clhter = if ter[d=1] ge 99 then clhcalipso
        let clhsic = if sic[d=1] ge 99 then clhcalipso
        let clhlic = if lic[d=1] ge 99 then clhcalipso

        reg/x=-180.:180./l=1/y=-90:90 !globe
eod

cat <<eod>| ZON.jnl
        go open.jnl
        go vardef.jnl
        let vv=$varmod[i=@ave,d=1]
        list/format=(f15.3,f15.3) y*(0*vv+1.),vv
        !list/noh/file=$varmod.ASCII/format=(f15.3,e15.3,e15.3) y*(0*vv[d=1]+1.),vv[d=1],vv[d=2]
eod

cat <<eod>| AMMACROSS.jnl
        go open.jnl
        go vardef.jnl
        let vv=$varmod[i=@ave,d=1]
        list/x=-10.:10./y=-5.:25./format=(f15.3,f15.3) y*(0*vv+1.),vv
        !list/noh/file=$varmod.ASCII/format=(f15.3,e15.3,e15.3) y*(0*vv[d=1]+1.),vv[d=1],vv[d=2]
eod

cat <<eod>| GLOB.jnl
        go open.jnl
        go vardef.jnl
        list/format=(f15.3) $varmod[i=@ave,j=@ave,d=1]
eod


        for d in $zone ; do
           outf=$orig/$sim/$d/$var
           echo $outf
           soutf=`wc -l ${outf} |awk '{print $1}'`
           if [ $force_create = 1 ] || [ "${soutf}" == "0" ]; then rm -f $outf ; fi
           if [ ! -f $outf ] ; then
               $ferret -batch tmp.ps -script $d.jnl | sed -e '/:/d' -e '/\*/d' >| $outf
           fi
        done

   done


done