#! /bin/bash -x

#PBS -o /lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img/run_plot_map-smna.out
#PBS -e /lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img/run_plot_map-smna.err
#PBS -l walltime=72:00:00
#PBS -l select=1:ncpus=1
#PBS -A CPTEC
#PBS -V
#PBS -S /bin/bash
#PBS -N PlotFigsSMNA
#PBS -q pesq
#PBS -J 0-8

module swap PrgEnv-cray/6.0.4 PrgEnv-gnu
module unload craype-x86-skylake

datai=#DATAI#

hours=$((${PBS_ARRAY_INDEX} * 24))
data=$(date -u +%Y%m%d%H -d "${datai:0:8} ${datai:8:2} +${hours} hours")  

bpath=/lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img

cd ${bpath}

#VarsAnl=("Topography" 
#         "Land sea mask"
#         "Surface pressure"
#         "Zonal wind (u)"
#         "Meridional wind (v)"
#         "Omega"
#         "Stream function"
#         "Velocity potential"
#         "Geopotential height"
#         "Pressure reduced to msl"
#         "Absolute temperature"
#         "Specific humidity"
#         "Inst. precipitable water"
#         "Surface temperature"
#         "var161 undefined"
#         "var199 undefined"  
#         "10 metre u-wind component"
#         "10 metre v-wind component" 
#         "Vertical dist total cloud cover"
#         "Time mean surface relative humidity")
VarsAnl=("Surface pressure"
         "Zonal wind (u)"
         "Meridional wind (v)"
         "Omega"
         "Stream function"
         "Velocity potential"
         "Geopotential height"
         "Pressure reduced to msl"
         "Absolute temperature"
         "Specific humidity"
         "Inst. precipitable water"
         "Surface temperature"
         "var161 undefined"
         "var199 undefined"  
         "10 metre u-wind component"
         "10 metre v-wind component" 
         "Vertical dist total cloud cover"
         "Time mean surface relative humidity")

Fcts=$(seq 0 24 264)

for fct in ${Fcts[@]}
do

  if [ ${fct} -eq 0 ]; then type="icn"; else type="fct"; fi

  datafct=$(date -u +%Y%m%d%H -d "${data:0:8} ${data:8:2} +${fct} hours")

  fname_grb=/lustre_xc50/ioper/models/SMNA-Oper/SMG/datainout/bam/pos/dataout/TQ0299L064/${data}/GPOSCPT${data}${datafct}P.${type}.TQ0299L064.grb
      
  for varanl in "${VarsAnl[@]}"
  do

    if [[ "${varanl}" == "Zonal wind (u)" ]] ||
       [[ "${varanl}" == "Meridional wind (v)" ]] || 
       [[ "${varanl}" == "Omega" ]] ||
       [[ "${varanl}" == "Stream function" ]] || 
       [[ "${varanl}" == "Velocity potential" ]] || 
       [[ "${varanl}" == "Geopotential height" ]] || 
       [[ "${varanl}" == "Absolute temperature" ]] || 
       [[ "${varanl}" == "Specific humidity" ]] ||
       [[ "${varanl}" == "Vertical dist total cloud cover" ]] ||
       [[ "${varanl}" == "Time mean surface relative humidity" ]]
    then 
      #Levels=(1000 985 975 960 950 930 925 915 900 880 875 850 825 800 775 750 725 700 675 650 600 500 300 250 200 150 100 70 50 30 20 10 3)
      Levels=(1000 925 850 500 250 20 10 3)
    else 
      Levels=(1000)
    fi
  
    for level in ${Levels[@]}
    do

      if [[ "${varanl}" == "Zonal wind (u)" ]] ||
         [[ "${varanl}" == "Meridional wind (v)" ]] || 
         [[ "${varanl}" == "Omega" ]] ||
         [[ "${varanl}" == "Stream function" ]] || 
         [[ "${varanl}" == "Velocity potential" ]] || 
         [[ "${varanl}" == "Geopotential height" ]] || 
         [[ "${varanl}" == "Absolute temperature" ]] || 
         [[ "${varanl}" == "Specific humidity" ]] ||
         [[ "${varanl}" == "Vertical dist total cloud cover" ]] ||
         [[ "${varanl}" == "Time mean surface relative humidity" ]]
      then 
        select_cmd="grb = grib.select(name='${varanl}',level=${level})[0]"
      else 
        select_cmd="grb = grib.select(name='${varanl}')[0]"
      fi

      echo "SMNA ${data} ${varanl} ${level} hPa ${fct} h"

      varanl_fmt=$(echo ${varanl} | sed "s, ,_,g" | sed "s,(,,g" | sed "s,),,g")

      fname_fig="${bpath}/imgs/SMNA/${data}/${varanl_fmt}/${level}/${fct}.png"

      if [ ! -e ${fname_fig} ]
      then

