#!/bin/sh
# $Id$
#
# Borrador para el ploteo de mapa de superficie
#
# Uso
# mapa_sfc -Jarg -Rarg fecha pasada proy ind [R|P]
#
# si se incluye el argumento R indica un mapa regional (mayor definicion)
# si se incluye el argumento P indica un mapa polar (menos meridianos)
PATH="${PATH}:/usr/local/bin"

if test ${#} -lt 6 ; then
  echo "Uso:"
  echo "mapa_sfc -Rarg -Jarg fecha pasada proy ind [R|P]"
  exit
fi

VIENTO="No"

# Mira el intervalo de presiones
if test ${#} -eq 7 ; then
  if test ${7} = "R" ; then
    #caso regional
    VIENTO="Si"
    REG=yes
    A=200
    dp=2
    dpl=4
    dm=5
    s1=10
    s2=9
    dx=2
    dy=1.5
  elif test ${7} = "P" ; then
    REG=no
    A=800
    dp=4
    dpl=8
    dm=15
    s1=10
    s2=9
    dx=10
    dy=5
  else
    REG=no
    A=800
    dp=4
    dpl=8
    dm=10
    s1=10
    s2=9
    dx=4
    dy=3
  fi
else
  REG=no
  A=800
  dp=4
  dpl=8
  dm=10
  s1=9
  s2=9
  dx=4
  dy=3
fi 

##### Preparacion de datos #######
#FECHA=`date -u +%Y%m%d`
FECHA=${3}
PASADA=${4}
PROY=${5}
LEVEL=SFC
LEVEL2=10M
if test ${5} -lt 100 ; then
 PROY0=0${5}
else
 PROY0=${5}
fi
PREFIX=${FECHA}${PASADA}H${PROY0}_${6}_
if test "${REG}" = "yes" ; then 
  PREFIX0=reg${FECHA}${PASADA}H${PROY}_
  tempo=reg
else
  PREFIX0=gen${FECHA}${PASADA}H${PROY}_
  tempo=gen
fi
PREFIX2=vie${FECHA}${PASADA}H${PROY}_
tempo2=vie
DIR_BASE=/home/gbv/gfs/latest_gfs_${PASADA}
TARGET=${PREFIX}${LEVEL}.ps
TARGET_PNG=${PREFIX}${LEVEL}.png
TARGET_GIF=${PREFIX}${LEVEL}.gif
TARGET_C=${PREFIX}${LEVEL}

# Hora_nivel_var
HNV_PRMSL=${PASADA}_MSL_PRMSL
HNV_PREC=${PASADA}_SFC_APCP
HNV_NUB=${PASADA}_TOT_TCDC
HNV_1000_HGT=${PASADA}_1000_HGT
HNV_500_HGT=${PASADA}_500_HGT
HNV_10M_U=${PASADA}_${LEVEL2}_UGRD
HNV_10M_V=${PASADA}_${LEVEL2}_VGRD

#CPT
DIR_CPT=/home/gbv/gfs/share
PREC_CPT=${DIR_CPT}/prec.cpt
NUB_CPT=${DIR_CPT}/nub.cpt
FRIO_CPT=${DIR_CPT}/frio.cpt
CALOR_CPT=${DIR_CPT}/calor.cpt


##geografia
#GEO="-Rg -JG${2}/${1}/12"
GEO="${1} ${2}" 

#limite de viento para pintar 
vlim=15.0

#precision en el resampleo de los grids
if test ${REG} = "yes" ; then
  resam=6m
else
  resam=15m
fi

### Ejemplos de proyeeciones
#GEO="-R-180/180/60/90 -JG0/90/16c 00 36" <- Polo norte           
#GEO="-R-40/15/70/70r -JG0/40/16c"  <- Templado norte
#GEO="-R-40/-35/40/35r -JG0/0/16c 00 36" <- Ecuador
#GEO="-R-95/-55/20/-20r -JG0/-40/16c 00 36" <- Templado Sur
#GEO="-R-180/180/-90/-60 -JG0/-90/16c 00 36" <- Polo Sur        
#GEO="-R-25/40/15/64r -JG-5/50/17c" <- United Kingdom 
#GEO="-R-20/25/20/55r -JG0/40/17c" <- Spain

# extrae los datos del grib y los pasa a binario.
# despues los pasa a grd y los resamplea a mayor precision
#### Presion al nivel del mar ##########
if test ! -e ${PREFIX0}prmsl.grd ; then 
  if test ! -e ${DIR_BASE}/gfs05_${FECHA}${HNV_PRMSL}_f${PROY}.grb ; then 
    echo " no existe ${DIR_BASE}/gfs05_${FECHA}${HNV_PRMSL}_f${PROY}.grb"
    exit
  fi
  wgrib ${DIR_BASE}/gfs05_${FECHA}${HNV_PRMSL}_f${PROY}.grb -d all -bin -nh \
    -o tmp_${tempo}.b
  xyz2grd tmp_${tempo}.b -Gtmp_prsml0_${tempo}.grd -I0.5 -R0/360/-90/90 \
    -Dgrad/grad/"Hp"/=/=/"Presion al Nivel del mar" -ZTLfx  
  # cambia de Pa a HPa
  grdmath tmp_prsml0_${tempo}.grd 0.01 MUL = tmp_prmsl0_${tempo}.grd
  # resample
  grdsample tmp_prmsl0_${tempo}.grd -I${resam} -Q -G${PREFIX0}prmsl.grd -Lx 
fi

#### Precipitacion ##########
if test ! -e ${PREFIX0}prec.grd ; then
  if test ! -e ${DIR_BASE}/gfs05_${FECHA}${HNV_PREC}_f${PROY}.grb ; then 
    exit
  fi
  wgrib ${DIR_BASE}/gfs05_${FECHA}${HNV_PREC}_f${PROY}.grb -d all -bin -nh \
    -o tmp1_${tempo}.b
  xyz2grd tmp1_${tempo}.b -Gtmp_prec0_${tempo}.grd -I0.5 -R0/360/-90/90 \
    -Dgrad/grad/"mm"/=/=/"Precipitacion en 6 horas" -ZTLfx  
  # resample
  grdsample tmp_prec0_${tempo}.grd -I${resam} -Q -G${PREFIX0}prec.grd -Lx 
  # primero selecciona la zona apreciable
  grdclip ${PREFIX0}prec.grd -G${PREFIX0}prec.grd -Sb0.22/NaN 
fi

#### Nubosidad total ##########
if test ! -e ${PREFIX0}nub.grd ; then
  if test ! -e ${DIR_BASE}/gfs05_${FECHA}${HNV_NUB}_f${PROY}.grb ; then 
    exit
  fi
  wgrib ${DIR_BASE}/gfs05_${FECHA}${HNV_NUB}_f${PROY}.grb -d all -bin -nh \
    -o tmp2_${tempo}.b
  xyz2grd tmp2_${tempo}.b -Gtmp_nub0_${tempo}.grd -I0.5 -R0/360/-90/90 \
     -Dgrad/grad/"%"/=/=/"Nubosidad total" -ZTLfx  
  # resample
  grdsample tmp_nub0_${tempo}.grd -I${resam} -Q -G${PREFIX0}nub.grd -Lx 
  # zona despejada se hace transparente
  grdclip ${PREFIX0}nub.grd -G${PREFIX0}nub.grd -Sb50/NaN 
fi

#### Espesor 500-1000 #########
if test ! -e ${PREFIX0}esp_500_1000.grd ; then
  if test ! -e ${DIR_BASE}/gfs05_${FECHA}${HNV_1000_HGT}_f${PROY}.grb ; then 
    exit
  fi
  wgrib ${DIR_BASE}/gfs05_${FECHA}${HNV_1000_HGT}_f${PROY}.grb -d all -bin -nh \
    -o tmp3_${tempo}.b
  if test ! -e ${DIR_BASE}/gfs05_${FECHA}${HNV_500_HGT}_f${PROY}.grb ; then 
    exit
  fi
  wgrib ${DIR_BASE}/gfs05_${FECHA}${HNV_500_HGT}_f${PROY}.grb -d all -bin -nh \
    -o tmp4_${tempo}.b
  xyz2grd tmp3_${tempo}.b -Gtmp5_1000_${tempo}.grd -I0.5 -R0/360/-90/90 \
    -Dgrad/grad/"mgp"/=/=/"Geopotencial 1000 mb" -ZTLfx  
  xyz2grd tmp4_${tempo}.b -Gtmp6_500_${tempo}.grd -I0.5 -R0/360/-90/90 \
    -Dgrad/grad/"mgp"/=/=/"Geopotencial 500 mb" -ZTLfx  
  #calcula el espesor 500-1000 y pasa a Dm
  grdmath tmp6_500_${tempo}.grd tmp5_1000_${tempo}.grd - 0.1 MUL = tmp7_${tempo}.grd
  # resample
  grdsample tmp7_${tempo}.grd -I${resam} -Q -G${PREFIX0}esp_500_1000.grd -Lx 
fi


if test "${VIENTO}" = "Si" ; then
 
#### Componente U del viento a 10 M ##########
if test ! -e ${PREFIX2}${LEVEL2}_U.grd ; then 
  if test ! -e ${DIR_BASE}/gfs05_${FECHA}${HNV_10M_U}_f${PROY}.grb ; then 
    exit
  fi
  wgrib ${DIR_BASE}/gfs05_${FECHA}${HNV_10M_U}_f${PROY}.grb -d all -bin -nh \
    -o tmp_${tempo2}.b
  xyz2grd tmp_${tempo2}.b -Gtmp_${LEVEL2}_tmp0_${tempo2}.grd -I0.5 -R0/360/-90/90 \
    -Dgrad/grad/"kt"/=/=/"Wind U (kt)" -ZTLfx  

  # convert to knots
  grdmath tmp_${LEVEL2}_tmp0_${tempo2}.grd 1.94 MUL = tmp_${LEVEL2}_tmp0_${tempo2}.grd

  # resample
  grdsample tmp_${LEVEL2}_tmp0_${tempo2}.grd -I${resam} -Q -G${PREFIX2}${LEVEL2}_U.grd -Lx 
fi

#### Componente V del viento a 10 M ##########
if test ! -e ${PREFIX2}${LEVEL2}_V.grd ; then 
  if test ! -e ${DIR_BASE}/gfs05_${FECHA}${HNV_10M_V}_f${PROY}.grb ; then 
    exit
  fi
  wgrib ${DIR_BASE}/gfs05_${FECHA}${HNV_10M_V}_f${PROY}.grb -d all -bin -nh \
    -o tmp1_${tempo2}.b
  xyz2grd tmp1_${tempo2}.b -Gtmp_${LEVEL2}_tmp2_${tempo2}.grd -I0.5 -R0/360/-90/90 \
    -Dgrad/grad/"kt"/=/=/"Wind V (Kt)" -ZTLfx

  # convert to knots
  grdmath tmp_${LEVEL2}_tmp2_${tempo2}.grd 1.94 MUL = tmp_${LEVEL2}_tmp2_${tempo2}.grd
  
  # resample
  grdsample tmp_${LEVEL2}_tmp2_${tempo2}.grd -I${resam} -Q -G${PREFIX2}${LEVEL2}_V.grd -Lx 
fi

#### Velocidad , para isotacas #####
if test ! -e ${PREFIX2}${LEVEL2}_vel.grd ; then 
  grdmath ${PREFIX2}${LEVEL2}_U.grd ${PREFIX2}${LEVEL2}_V.grd R2 SQRT = \
  ${PREFIX2}${LEVEL2}_vel.grd   

  # primero selecciona la zona apreciable
  grdclip ${PREFIX2}${LEVEL2}_vel.grd -G${PREFIX2}${LEVEL2}_vclip.grd -Sb${vlim}/NaN 

  # hace un grid con 0.0 y 1.0 dependiendo del limite 
  grdclip ${PREFIX2}${LEVEL2}_vel.grd -G${PREFIX2}${LEVEL2}_vclip0.grd -Sa${vlim}/1.0 -Sb${vlim}/NaN

  # multiplica el campo U  para no pintar flechas menores que el limite
  grdmath ${PREFIX2}${LEVEL2}_vclip0.grd ${PREFIX2}${LEVEL2}_U.grd MUL = ${PREFIX2}${LEVEL2}_U.grd

fi
fi


##### Comandos GMT #######
gmtset GRID_CROSS_SIZE 0 ANNOT_FONT_PRIMARY 0 ANNOT_FONT_SIZE_PRIMARY 9 \
 LABEL_FONT_SIZE 14 TICK_LENGTH -0.1c 

#### Comienza el ploteo 

#### continentes
pscoast ${GEO} -A${A} -N1 -W -G230/190/170 -S170/190/230 -K  > ${TARGET}

#### Contorneo (colores) de nubosidad
grdimage ${PREFIX0}nub.grd -Q -C${NUB_CPT} ${GEO} -K -O >> ${TARGET}

#### Precipitacion
# Ahora pinta la prec
grdimage ${PREFIX0}prec.grd -Q -C${PREC_CPT} ${GEO} -K -O >> ${TARGET}


if test "${VIENTO}" = "Si" ; then
#### pinta las flechas de viento
grdvector ${PREFIX2}${LEVEL2}_U.grd ${PREFIX2}${LEVEL2}_V.grd \
   -G150/120/90 ${GEO} -I${dx}/${dy} \
   -Q0.04c/0.1c/0.1c -E -S50c -T -K -O >> ${TARGET}
fi

#### Costas y fronteras solamente, de nuevo
pscoast ${GEO} -A${A} -N1 -W2 -K -O >> ${TARGET}

#### Mapa Base
psbasemap ${GEO} -B${dm}g${dm}SWEN -K -O >> ${TARGET}

#### Espesores
# por debajo de 552 Dm
grdcontour ${PREFIX0}esp_500_1000.grd ${GEO} -C${FRIO_CPT} \
-A6+g220/220/255+a0+p+kblue+s${s1}  -W+1p,- -K -O >> ${TARGET}

# por encima de 552 Dm
grdcontour ${PREFIX0}esp_500_1000.grd ${GEO} -C${CALOR_CPT} \
-A6+g255/255/200+a0+p+kred+s${s1}  -W+1p,- -K -O >> ${TARGET}

# Espesor de 552 Dm
grdcontour ${PREFIX0}esp_500_1000.grd ${GEO} -C552 \
-A552+g255/255/150+a0+p+kblue+s${s1} -W1p,blue -K -O >> ${TARGET}

if test "${VIENTO}" = "Si" ; then
grdcontour ${PREFIX2}${LEVEL2}_vel.grd ${GEO} -C${vlim} \
 -A${vlim}+p+a0+s${s2}+ukt -T1.0c/0.02c -W0.5p,-. -L${vlim}/100.0 \
   -K -O >> ${TARGET}
fi

#### contorneo de presiones
grdcontour ${PREFIX0}prmsl.grd ${GEO} -C${dp} \
 -A${dpl}+gwhite+p+a0+s${s2} -T1.0c/0.02c:LH -W0.5p \
   -K -O >> ${TARGET}

#### escalas 
if test "${VIENTO}" = "Si" ; then
  psscale -D17.9c/13.5c/6c/0.4c -C${NUB_CPT} -B:"CLOUD COVER":/:"%": -K -O >> ${TARGET}
  psscale -D17.9c/5c/8c/0.4c -C${PREC_CPT} -L -B:"PRECIPITATION":/:"mm": -O >> ${TARGET}
else
  psscale -D17.9c/13c/7c/0.4c -C${NUB_CPT} -B:"CLOUD COVER":/:"%": -K -O >> ${TARGET}
  psscale -D17.9c/4c/8c/0.4c -C${PREC_CPT} -L -B:"PRECIPITATION":/:"mm": -O >> ${TARGET}
fi

rm -f .gmtcommands4 .gmtdefaults4 tmp*_${tempo}* tmp*_${tempo2}* 

#### Conversion a Imagen 
if test "${VIENTO}" = "Si" ; then
  CONV2=/usr/local/bin/conv2v
else
  CONV2=/usr/local/bin/conv2
fi

convert -page A4 ${TARGET} -rotate 90 ${TARGET_GIF}
${CONV2} -i ${TARGET_C} -x 650 -y 595 --gif
rm ${TARGET}