mkdir -p ${bpath}/imgs/SMNA/${data}/${varanl_fmt}/${level}

cat << EOF > ${bpath}/plot_map-smna-${data}-${datafct}.py

#! /usr/bin/env python3

import os
import pygrib
import matplotlib.pyplot as plt
import cartopy, cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
import numpy as np
from cartopy.util import add_cyclic_point

grib = pygrib.open('${fname_grb}')

nx = 900
lon0 = 0.0
dx = 0.4
lons = lon0 + dx * np.arange(nx)

lats = np.array([
-89.69415, -89.29794, -88.89940, -88.50032, -88.10105, -87.70167, -87.30225, -86.90279,
-86.50331, -86.10381, -85.70430, -85.30479, -84.90526, -84.50574, -84.10621, -83.70667,
-83.30714, -82.90760, -82.50806, -82.10851, -81.70897, -81.30943, -80.90988, -80.51033,
-80.11079, -79.71124, -79.31169, -78.91214, -78.51259, -78.11304, -77.71349, -77.31394,
-76.91439, -76.51484, -76.11528, -75.71573, -75.31618, -74.91663, -74.51708, -74.11752,
-73.71797, -73.31842, -72.91886, -72.51931, -72.11976, -71.72020, -71.32065, -70.92110,
-70.52154, -70.12199, -69.72244, -69.32288, -68.92333, -68.52377, -68.12422, -67.72466,
-67.32511, -66.92556, -66.52600, -66.12645, -65.72689, -65.32734, -64.92778, -64.52823,
-64.12867, -63.72912, -63.32956, -62.93001, -62.53045, -62.13090, -61.73134, -61.33179,
-60.93223, -60.53268, -60.13312, -59.73357, -59.33401, -58.93446, -58.53490, -58.13535,
-57.73579, -57.33624, -56.93668, -56.53713, -56.13757, -55.73802, -55.33846, -54.93891,
-54.53935, -54.13980, -53.74024, -53.34069, -52.94113, -52.54157, -52.14202, -51.74246,
-51.34291, -50.94335, -50.54380, -50.14424, -49.74469, -49.34513, -48.94558, -48.54602,
-48.14646, -47.74691, -47.34735, -46.94780, -46.54824, -46.14869, -45.74913, -45.34958,
-44.95002, -44.55046, -44.15091, -43.75135, -43.35180, -42.95224, -42.55269, -42.15313,
-41.75358, -41.35402, -40.95446, -40.55491, -40.15535, -39.75580, -39.35624, -38.95669,
-38.55713, -38.15757, -37.75802, -37.35846, -36.95891, -36.55935, -36.15980, -35.76024,
-35.36069, -34.96113, -34.56157, -34.16202, -33.76246, -33.36291, -32.96335, -32.56380,
-32.16424, -31.76468, -31.36513, -30.96557, -30.56602, -30.16646, -29.76691, -29.36735,
-28.96779, -28.56824, -28.16868, -27.76913, -27.36957, -26.97002, -26.57046, -26.17090,
-25.77135, -25.37179, -24.97224, -24.57268, -24.17313, -23.77357, -23.37401, -22.97446,
-22.57490, -22.17535, -21.77579, -21.37623, -20.97668, -20.57712, -20.17757, -19.77801,
-19.37846, -18.97890, -18.57934, -18.17979, -17.78023, -17.38068, -16.98112, -16.58157,
-16.18201, -15.78245, -15.38290, -14.98334, -14.58379, -14.18423, -13.78468, -13.38512,
-12.98556, -12.58601, -12.18645, -11.78690, -11.38734, -10.98778, -10.58823, -10.18867,
-9.78912,  -9.38956, -8.99001,  -8.59045,  -8.19089,  -7.79134, -7.39178,  -6.99223,
-6.59267,  -6.19311,  -5.79356,  -5.39400,  -4.99445,  -4.59489,  -4.19534,  -3.79578,
-3.39622,  -2.99667,  -2.59711,  -2.19756,  -1.79800,  -1.39845,  -0.99889,  -0.59933,
-0.19978,   0.19978,   0.59933,   0.99889,   1.39845,   1.79800,   2.19756,   2.59711,
2.99667,   3.39622,  3.79578,   4.19534,   4.59489,   4.99445,   5.39400,   5.79356,
6.19311,   6.59267,   6.99223,   7.39178,   7.79134,   8.19089,   8.59045,   8.99001,
9.38956,   9.78912, 10.18867,  10.58823,  10.98778,  11.38734,  11.78690, 12.18645,
12.58601,  12.98556,  13.38512,  13.78468,  14.18423,  14.58379,  14.98334,  15.38290,
15.78245, 16.18201,  16.58157, 16.98112,  17.38068,  17.78023,  18.17979,  18.57934,
18.97890,  19.37846,  19.77801,  20.17757,  20.57712,  20.97668,  21.37623,  21.77579,
22.17535,  22.57490,  22.97446,  23.37401,  23.77357,  24.17313,  24.57268,  24.97224,
25.37179,  25.77135,  26.17090,  26.57046,  26.97002,  27.36957,  27.76913,  28.16868,
28.56824,  28.96779,  29.36735,  29.76691,  30.16646,  30.56602,  30.96557,  31.36513,
31.76468,  32.16424,  32.56380,  32.96335,  33.36291,  33.76246,  34.16202,  34.56157,
34.96113,  35.36069,  35.76024,  36.15980,  36.55935,  36.95891,  37.35846,  37.75802,
38.15757,  38.55713,  38.95669,  39.35624,  39.75580,  40.15535,  40.55491,  40.95446,
41.35402,  41.75358,  42.15313,  42.55269,  42.95224,  43.35180,  43.75135,  44.15091,
44.55046,  44.95002,  45.34958,  45.74913,  46.14869,  46.54824,  46.94780,  47.34735,
47.74691,  48.14646,  48.54602,  48.94558,  49.34513,  49.74469,  50.14424,  50.54380,
50.94335,  51.34291,  51.74246,  52.14202,  52.54157,  52.94113,  53.34069,  53.74024,
54.13980,  54.53935,  54.93891,  55.33846,  55.73802,  56.13757,  56.53713,  56.93668,
57.33624,  57.73579,  58.13535,  58.53490,  58.93446,  59.33401,  59.73357,  60.13312,
60.53268,  60.93223,  61.33179,  61.73134,  62.13090,  62.53045,  62.93001,  63.32956,
63.72912,  64.12867,  64.52823,  64.92778,  65.32734,  65.72689,  66.12645,  66.52600,
66.92556,  67.32511,  67.72466,  68.12422,  68.52377,  68.92333,  69.32288,  69.72244,
70.12199,  70.52154,  70.92110,  71.32065,  71.72020,  72.11976,  72.51931,  72.91886,
73.31842,  73.71797,  74.11752,  74.51708,  74.91663,  75.31618,  75.71573,  76.11528,
76.51484,  76.91439,  77.31394,  77.71349,  78.11304,  78.51259,  78.91214,  79.31169,
79.71124,  80.11079,  80.51033,  80.90988,  81.30943,  81.70897,  82.10851,  82.50806,
82.90760,  83.30714,  83.70667,  84.10621,  84.50574,  84.90526,  85.30479,  85.70430,
86.10381,  86.50331,  86.90279,  87.30225,  87.70167,  88.10105,  88.50032,  88.89940,
89.29794,  89.69415])

#grb = grib.select(name='${varanl}',level=${level})[0]
${select_cmd}

init  = str(grb.analDate)      # Init date / time
run   = str(grb.hour).zfill(2) # Run
valid = str(grb.validDate)     # Valid date / time

tmtmp = grb.values 

tmtmp_plot = np.flipud(tmtmp)

# remove linha branca vertical no meio da figura
data, lons = add_cyclic_point(tmtmp_plot, coord=lons)

plt.figure(figsize=(30,10))
ax = plt.axes(projection=ccrs.PlateCarree())

lon2d, lat2d = np.meshgrid(lons, lats)

#data_min = 600
#data_max = 1000
#interval = 50
#levels = np.arange(data_min, data_max+interval, interval)

vmin = np.nanmin(data)
vmax = np.nanmax(data)
interval = 50
levels = np.linspace(vmin, vmax, interval)

img1 = ax.contourf(lon2d, lat2d, data, cmap='jet', levels=levels, extend='both')

#img2 = ax.contour(lon2d, lat2d, data, colors='white', linewidths=0.3, levels=levels)
#ax.clabel(img2, inline=1, inline_spacing=15, fontsize=10, fmt='%1.0f', colors='white')

shapefile = list(shpreader.Reader('/lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img/br_unidades_da_federacao/BR_UF_2019.shp').geometries())
ax.add_geometries(shapefile, ccrs.PlateCarree(), edgecolor='black', facecolor='none', linewidth=0.3)

ax.coastlines(resolution='10m', color='black', linewidth=0.8)
ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=0.5)
gl = ax.gridlines(crs=ccrs.PlateCarree(), color='gray', alpha=1.0, linestyle='--', linewidth=0.25,
xlocs=np.arange(-180, 180, 5), ylocs=np.arange(-90, 90, 5), draw_labels=True)
gl.top_labels = False
gl.right_labels = False

plt.colorbar(img1, label='${varanl}', orientation='vertical', pad=0.01, fraction=0.05)

plt.title('SMNA: ${varanl} @ ${level} hPa, FCT ${fct} h', fontweight='bold', fontsize=12, loc='left')
plt.title('Valid: ' + valid, fontsize=12, loc='right')

fig_dir = os.path.join('${bpath}/imgs', 'SMNA', str(${data}), str(grb.name).replace(' ','_').replace('(','').replace(')',''), str(${level}))

os.makedirs(fig_dir, exist_ok=True)

#fname_fig = 'SMNA_' + str(grb.name).replace(' ','_').replace('(','').replace(')','') + '_level_' + str(${level}) + 'hPa_valid_for_' + str(valid).replace(' ','_').replace('-','_').replace(':00:00','Z') + '_fct_' + str(${fct}) + 'h'
#fname_fig = 'SMNA/' + str(grb.name).replace(' ','_').replace('(','').replace(')','') + '/' + str(${level}) + '/' + str(${fct})
#fname_fig = os.path.join('${bpath}/imgs', 'SMNA', str(grb.name).replace(' ','_').replace('(','').replace(')',''), str(${level}), str(${fct}))
fname_fig = os.path.join(fig_dir, str(${fct}))

#plt.savefig(fname_fig + '.png', dpi=150, bbox_inches='tight')
#plt.savefig(fname_fig + '.png', bbox_inches='tight')
plt.savefig(
    fname_fig + '.jpg',
    bbox_inches='tight',
    dpi=100,
    format='jpg',
    pil_kwargs={
        'quality': 70,
        'optimize': True,
        'progressive': True,
        'subsampling': 2
    }
)

grib.close()

EOF

      /lustre_xc50/carlos_bastarz/SMNAMonitoringApp/cron_scripts/anls_img/SMNAMonitoringApp/bin/python3 ${bpath}/plot_map-smna-${data}-${datafct}.py

      fi

    done

  done

done

echo "Rodou em: $(date) | UTC: $(date -u)" >> ${bpath}/cron-plot_map-smna-${data}-${datafct}.log

exit 0