% Jupyter Notebook code cell colors
\definecolor{nbsphinxin}{HTML}{307FC1}
\definecolor{nbsphinxout}{HTML}{BF5B3D}
\definecolor{nbsphinx-code-bg}{HTML}{F5F5F5}
\definecolor{nbsphinx-code-border}{HTML}{E0E0E0}
\definecolor{nbsphinx-stderr}{HTML}{FFDDDD}
% ANSI colors for output streams and traceback highlighting
\definecolor{ansi-black}{HTML}{3E424D}
\definecolor{ansi-black-intense}{HTML}{282C36}
\definecolor{ansi-red}{HTML}{E75C58}
\definecolor{ansi-red-intense}{HTML}{B22B31}
\definecolor{ansi-green}{HTML}{00A250}
\definecolor{ansi-green-intense}{HTML}{007427}
\definecolor{ansi-yellow}{HTML}{DDB62B}
\definecolor{ansi-yellow-intense}{HTML}{B27D12}
\definecolor{ansi-blue}{HTML}{208FFB}
\definecolor{ansi-blue-intense}{HTML}{0065CA}
\definecolor{ansi-magenta}{HTML}{D160C4}
\definecolor{ansi-magenta-intense}{HTML}{A03196}
\definecolor{ansi-cyan}{HTML}{60C6C8}
\definecolor{ansi-cyan-intense}{HTML}{258F8F}
\definecolor{ansi-white}{HTML}{C5C1B4}
\definecolor{ansi-white-intense}{HTML}{A1A6B2}
\definecolor{ansi-default-inverse-fg}{HTML}{FFFFFF}
\definecolor{ansi-default-inverse-bg}{HTML}{000000}

% Define an environment for non-plain-text code cell outputs (e.g. images)
\makeatletter
\newenvironment{nbsphinxfancyoutput}{%
    % Avoid fatal error with framed.sty if graphics too long to fit on one page
    \let\sphinxincludegraphics\nbsphinxincludegraphics
    \nbsphinx@image@maxheight\textheight
    \advance\nbsphinx@image@maxheight -2\fboxsep   % default \fboxsep 3pt
    \advance\nbsphinx@image@maxheight -2\fboxrule  % default \fboxrule 0.4pt
    \advance\nbsphinx@image@maxheight -\baselineskip
\def\nbsphinxfcolorbox{\spx@fcolorbox{nbsphinx-code-border}{white}}%
\def\FrameCommand{\nbsphinxfcolorbox\nbsphinxfancyaddprompt\@empty}%
\def\FirstFrameCommand{\nbsphinxfcolorbox\nbsphinxfancyaddprompt\sphinxVerbatim@Continues}%
\def\MidFrameCommand{\nbsphinxfcolorbox\sphinxVerbatim@Continued\sphinxVerbatim@Continues}%
\def\LastFrameCommand{\nbsphinxfcolorbox\sphinxVerbatim@Continued\@empty}%
\MakeFramed{\advance\hsize-\width\@totalleftmargin\z@\linewidth\hsize\@setminipage}%
\lineskip=1ex\lineskiplimit=1ex\raggedright%
}{\par\unskip\@minipagefalse\endMakeFramed}
\makeatother
\newbox\nbsphinxpromptbox
\def\nbsphinxfancyaddprompt{\ifvoid\nbsphinxpromptbox\else
    \kern\fboxrule\kern\fboxsep
    \copy\nbsphinxpromptbox
    \kern-\ht\nbsphinxpromptbox\kern-\dp\nbsphinxpromptbox
    \kern-\fboxsep\kern-\fboxrule\nointerlineskip
    \fi}
\newlength\nbsphinxcodecellspacing
\setlength{\nbsphinxcodecellspacing}{0pt}

% Define support macros for attaching opening and closing lines to notebooks
\newsavebox\nbsphinxbox
\newcommand{\nbsphinxstartnotebook}[1]{%
    \par
    % measure needed space
    \setbox\nbsphinxbox\vtop{{#1\par}}
    % reserve some space at bottom of page, else start new page
    \needspace{\dimexpr2.5\baselineskip+\ht\nbsphinxbox+\dp\nbsphinxbox}
    % mimic vertical spacing from \section command
      \addpenalty\@secpenalty
      \@tempskipa 3.5ex \@plus 1ex \@minus .2ex\relax
      \addvspace\@tempskipa
      {\Large\@tempskipa\baselineskip
             \advance\@tempskipa-\prevdepth
             \advance\@tempskipa-\ht\nbsphinxbox
             \ifdim\@tempskipa>\z@
               \vskip \@tempskipa
             \fi}
    \unvbox\nbsphinxbox
    % if notebook starts with a \section, prevent it from adding extra space
    \@nobreaktrue\everypar{\@nobreakfalse\everypar{}}%
    % compensate the parskip which will get inserted by next paragraph
    \nobreak\vskip-\parskip
    % do not break here
    \nobreak
}% end of \nbsphinxstartnotebook

\newcommand{\nbsphinxstopnotebook}[1]{%
    \par
    \setbox\nbsphinxbox\vbox{{#1\par}}
    \nobreak % it updates page totals
    \dimen@\pagegoal
    \advance\dimen@-\pagetotal \advance\dimen@-\pagedepth
    \advance\dimen@-\ht\nbsphinxbox \advance\dimen@-\dp\nbsphinxbox
    \ifdim\dimen@<\z@
      % little space left
      \unvbox\nbsphinxbox
      \kern-.8\baselineskip
      \nobreak\vskip\z@\@plus1fil
      \penalty100
      \vskip\z@\@plus-1fil
      \kern.8\baselineskip
    \else
    \fi
}% end of \nbsphinxstopnotebook

% Ensure height of an included graphics fits in nbsphinxfancyoutput frame
\newdimen\nbsphinx@image@maxheight % set in nbsphinxfancyoutput environment
\newcommand*{\nbsphinxincludegraphics}[2][]{%
    \gdef\spx@includegraphics@options{#1}%
    \setbox\spx@image@box\hbox{\includegraphics[#1,draft]{#2}}%
    \in@false
    \ifdim \wd\spx@image@box>\linewidth
      \g@addto@macro\spx@includegraphics@options{,width=\linewidth}%
      \in@true
    \fi
    % no rotation, no need to worry about depth
    \ifdim \ht\spx@image@box>\nbsphinx@image@maxheight
      \g@addto@macro\spx@includegraphics@options{,height=\nbsphinx@image@maxheight}%
      \in@true
    \fi
    \ifin@
      \g@addto@macro\spx@includegraphics@options{,keepaspectratio}%
    \fi
    \setbox\spx@image@box\box\voidb@x % clear memory
    \expandafter\includegraphics\expandafter[\spx@includegraphics@options]{#2}%
}% end of "\MakeFrame"-safe variant of \sphinxincludegraphics

\renewcommand*\sphinx@verbatim@nolig@list{\do\'\do\`}
\begingroup
\catcode`'=\active
\let\nbsphinx@noligs\@noligs
\g@addto@macro\nbsphinx@noligs{\let'\PYGZsq}
\endgroup
\renewcommand*\sphinxbreaksbeforeactivelist{\do\<\do\"\do\'}
\renewcommand*\sphinxbreaksafteractivelist{\do\.\do\,\do\:\do\;\do\?\do\!\do\/\do\>\do\-}
\fvset{codes*=\sphinxbreaksattexescapedchars\do\^\^\let\@noligs\nbsphinx@noligs}

A few executable scripts to perform low-level unit tasks bounded to croco

A few executable scripts to perform low-level unit tasks bounded to hycom3d

write_abfiles() (sloop.models.croco.atmfrc.AtmFrc static method)

write_ncfiles() (sloop.models.croco.atmfrc.AtmFrc static method)

get_constant_patterns() (in module sloop.models.hycom3d.__init__)

rename_vars() (sloop.models.hycom3d.atmfrc.AtmFrc static method)

write_abfiles() (sloop.models.hycom3d.atmfrc.AtmFrc static method)

write_ncfiles() (sloop.models.hycom3d.atmfrc.AtmFrc static method)

concat_field_ncfiles() (in module sloop.models.hycom3d.postprod)

get_cmdargs_nc_time_splitting() (in module sloop.models.hycom3d.postprod)

get_interpolation_clargs() (in module sloop.models.hycom3d.postprod)

get_tempconversion_clargs() (in module sloop.models.hycom3d.postprod)

postproc_time_splitting() (in module sloop.models.hycom3d.postprod)

preproc_time_splitting() (in module sloop.models.hycom3d.postprod)

platforms_nc2dataset_preprocess() (sloop.models.hycom3d.rivers.Rivers method)

platforms_nc2datasets() (sloop.models.hycom3d.rivers.Rivers method)

platforms_selection() (sloop.models.hycom3d.rivers.Rivers method)

usage: sloop-conda-install [-h] [--label LABEL] [-a] [-l] [-x EXPORT] [-n]
                           [-p PYTHON_VERSION] [-u USER] [--public]
                           [--profile PROFILE] [-i INSTALL_DIR]
                           [--exec-dir EXEC_DIR]
                           [--miniconda-installer MINICONDA_INSTALLER]
                           [-c CHANNEL] [--override-channels]
                           [--strict-channel-priority] [--no-pip]
                           [--env-dir ENV_DIR] [--env-file ENV_FILE]
                           [--no-env] [--mod-dir MOD_DIR] [-m MOD_NAME]
                           [--mod-as-template] [--no-mod] [-t TESTING_IMPORTS]
                           [--no-color]

Utility for installing conda and alot of useful packages for SLOOP users, locally or on known supercomputers

optional arguments:
  -h, --help            show this help message and exit
  --label LABEL         label used in install directory and modulefile names,
                        which defaults to "perso" or to the name of a meta-
                        package when specified alone (default: sloop)
  -a, --avail           show available meta-packages (default: False)
  -l, --list            only list packages to install (default: False)
  -x EXPORT, --export EXPORT
                        only installation specs as environment file (default:
                        None)
  -n, --dry-run         only show commands to be executed (default: False)
  -p PYTHON_VERSION, --python-version PYTHON_VERSION
                        python version like 3 or 3.6 (default: 3)
  -u USER, --user USER  user name (default: None)
  --public              install everything as a public installation with most
                        of the paths prefixed with "Public/" (default: False)
  --profile PROFILE     installation profile, whose default is detected; list
                        of profiles: path, datarmor, home, default,
                        meteofrance, localhome, where path is a file path with
                        a "/" inside (default: home)
  -i INSTALL_DIR, --install-dir INSTALL_DIR
                        installation directory, which depends on the
                        installation profile (default: None)
  --exec-dir EXEC_DIR   directory where to download the installer and execute
                        it, which defaults to home dir (default: ~)
  --miniconda-installer MINICONDA_INSTALLER
                        miniconda url or local file (default: https://repo.ana
                        conda.com/miniconda/Miniconda3-py37_4.8.2-Linux-x86_64
                        .sh)
  -c CHANNEL, --channel CHANNEL
                        choose channels manually, instead of default ones:
                        conda-forge (default: None)
  --override-channels   do not override user config conda channels (default:
                        False)
  --strict-channel-priority
                        do not pass option --strict-channel-priority (default:
                        False)
  --no-pip              do not install pip packages (default: False)
  --env-dir ENV_DIR     where to save the installed conda env (default: None)
  --env-file ENV_FILE   basename of the environment yml file (default:
                        {flat_install_dir}-{env_type}.yml)
  --no-env              do not save the installed environment (default: False)
  --mod-dir MOD_DIR     source directory for modulefiles (default: None)
  -m MOD_NAME, --mod-name MOD_NAME
                        modulefile name to create (default: conda3/sloop)
  --mod-as-template     save the module file as a template with the
                        {install_dir} pattern (default: False)
  --no-mod              do not create a modulefile (default: False)
  -t TESTING_IMPORTS, --testing-imports TESTING_IMPORTS
                        imports to check the installation as a comma separated
                        list of modules (default: None)
  --no-color            do not colorize outputs (default: False)

examples of usage::

    sloop-conda-install -n  # dry run means inform but do nothing, just to try
    sloop-conda-install

Installation profile: home
Chdir: /Users/mcaillau
Install in new directory: /Users/mcaillau/Soft/conda3-sloop
Download miniconda: https://repo.anaconda.com/miniconda/Miniconda3-py37_4.8.2-Linux-x86_64.sh
Install conda with: bash Miniconda3-py37_4.8.2-Linux-x86_64.sh -b -p /Users/mcaillau/Soft/conda3-sloop
Install conda packages with: /Users/mcaillau/Soft/conda3-sloop/bin/conda install -q -y --override-channels --strict-channel-priority -c conda-forge appdirs cfgrib colorlog dask esmpy=8.2.0 gitpython htop mpi4py mpich nose python=3.7.6 pyyaml six xarray=0.16.0 xesmf=0.6.0 xoa=0.4.0
Export env with: /Users/mcaillau/Soft/conda3-sloop/bin/conda env export -n sloop --file /Users/mcaillau/Soft/conda3-sloop/environment.yml
Create modulefile: /Users/mcaillau/Modulefiles/conda3/sloop
Installation successfull
############################################################
/Users/mcaillau/Soft/conda3-sloop/bin/python
/Users/mcaillau/Soft/conda3-sloop/bin/conda
module use /Users/mcaillau/Modulefiles
module load conda3/sloop

Traceback (most recent call last):
  File "/Users/mcaillau/Documents/ETUDES/PROJET_MED/SRC_CODE/sloop/bin/sloop", line 14, in <module>
    from sloop.cli.main import main
  File "/Users/mcaillau/Documents/ETUDES/PROJET_MED/SRC_CODE/sloop/sloop/cli/main.py", line 14, in <module>
    from ..io import CacheManager, SloopArgumentHelpFormatter
  File "/Users/mcaillau/Documents/ETUDES/PROJET_MED/SRC_CODE/sloop/sloop/io.py", line 16, in <module>
    from appdirs import AppDirs
ModuleNotFoundError: No module named 'appdirs'

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
CROCO utilities

"""
import logging
import os

from ...__init__ import SLOOP_SRC_DIR, SloopError

class SloopCrocoError(SloopError):
    pass

LOGGER = logging.getLogger(__name__)

_THISDIR = os.path.dirname(__file__)

#: CROCO vconf
CROCO_VCONF = ["manga", "med", "indien", "gijon"]

#: CROCO platform root dir
CROCO_DIR = os.path.join(SLOOP_SRC_DIR, 'CROCO')

#: CROCO GENERAL CONFI INI FILE
CROCO_CONF_INI_FILE = os.path.join(_THISDIR, "croco_config.ini")

#: CROCO initial and open boundary root dir
CROCO_IBC_DIR = os.path.join(CROCO_DIR, 'EXTRACT')

#: CROCO EXTRACT EXECUTABLE
CROCO_IBC_EXE = os.path.join(CROCO_IBC_DIR, 'extract_exe')

#: CROCO INPUTS DIR
CROCO_INPUTS_DIR='./FORCING'

#: CROCO_RESTART_NAME
CROCO_RST_NAME='croco_rst.nc'

#: CROCO EXTRACT TEMPLATE NAMELIST
CROCO_IBC_NAMELIST_TPL=os.path.join(CROCO_IBC_DIR,'IN/namelist.tpl')

#: CROCO EXTRACT INI FILE  
CROCO_IBC_INPUT_FILE=os.path.join(CROCO_INPUTS_DIR,'obc_collect/INPUT_EXTRACT.nc')

#: CROCO EXTRACT INI FILE
CROCO_EXTRACT_INI_FILE = os.path.join(_THISDIR, "extract.ini")

#: CROCO grid  file
CROCO_GRID_FILE = os.path.join(CROCO_INPUTS_DIR,'croco_grd.nc')

#: CROCO croco.in parameters file
CROCO_RUN_INPUT_FILE = "TPL/croco.in.tpl"

#CROCO FORCING FILE
CROCO_FRC_TIDE=os.path.join(CROCO_INPUTS_DIR,"croco_frc.nc")

#CROCO ATM FORCING FILE
CROCO_FRC_ATM=os.path.join(CROCO_INPUTS_DIR,"croco_atm.nc")

#CROCO_RUNOFF_FILE
CROCO_FRC_RIVER=os.path.join(CROCO_INPUTS_DIR,"croco_runoff.nc")

#CROCO CPPDEF FILE
CROCO_PREPROC_FILE="TPL/cppdefs.h.tpl"

#CROCO PARAM FILE
CROCO_PARAM_FILE="TPL/param.h.tpl"

#CROCO COMPILE FILE
CROCO_COMPILE_FILE="TPL/jobcomp.tpl"

#: CROCO 3D model directory
CROCO_MODEL_DIR = os.path.join(CROCO_DIR, 'Modele3d')

#: CROCO 3D model executable
CROCO_MODEL_EXE = os.path.join(
    CROCO_MODEL_DIR, 'croco')

#: CROCO rivers directory
CROCO_RIVERS_DIR = os.path.join(CROCO_DIR, 'Rivers')

#: CROCO rivers flow climatology
CROCO_RIVERS_FLOW_CLIM = os.path.join(CROCO_RIVERS_DIR, "hycom_rivers.d")

#: CROCO rivers hydro climatology
CROCO_RIVERS_HYDRO_CLIM = os.path.join(
    CROCO_RIVERS_DIR, "fleuves_climato_hydro")

CROCO_SIGMA_TO_STMT_FNS=""
CROCO_CONTEXT_TO_STMT_FNS=""
CROCO_IBC_GEN_INICON=""
CROCO_IBC_GEN_INICON_VCONF=""
CROCO_MODEL_DEFSTRETCH_FILES=""
CROCO_GRID_AFILE=""
CROCO_BLKDAT_INPUT_FILE=""

#: CROCO postproduction root dir
CROCO_POSTPROD_DIR = os.path.join(CROCO_DIR, 'PostProduction')

#: CROCO postproduction time filtering executable
CROCO_POSTPROD_FILTER = os.path.join(
    CROCO_POSTPROD_DIR, 'cdfTraitementTemporel3D')

#: CROCO postproduction file conversion executable
CROCO_POSTPROD_INTERPOLATION = os.path.join(
    CROCO_POSTPROD_DIR, 'cdfHycomTodataSHOMouSOAP')

#: CROCO postproduction temperature conversion executable
CROCO_POSTPROD_TEMPCONVERSION = os.path.join(
    CROCO_POSTPROD_DIR, 'temppot2tempis_postProd')

#: CROCO postproduction traductions noms longs
CROCO_POSTPROD_LONGNAMES = "PARAMETERS{rank}./traductions_noms_longs_{postprod}"

#: CROCO postproduction traductions noms standards
CROCO_POSTPROD_STDNAMES = "PARAMETERS{rank}./traductions_noms_standards_2"

#: CROCO postproduction traductions units
CROCO_POSTPROD_UNITS = "PARAMETERS{rank}./traductions_unites"

#: CROCO savefield.input parameters file
CROCO_SAVEFIELD_INPUT_FILE = "PARAMETERS{rank}./savefield.input"

#: CROCO blkdat_cmo.input parameters file
CROCO_SAVEFIELD_INPUT_FILES = {
            "spinup":["PARAMETERS{rank}./savefield_C.input",
                      "PARAMETERS{rank}./savefield_A.input"],
            "forecast":["PARAMETERS{rank}./savefield_C.input",
                        "PARAMETERS{rank}./savefield_P.input"],
            "spnudge_free":["PARAMETERS{rank}./savefield_Assim_SN.input"],
            "spnudge_relax":["PARAMETERS{rank}./savefield_C.input",
                             "PARAMETERS{rank}./savefield_A.input"],
            "free":["PARAMETERS{rank}./savefield_C.input",
                    "PARAMETERS{rank}./savefield_A.input"],
            }

#: CROCO MPI tiling file
CROCO_PATCH_INPUT_FILE = "SPLIT{rank}./patch.input"

#: CROCO run.input.tpl parameters file
CROCO_RUN_INPUT_TPL_FILE = "PARAMETERS{rank}./run.input.tpl"

#: CROCO post-production configuration files
CROCO_POSTPROD_CFG_FILE = "PARAMETERS{rank}./postprod.cfg"

#: CROCO post-production configuration specs file
CROCO_POSTPROD_INI_FILE = os.path.join(_THISDIR, "postprod.ini")

#: CROCO post-production configuration output files
CROCO_POSTPROD_CONFIG_FILE = "PARAMETERS{rank}./postprod_{postprod}.cfg"

#: CROCO check promises
CROCO_PROMISES = "PARAMETERS{rank}./promises.cfg"

#: CROCO rivers configuration specs file
CROCO_RIVERS_INI = os.path.join(_THISDIR, "rivers.ini")

#: CROCO rivers configuration file
CROCO_RIVERS_CFG = 'FORCING{rank}./nest/rivers.cfg'

#: CROCO spectral nudging parameters json file
CROCO_SPNUDGE_JSON_FILE = "PARAMETERS{rank}./spnudging_parameters.json"

#: From standard_name to name
CROCO_STANDARD_NAME_TO_NAME = {
    'sea_water_potential_temperature': 'temp',
    'sea_water_salinity': 'saln',
    'sea_surface_height_above_geoid': 'ssh',
    }

#: Required environment variables for each context
REQUIRED_ENV_VARIABLES = {
    # Pre- and post-processing tasks
    'prepost': [
        'FC', 'LDF', 'CC', 'FFLAGS', 'LDFLAGS', 'CPPFLAGS', 'MAKE',
        'MODEL_FFLAGS', 'MODEL_CFLAGS', 'MODEL_LDFLAGS'
        ],
    # Model
    'modelcompile': [
    
        'HOME_MODELE_SRC', 'MPIDIR','MPI_COMP',
        'FFLAGS', 'CPP', 'CPPFLAGS', 'LDFLAGS', 'CCFLAGS',
        'FC', 'CC', 'MAKE'
        ]
    ,
    'modelrun': [
        'HOME_MODELE_SRC', 'MPIDIR',
        ],
    }

REQUIRED_JOB_OPTIONS = {
        'submitcmd', 'queue', 'time'
        ],
        ],
    'recdata': [
        'submitcmd', 'queue', 'refill', 'time'
        ],
    }

#: Patterns of geographic constant files
CROCO_GEO_CONSTANT_PATTERNS = {
    "river": CROCO_FRC_RIVER,
    "tide": CROCO_FRC_TIDE, 
    "atm": CROCO_FRC_ATM,
    "grid":CROCO_GRID_FILE,
}

# Pattern of TPL files
CROCO_TPL_FILES= {
    "param": CROCO_PARAM_FILE,
    "preproc":CROCO_PREPROC_FILE,
    "input":CROCO_RUN_INPUT_FILE,
    "compile":CROCO_COMPILE_FILE,
}
"""""
    "run": [
        CROCO_RUN_INPUT_FILE,
        ],
    "postprod": [
        CROCO_POSTPROD_CFG_FILE,
        CROCO_POSTPROD_CONFIG_FILE.format(
            rank="{rank}", postprod="*")
    ],
    "savefield": (
        [CROCO_SAVEFIELD_INPUT_FILE] +
        list(set(p for pp in CROCO_SAVEFIELD_INPUT_FILES.values()
                 for p in pp))
    ),
    "split": [
        CROCO_PATCH_INPUT_FILE
    ],
    "promises": [
        CROCO_PROMISES
    ],
}
"""""
#: Patterns of non-geographic constant files
CROCO_NOT_GEO_CONSTANT_PATTERNS = {
    "naming": [
        "PARAMETERS{rank}./traductions_*",
    ]
}

#: CROCO specifications for the :mod:`xoa.cf` module use
CROCO_CF_SPECS = {
    "data_vars": {
        "temp": {
            "name": "temp"
        },
        "sal": {
            "name": "salt"
        },
        "ssh": {
            "name": "zeta"
        },
        "sigmatheta": {
            "name": "thdd",
        },
        "ps": {
            "alt_name": ["preatm", "prmsl"],
            "attrs": {
                "standard_name": ["air_pressure_at_mean_sea_level"],
                "long_name": ["Pressure reduced to MSL"],
                "header": "pression atmospherique"
            }
        },
        "tair": {
            "alt_name": ["airtmp", "t2m"],
            "attrs": {
                "units": "degree Celsius",
                "standard_name": [
                    "air_temperature",
                    "air_temperature_at_2m"
                ],
                "long_name": "2 metre temperature",
                "header": "temperature de surface"
            }
        },

    }
}
CROCO_PARAM = {

    "u10m": {
        "name": [
            "u10m",
            "u10",
            "u"
        ],
        "attrs": {
            "units": "m/s",
            "standard_name": "eastward_wind",
            "long_name": "Zonal wind speed at 10m",
            "header": "zonal wind"
        }
    },
    "v10m": {
        "name": [
            "v10",
            "v"
        ],
        "attrs": {
            "standard_name": "northward_wind",
            "long_name": "Meridional wind speed at 10m",
            "header": "Meridional wind"
        }
    },
    "pmer": {
        "name": [
            "patm",
            "pmsl",
            "preatm",
            "prmsl"
        ],
        "attrs": {
            "units": "Pa",
            "standard_name": "air_pressure_at_mean_sea_level",
            "long_name": "Pressure reduced to MSL",
            "header": "pression atmospherique"
        }
    },
    "t2m": {
        "name": [
            "airtmp",
            "t2m"
        ],
        "attrs": {
            "units": "degree Celsius",
            "standard_name": [
                "air_temperature",
                "air_temperature_at_2m"
            ],
            "long_name": "2 metre temperature",
            "header": "temperature de surface",
            "diurnal": "True"
        }
    },
    "flsolaire": {
        "name": [
            "flsolaire",
            "swhf",
            "shwflx",
            "ssr"
        ],
        "attrs": {
            "units": "W m-2 s-1",
                "short_wave_fluxes"
            ],
            "long_name": "surface net solar radiation",
            "header": "flux radiatifs Short Wave",
        }
    },
    "fltherm": {
        "name": [
            "fltherm",
            "lwhf",
            "radflx",
            "str",
       ],
        "attrs": {
                "radiative_fluxes"
            ],
            "long_name": "total radiative fluxes",
            "header": "flux radiatifs IR",
        }
    },
    "prate": {
        "name": [
            "eau",
            "prate",
            "rain"
        ],
        "attrs": {
            "units": "m s-1",
                "precipitations"
            ],
            "long_name": "Total Precipitation",
            "header": "precipitations"
        }
    },
    "hu2m": {
        "name": [
            "hu2m",
            "rh",
            "rh2m"
        ],
        "attrs": {
            "units": "frac",
                "relative_humidity"
            ],
            "long_name": "Relative humidity at 2m",
            "header": "Relative humidity"
        }
    },
    "wspd": {
        "name": [
            "wdspd",
        ],
        "attrs": {
                "wind_speed"
            ],
            "long_name": "10 metre wind speed",
            "header": "module du vent"
        }
    },
}

_CACHE = {}

def get_tpl_files():
    """Get dict of CROCO TPL files"""
    return CROCO_TPL_FILES

def get_constant_patterns(geo):
    """Get dict of CROCO constant file glob patterns for each pack

    A file pattern typically contains the `{rank}` subsitution pattern

    Parameters
    ----------
    geo: bool
        Bounded to a geograhpical area?
    """
    return CROCO_GEO_CONSTANT_PATTERNS if geo else CROCO_NOT_GEO_CONSTANT_PATTERNS

def get_cf_specs():
    """Get a :class:`xoa.cf.CFSpecs` object for Croco"""
    if "cf_specs" not in _CACHE:
        import xoa.cf
        _CACHE["cf_specs"] = xoa.cf.CFSpecs(CROCO_CF_SPECS)
    return _CACHE["cf_specs"]

def format_ds(ds):
    """Format all variables in a dataset or dataarray to be Croco compliant"""
    # Time units
    time = get_cf_specs().coords.get(ds, "time")
    if time is not None:
        ds[time.name].encoding["units"] = "seconds since 1900-01-01"

    # Other encodings
    for vname, var in ds.items():
        for attr in 'dtype', 'scale_factor', 'add_offset':
            if attr in ds[vname].encoding:
                del ds[vname].encoding[attr]
        if ds[vname].dtype.kind in "fd":
            ds[vname].encoding['_FillValue'] = -1e34

    return ds

"""
Atmospheric Forcing preprocessings
"""

import numpy as np
from ...log import LoggingLevelContext
import pandas as pd
import glob

with LoggingLevelContext(logging.WARNING):
    import xarray as xr
 
from xoa.filter import erode_mask

from ...times import convert_to_julian_day
from ...interp import interp_time, Regridder
from ...phys import (
    windstress,
    radiativeflux,
    celsius2kelvin,
    watervapormixingratio,
)
from ...filters import erode_mask_vec
from .__init__ import CROCO_PARAM
from .io import get_ab_file_names

class AtmFrc(object):
    """ Atmospheric Forcing preprocessings """

    def __init__(
        self, insta_files, cumul_files,
    ):
        """
        Parameters
        ----------
        insta_files: list
            List of instantaneous grib files::

                insta_files = [grib_time0, grib_time1, grib_time2, ...]

        cumul_files: list(list)
            List of list of cumulated grib files.
            The outer level is for dates and the inner one for terms::

                cumul_files = [[grib_date0_term0, grib_date0_term1],
                               [grib_date1_term0, grib_date1_term1],
                               [grib_date2_term0, grib_date2_term1]]

            Note that the cumul is reset at each new date.

        """
        self.insta = insta_files
        self.cumul = cumul_files
        LOGGER.debug("insta files: {self.insta}")
        LOGGER.debug("cumul files: {self.cumul}")

    def grib2dataset(self):
        LOGGER.info("Convert grib format to dataset")
        dsc = []
        LOGGER.info("Cumulated data processing")
        for grib_at_terms in self.cumul:
            LOGGER.debug(grib_at_terms)
            ds = []
            for data in grib_at_terms:
                ds.append(
                    xr.open_dataset(
                        data, engine="cfgrib", backend_kwargs={"indexpath": ""}
                    )
                )
            ds = xr.combine_nested(
                ds, concat_dim="valid_time", combine_attrs="override"
            )
            ds = ds.drop(["time", "surface"])
            ds = ds.rename(
                {"longitude": "lon", "latitude": "lat", "valid_time": "time"}
            )
            ds = ds.sortby(ds["time"])
            ds = ds.astype("float64")
            ds = self.decumul(ds)
            dsc.append(ds)
        ds_cumul = xr.combine_nested(
            dsc, concat_dim="time", combine_attrs="override"
        )
        ds_cumul = ds_cumul.sortby(ds_cumul["time"])
        LOGGER.info("Instantaneous data processing")
        dsi = []
        for data in self.insta:
            ds = []
            for level in [0, 2, 10]:
                        data, engine="cfgrib",
                        backend_kwargs={
                            "indexpath": "",
                            "filter_by_keys": {"edition": 2, "level": level},
                        },
                        drop_variables=[
                            "heightAboveGround",
                            "meanSea",
                            "surface",
                        ],
                    )
                )
            dsi.append(xr.merge(ds, combine_attrs="override"))
        ds_insta = xr.combine_nested(
            dsi, concat_dim="valid_time", combine_attrs="override"
        )
        ds_insta = ds_insta.drop(["step", "time"])
        ds_insta = ds_insta.rename(
            {"longitude": "lon", "latitude": "lat", "valid_time": "time"}
        )
        ds_insta = ds_insta.sortby(ds_insta["time"])
        return ds_cumul, ds_insta

    @staticmethod
    def decumul(ds_cumul):
        LOGGER.debug("start decumul")
        if ds_cumul.dims["time"] == 1:
            LOGGER.debug("one time step")
            diff_time = ds_cumul["step"]
            ds_decumul = ds_cumul/(diff_time/np.timedelta64(1, "s"))
        else:
            LOGGER.debug("more than one time step")
            ds_decumul = xr.concat([ds_cumul.isel(time=0), ds_cumul.diff("time")], "time")
            diff_time = ds_decumul.step.values/np.arange(1, ds_decumul.dims["time"]+1, 1)
            ds_decumul = ds_decumul/(diff_time/np.timedelta64(1, "s"))
        ds_decumul["time"] = ds_decumul.time - 0.5 * diff_time
        if "step" in ds_decumul:
            ds_decumul = ds_decumul.drop(["step"]) 
        return ds_decumul

    def interp_time(self, time, nc_out):
        ds_cumul, ds_insta = self.grib2dataset()
        diurnal_vars = []
        for pv in CROCO_PARAM.values():
            if "diurnal" in pv["attrs"]:
                if pv["attrs"]["diurnal"] is "True":
                    diurnal_vars.extend(pv["name"])
        LOGGER.debug(f"Diurnal vars : {diurnal_vars}")
        LOGGER.info("time interpolation of instantaneous data")
        ds_insta = interp_time(ds_insta, time, diurnal_vars=diurnal_vars)
        LOGGER.info("time interpolation of cumulated data")
        ds_cumul = interp_time(ds_cumul, time, diurnal_vars=diurnal_vars)
        ds_atmfrc = xr.merge(
            [ds_cumul, ds_insta], combine_attrs="override", compat='no_conflicts',
        )
        LOGGER.info(f"Store results in {nc_out}")
        ds_atmfrc.to_netcdf(nc_out)
         
    def parameters(nc_in, nc_out):
        ds = xr.open_dataset(nc_in)
        ds = celsius2kelvin(ds)
        if not set(["ewss", "nsss"]).issubset([n for n in ds.variables]):
            ds = windstress(ds, method="Speich")
        ds = radiativeflux(ds)
        ds.drop("str")
        ds = watervapormixingratio(ds)
        ds["tp"] *= 0.001  # [kg.m-2.s-1]-->[m.s-1]
        ds["tp"] = ds["tp"].where(ds["tp"]>=0, 0)
        ds.to_netcdf(nc_out)

    def regridmask(nc_mask, nc_in, nc_out, weightsfile):
        """Add an atmos mask to the dataset"""
        LOGGER.info("Add atmospheric mask to the dataset")
        ds = ds.reindex({'lat': ds.lat[::-1]})
        ds_mask = xr.open_dataset(nc_mask)
        ds = ds.sel(lon=slice(ds_mask.lon.min(), ds_mask.lon.max(),1),
                    lat=slice(ds_mask.lat.min(), ds_mask.lat.max(),1))
        regridder = Regridder(
            ds_mask, ds, regridder=None, filename=weightsfile)
        ds = regridder.regrid(ds_mask)
        ds["mask"] = (ds["mask"] <= 0.5).astype('i')
        ds["mask"].attrs.update(
            long_name="Land-sea mask with 1 on the ocean and 0 on land")
        LOGGER.debug(ds)
         
    def regridvar(nc_in, nc_out, hycom_grid, weightsfile):
        """Regrid to ocean grid with appropriate handling of coasts"""
        kernel = {"lon": 3, "lat": 3}
        ds = erode_mask_vec(ds, kernel, param="wind", until=5)
        scalar_variables = [
            "t2m",
            "r2",
            "prmsl",
            "ssr",
            "tp",
            "si10",
            "vapmix",
            "radflx",
        ]
        ds = ds[scalar_variables+["ewss", "nsss", "mask"]]
        for var in scalar_variables:
            LOGGER.info(f"Erode mask on {var}")
            ds[var] = ds[var].where(ds.mask == 1)
            ds[var] = erode_mask(ds[var], kernel=kernel, until=5)
        ds = ds.transpose('time', 'lat', 'lon')
        LOGGER.info("Build regridder")
            ds, hycom_grid, regridder=None, filename=weightsfile
        )
        LOGGER.info("Redrig the dataset")
        ds = regridder.regrid(ds)
        LOGGER.info(f"Store dataset in NetCDF file {nc_out}")

    def rename_vars(ds):
        LOGGER.info("Rename atmfrc variables as expected by Hycom")
        for hyparam in list(CROCO_PARAM.keys()):
            if len(CROCO_PARAM[hyparam]["name"]) > 1:
                CROCO_PARAM[hyparam]["name"].remove(hyparam)
                if CROCO_PARAM[hyparam]["name"][0] in ds:
                    ds = ds.rename({CROCO_PARAM[hyparam]["name"][0]: hyparam})
            ds[hyparam].attrs = CROCO_PARAM[hyparam]["attrs"]
        for dsi in list(ds.keys()):
            if dsi not in list(CROCO_PARAM.keys()):
                ds = ds.drop(dsi)
        return ds

    def write_abfiles(ds, ab_out="forcing.{var_name}", **kwargs):
        LOGGER.info("Write Hycom [a/b] files")
        ds = ds.fillna(0)
        freq = kwargs.get("freq", None)
        fa_pattern, fb_pattern = get_ab_file_names(ab_out)
        header_pattern = (
            "{var_header} a {freq} heures \n\n"
            "ORIGINE:unknow \n"
            " II   JJ    XLONE    XLONW    XLATS    XLATN \n"
            "{ds.lon.shape[1]}  {ds.lon.shape[0]}    "
            "{lonmax}  {lonmin}   "
            "{latmin}   {latmax} \n"
        )
        line_pattern = (
            "{var_name}: date,range =     {date}     {varmin}     {varmax}\n"
        )
        rec = int(np.ceil(float(ds.lon.size * 4) / 2 ** 14) * 2 ** 14)
        valformat = "{:.3f}"
        lonmin = valformat.format(ds.lon.values.min())
        lonmax = valformat.format(ds.lon.values.max())
        latmin = valformat.format(ds.lat.values.min())
        latmax = valformat.format(ds.lat.values.max())
        for var_name in list(ds.keys()):
            LOGGER.info(f"{var_name}")
            var_header = ds[var_name].attrs["header"]
            header = header_pattern.format(**locals())
            valformat = "{:.8E}"
            with open(fa_pattern.format(**locals()), "wb") as fa:
                with open(fb_pattern.format(**locals()), "w") as fb:
                    LOGGER.debug(header)
                    fb.writelines(header)
                    for itime in range(len(ds.time)):
                        ds[var_name][itime, :, :].values.astype(">f4").tofile(
                            fa
                        )
                        fa.seek(rec * (itime + 1))
                        date = valformat.format(
                            convert_to_julian_day(ds.time.values[itime])
                        )
                        varmin = valformat.format(
                            ds[var_name][itime, :, :].values.min()
                        )
                        varmax = valformat.format(
                            ds[var_name][itime, :, :].values.max()
                        )
                        fb.write(line_pattern.format(**locals()))

    def write_ncfiles(ds, nc_out="forcing.{var_name}.nc", **kwargs):
            LOGGER.info("Create {}".format(
                nc_out.format(**locals())))
            ds[var_name].to_netcdf(nc_out.format(**locals()))

    def outgen(nc_in, **kwargs):
        ds = AtmFrc.rename_vars(ds)
        AtmFrc.write_abfiles(ds, freq=freq)
        AtmFrc.write_ncfiles(ds)

meteo_models={
    "arpege-hr": {
        "data_dir":"/home/ref-oc-intranet-restricted/modeles_mf/arpege_hr/best_estimate",
        "prefix":"METEOFRANCE_ARPEGE-HR_",
        "suff_inst":".nc",
        "suff_flux":"_CUMUL.nc",
        "freq":"1H",
        "varnames":{
            'lwhf':'fltherm',
            'swhf':'flsolaire',
            'prate':'eau',
            'u10':'u10m',
            'v10':'v10m',
            'pressure':'pmer'

        }
    }
}
def decumul_flux(ds_cumul,model):
    """Decumul Atm fluxes"""
    LOGGER.info("Decumul flux files")
    #get the data frequency (cumulated period)
    mod=meteo_models[model]
    step=mod['freq']
    #compute time diff in seconds
    diff_time=pd.Timedelta(step).total_seconds()
    #decumul data
    ds_decumul=ds_cumul/diff_time
    #shift time to match the average period of decumulated data
    ds_decumul["time"] = ds_decumul.time - 0.5 * pd.Timedelta(step) 
    #now interp data on the original time axis
    ds_decumul=ds_decumul.interp(time=ds_cumul.time)
    return ds_decumul

def concat_met_files(ds_inst,ds_cumul):
    """Combine meteo files in one Dataset"""
    LOGGER.info("Concat inst and flux files")
    ds_atmfrc = xr.merge(
            [ds_cumul, ds_inst], combine_attrs="override", compat='no_conflicts')
    return ds_atmfrc

def rename_vars(ds):
    """Rename meteo variable to match Croco model"""
    LOGGER.info("Rename atmfrc variables as expected by Croco")
    for var in list(ds.keys()):
     #test if the variable is already an hyperparam => no rename , just add attrs
     if var in list(CROCO_PARAM.keys()):
       ds[var].attrs = CROCO_PARAM[var]["attrs"]
     else:
       #need to browse every hyperparam's names list  
       for hyparam in list(CROCO_PARAM.keys()):
         if var in CROCO_PARAM[hyparam]["name"]:
           ds = ds.rename({var: hyparam})
           ds[hyparam].attrs = CROCO_PARAM[hyparam]["attrs"]
           break
    #remove unused variables
    for dsi in list(ds.keys()):
       if dsi not in list(CROCO_PARAM.keys()):
          ds = ds.drop(dsi)
    return ds

def format_to_croco(ds,model):
    """Format variables as expected by Croco"""
    LOGGER.info("Format to CROCO")
    #get the model params
    #change sign of long wave flux (positive downard)
    lwhf_name=mod['varnames']['lwhf']
    ds[lwhf_name]=ds[lwhf_name]*-1.0
    #rename vars
    ds=rename_vars(ds)
    #rename longitude,latitude to lon,lat
    ds = ds.rename(
                {"longitude": "lon", "latitude": "lat"}
            )
    #look for the name of rain in CROCO_FORMAT
    meteo_names=list(CROCO_PARAM.keys())
    pnames=["eau","prate","rain"]
    pname=set(pnames).intersection(set(meteo_names))
    pname=''.join(pname)
    #if prate is present put it to 0. else create prate
    if pname in list(ds.keys()):
        ds[pname]=ds[pname]*0.0
    else:
        ds[pname]=ds[list(ds.keys())[0]]*0.0
    return ds

def build_meteo_list(date1,date2,model):
    """Build the list of meteo files from given dates"""
    LOGGER.info("Get list of meteo files")
    #retrieve params of meteo model
    rac=mod["prefix"]
    data_dir=mod['data_dir']
    suff_inst=mod['suff_inst']
    suff_flux=mod['suff_flux']
    #build list of date
    dates_list=pd.date_range(start=date1,end=date2,freq=mod["freq"])
    #build list of names
    inst_list=[f"{rac}{date.strftime('%Y%m%dT%HZ')}{suff_inst}" for date in dates_list]
    flux_list=[f"{rac}{date.strftime('%Y%m%dT%HZ')}{suff_flux}" for date in dates_list]

    #finally build list of files
    inst_list=[glob.glob(f"{data_dir}/**/{file}")[0] for file in inst_list]
    flux_list=[glob.glob(f"{data_dir}/**/{file}")[0] for file in flux_list]

    return inst_list,flux_list

def read_atm_data(inst_list,cumul_list,lon_bnds,lat_bnds):
    """Read list of meteo files and reduce it to the area of interest"""
    from functools import partial
    def _preprocess(x, lon_bnds, lat_bnds):
        return x.sel(longitude=slice(*lon_bnds), latitude=slice(*lat_bnds))
    partial_func = partial(_preprocess, lon_bnds=lon_bnds, lat_bnds=lat_bnds)
    ds_inst=xr.open_mfdataset(inst_list,concat_dim="time",combine="nested",preprocess=partial_func)
    ds_cumul=xr.open_mfdataset(cumul_list,concat_dim="time",combine="nested",preprocess=partial_func)
    ds_inst=ds_inst.squeeze(dim='height')
    ds_cumul=ds_cumul.squeeze(dim='height')
    return ds_inst,ds_cumul

"""
In/output routines

Source for .a/.b handling: http://code.env.duke.edu/projects/mget/wiki/HYCOM

"""

import re
import math
import os
import shutil

    import pandas as pd
import scipy.io as sio

from .__init__ import (
    get_cf_specs,
    SloopCrocoError,
    CROCO_SIGMA_TO_STMT_FNS,
    CROCO_CONTEXT_TO_STMT_FNS,
    CROCO_IBC_GEN_INICON,
    CROCO_IBC_GEN_INICON_VCONF,
    CROCO_MODEL_DEFSTRETCH_FILES
)
from sloop.times import convert_from_julian_day, convert_to_julian_day

BIN2HYCOM_SPECS = {
    "h": {
        "header": "epaisseurs de couches deduites des donnees Mercator",
        "positive": True,
    },
    "temp": {
        "header": "temperatures deduites des donnees Mercator",
        "positive": False,
    },
    "saln": {
        "header": "salinites deduites des donnees Mercator",
    },
    "u": {
        "header": "vitesses zonales deduites des donnees Mercator",
    },
    "v": {
        "header": "vitesses meridiennes deduites des donnees Mercator",
    },
}

def get_ab_file_names(path):
    """Get .a and .b file name from a single path

    path: str
        Root, .a file or .b file

    Return
    ------
    str: afile
    str: bfile

    """
    m = re.match(r"(?P<root>.+)\.[ab]$", path)
    if m is not None:
        root = m.group(1)
    else:
        root = path
        if root.endswith("."):
            root = root[:-1]
    afile = root + ".a"
    bfile = root + ".b"
    return afile, bfile

def read_regional_grid_b(grid_b_file):
    """Read the regional.grid.b file

    grid_b_file: str

    Return
    ------
    dict
    """
    _, grid_b_file = get_ab_file_names(grid_b_file)

    grid_specs = {}
    with open(grid_b_file, "rU") as f:
        grid_specs["idm"] = int(f.readline().split()[0])  # rows
        grid_specs["jdm"] = int(f.readline().split()[0])  # cols
        grid_specs["mapflg"] = int(f.readline().split()[0])
        for line in f:
            if line:
                sline = line[:-1].split()
                key = sline[0][:-1]
                grid_specs[key + "_min"] = float(sline[-2])
                grid_specs[key + "_max"] = float(sline[-1])

    return grid_specs

def read_regional_grid(grid_file, grid_loc=None):
    """Read the coordinates stored in a regional grid

    grid_file: str
        Grid .a or .b file name
    grid_loc: None or str
        Grid to read as one or several of "p", "q", "u" or "v"

    Return
    ------
    xarray.Dataset
        Dataset that contains the coordinates
    """

    allowed_grid_locs = "pquv"

    grid_a_file, grid_b_file = get_ab_file_names(grid_file)

    grid_specs = read_regional_grid_b(grid_b_file)
    rows = grid_specs["idm"]
    cols = grid_specs["jdm"]

    # Compute the number of bytes of a variable stored in a HYCOM .a
    # file for a grid of this size.

    varByteLength = int(math.ceil(float(cols * rows * 4) / 16384) * 16384)

    # Read the longitude and latitude grids as numpy arrays.

    coords = {}

    with open(grid_a_file, "rb") as f:

        for i, gl in enumerate(allowed_grid_locs):

            if grid_loc is not None and gl not in grid_loc:
                continue

            #            if grid_loc == 'p':
            #                varsToSkip = 0
            #            elif grid_loc == 'q':
            #                varsToSkip = 2
            #            elif grid_loc == 'u':
            #                varsToSkip = 4
            #            elif grid_loc == 'v':
            #                varsToSkip = 6

            f.seek(i * 2 * varByteLength)
            lon = np.ascontiguousarray(
                np.cast["float"](np.fromfile(f, ">f4", rows * cols)).reshape(
                    (cols, rows)
                )
            )
            f.seek((i * 2 + 1) * varByteLength)
            lat = np.ascontiguousarray(
                )
            )

            # I was told by Alan Wallcraft that the grid coordinates are only
            # accurate to four decimal digits (Fortran REAL*4 precision).
            # Round the coordinates to four digits.

            lon = lon.round(4)
            lat = lat.round(4)

            xname = gl + "lon"
            yname = gl + "lat"

            # Verifications
            assert (
                lon.min() == np.round(grid_specs[xname + "_min"], 4)
                and lon.max() == np.round(grid_specs[xname + "_max"], 4)
                and lat.min() == np.round(grid_specs[yname + "_min"], 4)
                and lat.max() == np.round(grid_specs[yname + "_max"], 4)
            ), "grid a and b files are inconsistant"

            if gl != "p":
                sn_suffix = "_at_{}_location".format(gl)
                ln_suffix = " at {} location".format(gl.upper())
            else:
                sn_suffix = ln_suffix = ""

            if grid_loc:
                xname = "lon"
                yname = "lat"

            coords[xname] = xr.DataArray(
                lon,
                name=xname,
                dims=("jdm", "idm"),
                attrs=dict(
                    long_name="Longitude" + ln_suffix,
                    standard_name="longitude" + sn_suffix,
                    units="degree_east",
                ),
            )
            coords[yname] = xr.DataArray(
                lat,
                name=yname,
                    long_name="Latitude" + ln_suffix,
                    standard_name="latitude" + sn_suffix,
                    units="degree_north",
                ),
            )

    # Create a dataset with coords only
    ds = xr.Dataset(coords=coords)
    ds.attrs.update(grid_specs)
    return ds

def read_regional_depth(depth_a_file, lon, lat):
    """Read the regional depth binary file

    depth_a_file: str
        Binary bathymetry file
    lon: xarray.DataArray
        Longitudes 2D array
    lat: xarray.DataArray

    Return
    ------
    xarray.DataArray
        Bathymetry
    """
    grid_a_file, _ = get_ab_file_names(depth_a_file)

    jdm, idm = lon.shape
    file = open(depth_a_file, mode="rb")
    field = np.reshape(
        np.fromfile(file, dtype="float32", count=idm * jdm).byteswap(),
        (jdm, idm),
    )
    file.close()
    field[field > 2 ** 99] = np.nan

    return xr.DataArray(
        field,
        name="bathymetry",
        dims=lon.dims,
        coords={"plon": lon, "plat": lat},
        attrs=dict(
            long_name="Bathymetry",
            units="m",
            standard_name="model_sea_floor_" "depth_below_sea_level",
        ),
    )

def write_res(ds, outfile_pattern="{field}.res"):
    """Write a Dataset to an unformatted fortran files

    One file per variable.

    ds: xarray.Dataset
    outfile_pattern: str
        File pattern that must contain ``{field}``

    Return
    ------
    dict
        Dict of lists of written files, where keys are fields

    """
    logger = logging.getLogger(__name__)
    cfs = get_cf_specs()

    # Coordinates
    lon = cfs.coords.get(ds, "lon").values.astype(">f4")
    lat = cfs.coords.get(ds, "lat").values.astype(">f4")
    if lon.ndim == 1:
        lon, lat = np.meshgrid(lon, lat)
    # lon = np.asfortranarray(lon.T)
    # lat = np.asfortranarray(lat.T)
    depth = cfs.get(ds, "depth").values.astype(">f4")
    time = (
        (cfs.coords.get(ds, "time").values - np.datetime64("1950-01-01"))
        / np.timedelta64(1, "D")
    ).astype(">f4")

    # Loop on time
    outfiles = {}

    for cf_var_name in ["temp", "sal", "ssh"]:

        var = cfs.get(ds, cf_var_name).fillna(-1e34)
        field = cfs.data_vars.get_name(cf_var_name, specialize=True)
        shape = list(var.shape[::-1])
        if var.ndim == 3:
            shape.insert(2, 1)
            dep = depth[0]
        else:
            dep = depth
        shape = np.array(shape, ">i4")

        outfile = outfile_pattern.format(**locals())
        with sio.FortranFile(outfile, "w", ">u4") as f:
            f.write_record(shape)
            f.write_record(lon)
            f.write_record(lat)
            f.write_record(dep)
            f.write_record(time)
            f.write_record(var.values.astype(">f4"))
        outfiles[field] = outfile
        logger.debug("Created .res file: " + outfile)

    return outfiles

def nc_to_res(
    ncfiles, outfile_pattern="{field}.res{ifile:03d}", preproc=None
):
    """Convert netcdf files to .res files

    There is one .res file per variable and per netcdf file.

    ncfiles: list
        List of netcdf files
        Pattern for .res file names. `{field}` refers to the
        Hycom variable name and `{ifile}` refers to the file index
        starting from 1.
    preproc: callable
        Callable that operates on the input datasets after reading it

    Return
    ------
    dict
        List of output files
    """
    resfiles = {}
    for ifile, ncfile in enumerate(ncfiles):
        ds = xr.open_dataset(ncfile)
        if preproc:
            ds = preproc(ds)
        field = "{field}"
        outfile_pattern = outfile_pattern.format(**locals())
        for vname, resfile in write_res(ds, outfile_pattern).items():
            resfiles.setdefault(vname, []).append(resfile)
        ds.close()
    return resfiles

def check_grid_dimensions(ddims, ds):
    """Check the consistency of grid dimensions between config and grid file

    conf: dict or vortex config
        Results from :func:`get_dimensions_dict_from_conf`
    ds: xarray.Dataset or str
        Result from :func:`read_regional_grid`
    """

    # # From conf
    # if not isinstance(ddims, dict):
    #     ddims = get_dimensions_dict_from_conf(ddims)

    # From file
    if isinstance(ds, str):
        ds = read_regional_grid(ds)

    # Check
    if "itdm0" not in ddims:
        ddims["itdm0"] = ds.plon.shape[1]
    if "jtdm0" not in ddims:
        ddims["jtdm0"] = ds.plon.shape[0]
    if (
        int(ddims["itdm0"]) != ds.plon.shape[1]
        or int(ddims["jtdm0"]) != ds.plon.shape[0]
    ):
        raise SloopCrocoError(
            "Dimensions inconsistency between "
            "config ({}, {}) "
            "and regional_grid file ({}, {}) ".format(
                ddims["itdm0"],
                ddims["jtdm0"],
                ds.plon.shape[1],
                ds.plon.shape[0],
            )
        )

def setup_stmt_fns(sigma, context):
    """Place the stmt_fns.h file depending on the context

    sigma: str, int
        Sigma value like 0 or 2
    context: {'ibc', 'model3d'}
    """
    sigma = str(sigma)
    if sigma not in CROCO_SIGMA_TO_STMT_FNS:
            "Invalid sigma value. Please choose one of: "
            + " ".join(CROCO_SIGMA_TO_STMT_FNS)
        )
    if context not in CROCO_CONTEXT_TO_STMT_FNS:
            "Invalid context name. Please choose one of: "
            + " ".join(CROCO_CONTEXT_TO_STMT_FNS)
        )

    stmt_fns_in = CROCO_SIGMA_TO_STMT_FNS[sigma]
    stmt_fns_out = CROCO_CONTEXT_TO_STMT_FNS[context]

    if os.path.exists(stmt_fns_out):
        os.remove(stmt_fns_out)
    shutil.copy(stmt_fns_in, stmt_fns_out)

def setup_gen_inicon(vconf):
    """Place the gen_inicon_hycom.F file depending on the vconf
    
    vconf: str
        Vconf value like manga, med, indien, gijon
    """
    gen_inicon_out = CROCO_IBC_GEN_INICON
    if vconf in CROCO_IBC_GEN_INICON_VCONF:
        gen_inicon_in = CROCO_IBC_GEN_INICON+f".{vconf}"
    else:
        gen_inicon_in = CROCO_IBC_GEN_INICON+".generic"
    logger.info(f"{gen_inicon_in} will be used as {gen_inicon_out}")
    if os.path.exists(gen_inicon_out):
        os.remove(gen_inicon_out)
    shutil.copy(gen_inicon_in, gen_inicon_out)
    
def setup_defstretch(vconf):
    """Place the defstretch.f/F file depending on the vconf
    
    vconf: Str
    """
    if vconf not in CROCO_MODEL_DEFSTRETCH_FILES.keys():
        vconf = "generic"
    defstretch_in = CROCO_MODEL_DEFSTRETCH_FILES[vconf]
    defstretch_out = defstretch_in.replace(f".{vconf}", "")
    logger.info(f"{defstretch_in} will be used as {defstretch_out}")
    defstretch_path = os.path.dirname(defstretch_out)
    for fmt in [".F", ".f"]:
        defstretch_file = os.path.join(defstretch_path,
                                       "defstretch"+fmt)
        if os.path.exists(defstretch_file):
            os.remove(defstretch_file)
    shutil.copy(defstretch_in, defstretch_out)
    
RE_SCALAR_DECLARATION_MATCH = re.compile(
    r"^\s*(\S+)\s*'(\w+)\s*'\s+= (.+)$"
).match

def parse_scalar_declaration(line):
    """Parse a string that declare a scalar

    Format: ``  720      'idm   ' = longitudinal array size``

    line: str

    Return
    ------
    tuple
        ``(name, value, description)`` where ``value`` is converted to
        its numerical form with eval.
    """
    groups = RE_SCALAR_DECLARATION_MATCH(line).groups()
    return groups[1], eval(groups[0]), groups[2]

def read_blkdat_input(blkdat_input_file):
    """Read and parse the blkdat.input file

    blkdat_input_file: str

    Return
    ------
        All scalars are stored to a 0D array with a long_name attribute.
        Sigma value is stored as an 1D array of dim "lev" = kdm.
    """

    # Parse the file and get a dict

    content = {}
    sigmas = []
    header = ""
    with open(blkdat_input_file) as f:
        for i in range(3):
            header += f.readline()[:-1]
        f.readline()
            if line:
                name, value, long_name = parse_scalar_declaration(line[:-1])
                if name == "sigma":
                    sigmas.append(value)
                else:
                    content[name] = xr.DataArray(
                        value, name=name, attrs={"long_name": long_name}
                    )

    # Sigma as array
    nlayer = len(sigmas)
    content["sigma"] = xr.DataArray(
        sigmas,
        coords=[
            ("lev", np.arange(1, nlayer + 1), {"long_name": "Layer level"})
        ],
    )

    return xr.Dataset(content, attrs={"header": header})

def get_restart_date(ncname):
     ds=xr.open_dataset(ncname,decode_times=False)
     units, reference_date = ds.time.attrs['units'].split('since')
     time=pd.to_datetime(ds.scrum_time.values,unit='s',origin=pd.Timestamp(reference_date.strip()))
     return time[0].strftime('%Y%m%dT%H%M')

    """Format variable in a dataset or dataarray to be hycom compliant"""
        ds[time.name].encoding["units"] = "days since 1950-01-01"

        for attr in "dtype", "scale_factor", "add_offset":
            ds[vname].encoding["_FillValue"] = -1e34

    return ds

def bin2hycom(
    resfile, nx, ny, nz, resdates, positive, fileab, header, field
):
    """Convert .res files to .a,.b files

    This does the job of program `bin2hycom`

    resfile: str
    nx: int
    ny: int
    nz: int
    resdates: array_like
    fileab: str
    header: str
    field: str
    positive: bool
    """
    # Inits
    n2drec = int(((nx * ny + 4095) // 4096) * 4096)
    field1d = np.zeros(n2drec, dtype=">f4")
    filea, fileb = get_ab_file_names(fileab)

    # Read dates
    with open(resdates, "rb") as fd:
        nt = int(np.fromfile(fd, ">i4", 1))
        fd.seek(4)
        dates = np.empty(nt, dtype=">f4")
        for it in range(nt):
            fd.seek(4*(it+1))
            dates[it] = np.fromfile(fd, ">f4", 1)

    # Write header
    fb = open(fileb, "w")
    for line in header, "", "", "", f"i/jdm =  {nx} {ny}":
        fb.write("{:<79}\n".format(line))

    # Loop on .res file records
    with open(resfile, "rb") as fr, open(filea, "wb") as fa:
        for it, date in enumerate(dates):
            for iz in range(nz):
                irec = it * nz + iz

                # Read .res
                fr.seek(irec * nx * ny * 8)
                field2d = np.fromfile(fr, ">f8", nx * ny)

                # Write to .a
                field1d[: nx * ny] = field2d.ravel()
                if positive:
                    field1d[field1d < 0] = 0
                fa.seek(irec * n2drec * 4)
                field1d.tofile(fa)

                # Write to .b
                vmin = np.nanmin(field2d)
                vmax = np.nanmax(field2d)
                # (a21,1x,E18.8,1x,E18.8,1x,E18.8)
                fb.write(
                    f"{field:>7}: date,range = "
                    f"{date:18e} {vmin:18e} {vmax:18e}\n"
                )

def run_bin2hycom(
    nx,
    ny,
    nz,
    inpat="{field}.res",
    resdates="time_values.res",
    outpat="{short_name}-nest",
):
    """Run the :func:`binhycom` function with the default file names

    Specifications are provided by the :data:`BYN2HYCOM_SPECS` content.

    inpat: str
        Input .res file name pattern with "{field}" pattern inside.
    resdates: str
        Input .res file that contains dates.
    outpat: str
        Output .a and .b file name patterns with "{short_name}
        pattern inside and without extension.

    Return
    ------
    list
        The effective list of output files.
    """
    outfiles = []
    for field, specs in BIN2HYCOM_SPECS.items():
        resfile = inpat.format(**locals())
        short_name = field[:1]
        fileab = outpat.format(**locals())
        bin2hycom(
            resfile,
            nx,
            ny,
            nz,
            resdates,
            specs["positive"],
            fileab,
            specs["header"],
            field,
        )

        outfiles.append(fileab + ".a")
        outfiles.append(fileab + ".b")

def rest_head(date, outfile="rest_new_head", xcycles=1):
    """Write a binary date file for restart

    ==========
    date:
        A pandas compatible date that is rounded to midnight.
        It is typically at J-1.
    outfile: str
        Output file
    xcycles: int
        1 for restart from an external source else > 1
    """
    date = pd.to_datetime(date).floor("D")
    date = (date - pd.to_datetime("1950-01-01")) / pd.to_timedelta("1D")
    with open(outfile, "wb") as f:
        np.array([xcycles, date], dtype=">f8").tofile(f)

def read_patch_input(patch_input):
    """Read the patch.input MPI tiling file

    patch_input: str

    Return
    ------
        2D arrays contains the start and end indices along i
        and 1D arrays contains the start and end indices along j.
        Other numbers are stored as global attributes.
    """
    with open(patch_input) as f:

        # Read sizes
        names = f.readline()[:-1].split()
        numbers = [int(n) for n in f.readline()[:-1].split()]
        sizes = dict(zip(names, numbers))
        print(sizes)

        # Init arrays
        ispt = xr.DataArray(
            np.zeros((sizes["mpe"], sizes["npe"]), dtype='i'),
            dims=("mpe", "npe"),
            attrs={"long_name": "i index start"})
        iipe = ispt.copy()
        iipe.attrs.update(long_name="i index end")
        jspt = xr.DataArray(
            np.zeros(sizes["mpe"], dtype='i'), dims="mpe",
            attrs={"long_name": "j index start"})
        jjpe = jspt.copy()
        jjpe.attrs.update(long_name="j index end")

        # Read i indices
        for j in range(sizes["mpe"]):
            ispt[j] = f.readline()[11:-1].split()
            iipe[j] = f.readline()[11:-1].split()

        # Read j indices
        jspt[:] = f.readline()[11:-1].split()
        jjpe[:] = f.readline()[11:-1].split()

        # Create dataset
        ds = xr.Dataset(
            {"ispt": ispt, "iipe": iipe, "jspt": jspt, "jjpe": jjpe})
        for name in set(names).difference({"npe", "mpe"}):
            ds.attrs[name] = sizes[name]
        return ds

def hycom2nc(nx, ny, nz,
             grid_file,
             hycom_file,
             outpat="{field}-nest.nc",
    ):
    """func hycomnc function 
    
        Output nc file name

    Return
    ------
    str
        The output file
    """
     
    ds = read_regional_grid(grid_file, grid_loc='p')
    depth = np.arange(0, nz, 1)
    filea, fileb = get_ab_file_names(hycom_file)
    with open(fileb, 'r') as fb:
        series = fb.readlines()[5:]
    time = []
    for itime, iseries in enumerate(np.arange(0, len(series), nz)):
        matches  = re.search(r"\s*(?P<field>\w*)\W\s*date,range\s=\s*(?P<time>\d.\d*\w\W\d*)\s*(?P<min>\d.\d*\w\W\d*)\s*(?P<max>\d.\d*\w\W\d*)", 
                             series[iseries], re.DOTALL)
        time.append(convert_from_julian_day(float(matches['time'])*86400, unit='s'))
        field = matches['field']
    field_values = np.zeros((len(time), nz, ny, nx))
    fa = open(filea, 'r')
    for it in range(len(time)):
        for iz in range(len(depth)):
            irec = it * nz + iz
            fa.seek(irec * n2drec * 4)
            field_values[it, iz, :, :] = np.reshape(np.fromfile(fa, ">f4", nx * ny), (ny,nx))
    fa.close()
    da = xr.DataArray(field_values, name=field, 
                      dims=('time', 'lev', 'y', 'x'),
                      )
    da = da.where(da!=0)
    da.encoding["_FillValue"] = -1.e34
    ds = xr.merge([ds, da])
    ds = ds.assign_coords(time=time, depth=depth)
    ds = format_ds(ds) 
    ncout = outpat.format(**locals())
    ds.to_netcdf(ncout)
    return ncout

def nc2hycom(varname, nc_in, header):
    varname = varname[0] if len(varname)>1 else varname
    ds = xr.open_dataset(nc_in)
    ds = ds.fillna(0)
    nt, nz, ny, nx = np.shape(ds[varname])
    
    filea, fileb = get_ab_file_names(nc_in.split('.')[0])

     # Write header

    with open(filea, "wb") as fa:
        for it, date in enumerate(ds.time.values):
                field1d[: nx * ny] = ds[varname][it, iz, :, :].values.ravel()
                date = convert_to_julian_day(ds.time.values[it]) 
                vmin = ds[varname][it, iz, :, :].values.min()
                vmax = ds[varname][it, iz, :, :].values.max()
                fb.write(
                    f"{varname:>7}: date,range = "
                )
    fb.close()

"""
Open Boundary condition Forcing preprocessings
"""

import os

)

obc_models={
    "glorys12v1": {
        "data_dir":"/home/ref-ocean-reanalysis/global-reanalysis-phy-001-030-daily/",
        "prefix":"mercatorglorys12v1_gl12_mean_",
        "suff":".nc",
        "freq":"1d",
        "temp_in_kelvin":False,
            'depth':"depth",
            'longitude':'nav_lon',
            'latitude':'nav_lat',
            'time':'time_counter',
            'XE':'zos',
            'TEMP':'thetao',
            'SAL':'so',
            'UZ':'uo',
            'VZ':'vo',
        }
    }
}
def rename_vars(ds,mod):
    """Rename variables as expected by Extract"""
    LOGGER.info("Rename obc variables as expected by Extract program")
    obc_names=mod['varnames']
    list_nc_vars=list(ds.variables.keys())
    list_nc_dims=list(ds.dims.keys())
    for variable in list_nc_vars:
      for key,value in obc_names.items():
        if variable == value and value != key:
          ds = ds.rename_vars({variable: key})
          if variable in list_nc_dims:
           print("variable also a dim")
           ds = ds.rename_dims({variable: key})
    for dsi in list(ds.variables.keys()): 
       if dsi not in obc_names: 
    return ds

def format_obc_to_extract(ds,model):
    """Format OBC file for the Extract program"""
    LOGGER.info("Format to Extract")
    mod=obc_models[model]
    ds=rename_vars(ds,mod)
    #change temperature units if necessary
    if mod['temp_in_kelvin']:
        ds.variables['TEMP']=ds.variables['TEMP']-273.15
    #set time in seconds
    ds.time.encoding.update({'units':"seconds since 1900-01-01",'dtype':float})
    # unpack data
    LOGGER.info('unpack and set valid min & max')
    for var in list(ds.variables.keys()):
      if var not in list(ds.dims.keys()):
        ds[var].encoding.update({'_FillValue':-999.,'scale_factor':1.0,'add_offset':0.0,'dtype':float})
        #set_validminmax(ds[var])
    return ds

def reset_validminmax(fname):
   """reset valid min and max for a dataset"""
   LOGGER.info('reset valid min & max')
   ds=xr.open_mfdataset(fname)
   ds2=ds.copy()
   ds.close()
   for var in list(ds2.variables.keys()):
      if var not in list(ds2.dims.keys()): 
       set_validminmax(ds2[var])
   os.remove(fname)
   ds2.to_netcdf(fname,mode='w')
     
def set_validminmax(data):
    """set valid min & max of data variable"""
    vmin=float(data.min().to_numpy())
    vmax=float(data.max().to_numpy())
    data.attrs.update(valid_min=vmin,valid_max=vmax)
    return data

def get_dates_from_file(fname,tname="time"):
    print("open",fname)
    ds=xr.open_mfdataset(fname)
    time=ds.indexes[tname]
    d1=time[0].strftime('%Y%m%dT%H%M')
    d2=time[-1].strftime('%Y%m%dT%H%M')
    return d1,d2

def build_obc_list(date1,date2,model):
    """Build obc list of files from OGCM model"""
    LOGGER.info("Get list of obc files")
    suff=mod['suff']
    inst_list=[f"{rac}{date.strftime('%Y%m%d')}" for date in dates_list]

    inst_list=[glob.glob(f"{data_dir}/*/*/{file}*.nc")[0] for file in inst_list]

    return inst_list

def read_obc_data(inst_list,lon_bnds,lat_bnds):
    """Read the list of obc files and reduce it to area of interest"""
    return ds_inst

"""
HYCOM utilities

"""
import os

class SloopHycomError(SloopError):
    pass

#: HYCOM vconf
HYCOM3D_VCONF = ["manga", "med", "indien", "gijon"]

#: HYCOM platform root dir
HYCOM3D_DIR = os.path.join(SLOOP_SRC_DIR, 'HYCOM')

#: HYCOM initial and open boundary root dir
HYCOM3D_IBC_DIR = os.path.join(HYCOM3D_DIR, 'ChaineFrontieresOuvertes')

#: HYCOM gen_inicon_hycom.F file
HYCOM3D_IBC_GEN_INICON = os.path.join(
    HYCOM3D_IBC_DIR, 'INTERP-VER/gen_inicon_hycom.F')

#: HYCOM vconf which needs particular gen_inicon_hycom.F file
HYCOM3D_IBC_GEN_INICON_VCONF = ["med", "indien"]

#: HYCOM inicon executable
HYCOM3D_IBC_INICON = os.path.join(HYCOM3D_IBC_DIR, 'INTERP-VER/inicon')

#: HYCOM regridcdf executable
HYCOM3D_IBC_REGRIDCDF = os.path.join(HYCOM3D_IBC_DIR, 'INTERP-HOR/regridcdf')

#: HYCOM 3D model directory
HYCOM3D_MODEL_DIR = os.path.join(HYCOM3D_DIR, 'Modele3D')

#: HYCOM 3D model executable
HYCOM3D_MODEL_EXE = os.path.join(
    HYCOM3D_MODEL_DIR, 'MAKELIB.', 'main')

#: HYCOM 3D model dimensions.h.template
HYCOM3D_MODEL_DIMENSIONSH_TEMPLATE = os.path.join(
    HYCOM3D_MODEL_DIR, 'COMMUNS./dimensions.h.template')

#: HYCOM rivers directory
HYCOM3D_RIVERS_DIR = os.path.join(HYCOM3D_DIR, 'Rivers')

#: HYCOM rivers flow climatology
HYCOM3D_RIVERS_FLOW_CLIM = os.path.join(HYCOM3D_RIVERS_DIR, "hycom_rivers.d")

#: HYCOM rivers hydro climatology
HYCOM3D_RIVERS_HYDRO_CLIM = os.path.join(
    HYCOM3D_RIVERS_DIR, "fleuves_climato_hydro")

#: HYCOM 3D model defstretch.f/F
HYCOM3D_MODEL_DEFSTRETCH_FILES = {
    "generic": os.path.join(
        HYCOM3D_MODEL_DIR, 'UTIL./defstretch.f.generic'),
    "med": os.path.join(
        HYCOM3D_MODEL_DIR, 'UTIL./defstretch.F.med')
    }

#: HYCOM spnudge root dir
HYCOM3D_SPNUDGE_DIR = os.path.join(HYCOM3D_DIR, 'Assimilation', 'Spectral_Nudging')

#: HYCOM demerliac executable
HYCOM3D_SPNUDGE_DEMERLIAC = os.path.join(
    HYCOM3D_SPNUDGE_DIR, 'demerliac')

#: HYCOM filter executable
HYCOM3D_SPNUDGE_SPECTRAL = os.path.join(
    HYCOM3D_SPNUDGE_DIR, 'lissage_spectral_fft')

#: HYCOM postproduction root dir
HYCOM3D_POSTPROD_DIR = os.path.join(HYCOM3D_DIR, 'PostProduction')

#: HYCOM postproduction time filtering executable
HYCOM3D_POSTPROD_FILTER = os.path.join(
    HYCOM3D_POSTPROD_DIR, 'cdfTraitementTemporel3D')

#: HYCOM postproduction file conversion executable
HYCOM3D_POSTPROD_INTERPOLATION = os.path.join(
    HYCOM3D_POSTPROD_DIR, 'cdfHycomTodataSHOMouSOAP')

#: HYCOM postproduction temperature conversion executable
HYCOM3D_POSTPROD_TEMPCONVERSION = os.path.join(
    HYCOM3D_POSTPROD_DIR, 'temppot2tempis_postProd')

#: HYCOM postproduction traductions noms longs
HYCOM3D_POSTPROD_LONGNAMES = "PARAMETERS{rank}./traductions_noms_longs_{postprod}"

#: HYCOM postproduction traductions noms standards
HYCOM3D_POSTPROD_STDNAMES = "PARAMETERS{rank}./traductions_noms_standards_2"

#: HYCOM postproduction traductions units
HYCOM3D_POSTPROD_UNITS = "PARAMETERS{rank}./traductions_unites"

#: HYCOM sigma to stmt_fns*.h files
HYCOM3D_SIGMA_TO_STMT_FNS = {
    '0': os.path.join(HYCOM3D_MODEL_DIR, 'COMMUNS./stmt_fns_SIGMA0_7term.h'),
    '2': os.path.join(HYCOM3D_MODEL_DIR, 'COMMUNS./stmt_fns_SIGMA2_17term.h'),
}

#: HYCOM mapping for stmt_fns.h files
HYCOM3D_CONTEXT_TO_STMT_FNS = {
    'ibc_hor': os.path.join(HYCOM3D_IBC_DIR, 'INTERP-HOR/stmt_fns.h'),
    'ibc_ver': os.path.join(HYCOM3D_IBC_DIR, 'INTERP-VER/stmt_fns.h'),
    'model': os.path.join(HYCOM3D_MODEL_DIR, 'COMMUNS./stmt_fns.h')
    }

#: HYCOM mask file
HYCOM3D_MASK_FILE = 'FORCING{rank}./regional.mask.nc'

#: HYCOM grid a file
HYCOM3D_GRID_AFILE = 'FORCING{rank}./regional.grid.a'

#: HYCOM grid b file
HYCOM3D_GRID_BFILE = 'FORCING{rank}./regional.grid.b'

#: HYCOM depth a file
HYCOM3D_DEPTH_AFILE = "FORCING{rank}./regional.depth.a"

#: HYCOM depth b file
HYCOM3D_DEPTH_BFILE = "FORCING{rank}./regional.depth.b"

#: HYCOM isosigma a file
HYCOM3D_ISOSIGMA_AFILE = 'FORCING{rank}./iso.sigma.a'

#: HYCOM isosigma b file
HYCOM3D_ISOSIGMA_BFILE = 'FORCING{rank}./iso.sigma.b'

#: HYCOM regricdf tabindicesT
HYCOM3D_REGRIDCDF_TABINDICEST = 'FORCING{rank}./tabindicesT'

#: HYCOM blkdat.input parameters file
HYCOM3D_BLKDAT_INPUT_FILE = "PARAMETERS{rank}./blkdat.input"

HYCOM3D_BLKDAT_CMO_INPUT_FILE = "PARAMETERS{rank}./blkdat_cmo.input"

#: HYCOM blkdat_cmo.input parameters file
HYCOM3D_BLKDAT_CMO_INPUT_FILES = {
            "spinup":"PARAMETERS{rank}./blkdat_cmo_A.input",
            "forecast":"PARAMETERS{rank}./blkdat_cmo_P.input",
            "spnudge_free":"PARAMETERS{rank}./blkdat_cmo_A.input",
            "spnudge_relax":"PARAMETERS{rank}./blkdat_cmo_Assim.input",
            "free": "PARAMETERS{rank}./blkdat_cmo_A.input",
            }

#: HYCOM savefield.input parameters file
HYCOM3D_SAVEFIELD_INPUT_FILE = "PARAMETERS{rank}./savefield.input"

HYCOM3D_SAVEFIELD_INPUT_FILES = {
            }

#: HYCOM defstrech.input parameters file
HYCOM3D_DEFSTRECH_INPUT_FILE = "PARAMETERS{rank}./defstrech*.input"

#: HYCOM MPI tiling file
HYCOM3D_PATCH_INPUT_FILE = "SPLIT{rank}./patch.input"

#: HYCOM ports.input parameters file
HYCOM3D_PORTS_INPUT_FILE = "PARAMETERS{rank}./ports.input"

#: HYCOM run.input.tpl parameters file
HYCOM3D_RUN_INPUT_TPL_FILE = "PARAMETERS{rank}./run.input.tpl"

#: HYCOM post-production configuration files
HYCOM3D_POSTPROD_CFG_FILE = "PARAMETERS{rank}./postprod.cfg"

#: HYCOM post-production configuration specs file
HYCOM3D_POSTPROD_INI_FILE = os.path.join(_THISDIR, "postprod.ini")

#: HYCOM post-production configuration output files
HYCOM3D_POSTPROD_CONFIG_FILE = "PARAMETERS{rank}./postprod_{postprod}.cfg"

#: HYCOM check promises
HYCOM3D_PROMISES = "PARAMETERS{rank}./promises.cfg"

#: HYCOM rivers configuration specs file
HYCOM3D_RIVERS_INI = os.path.join(_THISDIR, "rivers.ini")

#: HYCOM rivers configuration file
HYCOM3D_RIVERS_CFG = 'FORCING{rank}./nest/rivers.cfg'

#: HYCOM spectral nudging parameters json file
HYCOM3D_SPNUDGE_JSON_FILE = "PARAMETERS{rank}./spnudging_parameters.json"

HYCOM3D_STANDARD_NAME_TO_NAME = {
    }

        ],
    # Model
        ],
        ],
    }

        ],
        ],
        ],
    }

HYCOM3D_GEO_CONSTANT_PATTERNS = {
    "nest": [
        "FORCING{rank}./nest/debitmeso.list",
        HYCOM3D_RIVERS_CFG,
        "FORCING{rank}./nest/rivers.list",
        "FORCING{rank}./nest/rmu.[ab]"
    ],
    "tide": [
        "FORCING{rank}./tide/tide.list",
        "FORCING{rank}./tide/forcing.*tide.[ab]"
    ],
    "regional": [
        HYCOM3D_GRID_AFILE,
        HYCOM3D_GRID_BFILE,
        HYCOM3D_DEPTH_AFILE,
        HYCOM3D_DEPTH_BFILE,
        HYCOM3D_MASK_FILE,
        HYCOM3D_REGRIDCDF_TABINDICEST,
        HYCOM3D_ISOSIGMA_AFILE,
        HYCOM3D_ISOSIGMA_BFILE
    ],
    "run": ([
        HYCOM3D_BLKDAT_INPUT_FILE,
        HYCOM3D_BLKDAT_CMO_INPUT_FILE,
        HYCOM3D_DEFSTRECH_INPUT_FILE,
        HYCOM3D_PORTS_INPUT_FILE,
        HYCOM3D_RUN_INPUT_TPL_FILE,
        HYCOM3D_SPNUDGE_JSON_FILE
        ]+list(HYCOM3D_BLKDAT_CMO_INPUT_FILES.values())),
        HYCOM3D_POSTPROD_CFG_FILE,
        HYCOM3D_POSTPROD_CONFIG_FILE.format(
    ],
        [HYCOM3D_SAVEFIELD_INPUT_FILE] +
        list(set(p for pp in HYCOM3D_SAVEFIELD_INPUT_FILES.values()
    ),
    "split": [
        HYCOM3D_PATCH_INPUT_FILE
    ],
        HYCOM3D_PROMISES
    ],
}

HYCOM3D_NOT_GEO_CONSTANT_PATTERNS = {
    ]
}

#: Hycom3d specifications for the :mod:`xoa.cf` module use
HYCOM3D_CF_SPECS = {
        "temp": {
        },
        "sal": {
            "name": "saln"
        },
        "ssh": {
            "name": "ssh"
        },
        },
        "ps": {
            "attrs": {
            }
        },
        "tair": {
            "attrs": {
                ],
            }
        },

    }
}

HYCOM3D_PARAM = {
    "preatm": {
        "name": [
            "preatm",
            "prmsl"
        ],
        "attrs": {
        }
    },
    "airtmp": {
        "name": [
            "airtmp",
            "t2m"
        ],
        "attrs": {
            ],
        }
    },
    "vapmix": {
        "name": [
            "vapmix"
        ],
        "attrs": {
            "units": "kg kg-1",
                "water_vapor_mixing_ratio"
            ],
            "long_name": "water vapor mixing ratio",
            "header": "vapeur"
        }
    },
    "shwflx": {
        "name": [
            "shwflx",
            "ssr"
        ],
        "attrs": {
            ],
        }
    },
    "radflx": {
        "name": [
            "radflx",
            "str",
       ],
        "attrs": {
            ],
        }
    },
    "precip": {
        "name": [
            "precip",
            "tp"
        ],
        "attrs": {
            ],
        }
    },
    "wndspd": {
        "name": [
            "wndspd",
            "si10"
        ],
        "attrs": {
            ],
        }
    },
    "tauewd": {
        "name": [
            "tauewd",
            "ewss"
        ],
        "attrs": {
            "units": "N m-2",
                "zonal_wind_stress"
            ],
            "long_name": "zonal wind stress",
            "header": "vent - composante zonale"
        }
    },
    "taunwd": {
        "name": [
            "taunwd",
            "nsss"
        ],
        "attrs": {
                "meridional_wind_stress"
            ],
            "long_name": "meridional wind stress",
            "header": "vent - compostante meridienne"
        }
    }
}

_CACHE = {}

    """Get dict of HYCOM3D constant file glob patterns for each pack

    geo: bool
    """
    return HYCOM3D_GEO_CONSTANT_PATTERNS if geo else HYCOM3D_NOT_GEO_CONSTANT_PATTERNS

    """Get a :class:`xoa.cf.CFSpecs` object for Hycom3D"""
        _CACHE["cf_specs"] = xoa.cf.CFSpecs(HYCOM3D_CF_SPECS)

    """Format all variables in a dataset or dataarray to be hycom compliant"""

    return ds

"""
"""

)
from .__init__ import HYCOM3D_PARAM

    ):
        """

        """

        dsc = []
            ds = []
                    )
                )
            )
            )
        )
        dsi = []
            ds = []
                        },
                        ],
                    )
                )
        )
        )

        else:

        for pv in HYCOM3D_PARAM.values():
        )
         
            "t2m",
            "r2",
            "prmsl",
            "ssr",
            "tp",
            "si10",
            "vapmix",
            "radflx",
        ]
        )

        for hyparam in list(HYCOM3D_PARAM.keys()):
            if len(HYCOM3D_PARAM[hyparam]["name"]) > 1:
                HYCOM3D_PARAM[hyparam]["name"].remove(hyparam)
                if HYCOM3D_PARAM[hyparam]["name"][0] in ds:
                    ds = ds.rename({HYCOM3D_PARAM[hyparam]["name"][0]: hyparam})
            ds[hyparam].attrs = HYCOM3D_PARAM[hyparam]["attrs"]
            if dsi not in list(HYCOM3D_PARAM.keys()):
        return ds

        )
        )
                            fa
                        )
                        )
                        )
                        )

"""

"""

import re
import os

    SloopHycomError,
    HYCOM3D_SIGMA_TO_STMT_FNS,
    HYCOM3D_CONTEXT_TO_STMT_FNS,
    HYCOM3D_IBC_GEN_INICON,
    HYCOM3D_IBC_GEN_INICON_VCONF,
    HYCOM3D_MODEL_DEFSTRETCH_FILES
)

    "h": {
    },
    "temp": {
    },
    "saln": {
    },
    "u": {
    },
    "v": {
    },
}

    path: str

    Return
    ------
    str: afile
    str: bfile

    """
    else:
        root = path

    Return
    ------
    dict
    """

            if line:

    Return
    ------
    """

    coords = {}

                continue

                )
            )
                )
            )

            assert (

            else:

                lon,
                ),
            )
                lat,
                ),
            )

    return ds

    Return
    ------
    """

    )

        field,
        ),
    )

    Return
    ------
    dict

    """

    time = (

        else:
            dep = depth

):

    Return
    ------
    dict
    """

    """

    # From file

    # Check
    if (
    ):
        raise SloopHycomError(
            )
        )

    """
    if sigma not in HYCOM3D_SIGMA_TO_STMT_FNS:
            + " ".join(HYCOM3D_SIGMA_TO_STMT_FNS)
        )
    if context not in HYCOM3D_CONTEXT_TO_STMT_FNS:
            + " ".join(HYCOM3D_CONTEXT_TO_STMT_FNS)
        )

    stmt_fns_in = HYCOM3D_SIGMA_TO_STMT_FNS[sigma]
    stmt_fns_out = HYCOM3D_CONTEXT_TO_STMT_FNS[context]

    vconf: str
    """
    gen_inicon_out = HYCOM3D_IBC_GEN_INICON
    if vconf in HYCOM3D_IBC_GEN_INICON_VCONF:
        gen_inicon_in = HYCOM3D_IBC_GEN_INICON+f".{vconf}"
    else:
        gen_inicon_in = HYCOM3D_IBC_GEN_INICON+".generic"
    
    vconf: Str
    """
    if vconf not in HYCOM3D_MODEL_DEFSTRETCH_FILES.keys():
    defstretch_in = HYCOM3D_MODEL_DEFSTRETCH_FILES[vconf]
    
).match

    line: str

    Return
    ------
    tuple
    """

    Return
    ------
    """

    sigmas = []
    header = ""
            if line:
                else:
                    )

        sigmas,
        coords=[
        ],
    )

    return ds

):

    nx: int
    ny: int
    nz: int
    field: str
    """
    # Inits

                # Read .res

                fb.write(
                )

    nx,
    ny,
    nz,
):

    inpat: str

    Return
    ------
    list
    """
            resfile,
            nx,
            ny,
            nz,
            resdates,
            fileab,
            field,
        )

    date:
    """

    Return
    ------
    """

        return ds

    ):
    
    Return
    ------
    str
    """
     
    time = []
                      )

                fb.write(
                )

"""
Hycom3d post-production
"""
import xarray as xr
import xesmf
import os, shutil, re, logging
import multiprocessing as mp

from ...interp import Regridder

FILTER_map = {
    "mean":"1",
    "demerliac":"2",
    "godin":"3"
    }

def preproc_time_splitting(ds):
    return ds.isel({"time":slice(None,-1)})

def postproc_time_splitting(ds):
    
            if vname=="bathymetry":
                if "time" in ds[vname].dims:
                    ds[vname] = ds[vname].isel(time=0)
    return ds

def get_cmdargs_nc_time_splitting(ncins, fields, ncout="{field}_{date}.nc", freq="PT1H"):
    """Split a netcdf file into several netcdf files at a specific frequency"""
    cmdargs_list = []
    for field in fields:
        for ncin in ncins:
            if ncin.startswith(field):
                if ("spinup" in ncin) or ("spnudge1" in ncin):
                    cmdargs = [ncin, freq, ncout.format(field=field, date="{date}"), preproc_time_splitting, postproc_time_splitting]
                else:
                    cmdargs = [ncin, freq, ncout.format(field=field, date="{date}"), None, postproc_time_splitting]
                cmdargs_list.append(cmdargs)
    return cmdargs_list

def write_verticalgrid(depths, filename="verticalgrid.txt"):
    """Write the text file containing the depth for vertically interpolation"""
    
    f = open(filename, "w+")
    f.write("{}\n".format(int(len(depths))))
    for d in depths:
        f.write(f"{d}\n")
    f.close()
    return filename

def get_fields_from_cfg(cfg_var, dim=None):
    """Get the hycom3d variable names with respect to post-prod name"""

    hycom3d_var = []
    areas = cfg_var.keys()
    areas.remove("common")
    srcs = []
    for area in areas:
        src = cfg_var[area]["src"]
        if src not in srcs:
            srcs.append(src)
    for src in srcs:
        hycom3d_var.append(f"h_3D_{src}")
        if (dim=="2D") or (dim==None):
            if "none" not in cfg_var[area]["variables"]["2D"]:
                if "ssh" in cfg_var[area]["variables"]["2D"]:
                    hycom3d_var.append(f"ssh_2D_{src}")
                if "sst" in cfg_var[area]["variables"]["2D"]:
                    hycom3d_var.append(f"sst_2D_{src}")
                if "sss" in cfg_var[area]["variables"]["2D"]:
                    hycom3d_var.append(f"sss_2D_{src}")
                if "u" in cfg_var[area]["variables"]["2D"]:
                    hycom3d_var.append(f"u_2D_{src}")
                if "v" in cfg_var[area]["variables"]["2D"]:
                    hycom3d_var.append(f"v_2D_{src}")
        if (dim=="3D") or (dim==None):            
            if "none" not in cfg_var[area]["variables"]["3D"]:
                for var in cfg_var[area]["variables"]["3D"]:
                    hycom3d_var.append(f"{var}_3D_{src}")
    return hycom3d_var
  
def get_postprod_term(specs, filtering=None):
    """Get the list of post-production dates from configuration"""
    
    bounds_term = specs["common"]["bounds_term"]
    start = bounds_term["min"]
    end = bounds_term["max"]
    term_postprod = np.arange(start, end+24, 24)
    if (filtering=="demerliac") or (filtering=="godin"):
        return term_postprod[1:-2]
    elif filtering=="mean":
        return term_postprod[:-1]
    else:
        return term_postprod

def get_filter_clargs(ncins, specs, ncout_patt):
    """Get opts for filtering executable"""
    
    terms_postprod = get_postprod_term(specs)
    opts_list = []
    srcs = []
    areas = specs.keys()
        src = specs[area]["src"]
        for dim in ["3D", "2D"]:
            for filter in specs[area]["filters"][dim]:
                if not filter=="none":
                    fields = specs[area]["variables"][dim]
                    if dim=="3D":
                            fields.append("h")
                    elif dim=="2D":
                        if filter not in specs[area]["filters"]["3D"]:
                    for term_postprod in get_postprod_term(specs, filtering=filter):
                        term_before = term_postprod-24
                        term_after = term_postprod+24
                        for field in fields:
                            if field=="h":
                                filling = "3D"
                            else:
                                filling = dim
                            if not field=="none":
                                opts = {"filter": FILTER_map[filter],
                                    "ncout": f"{field}_{filling}_{src}_{term_postprod:04d}_{filter}.nc" if term_postprod>=0 else f"{field}_{filling}_{src}_{term_postprod:05d}_{filter}.nc",
                                    "ncins":{
                                        "ncin_for": f"{field}_{filling}_{src}_{term_after:04d}.nc" if term_after>=0 else f"{field}_{filling}_{src}_{term_after:05d}.nc",
                                        "ncin_mid": f"{field}_{filling}_{src}_{term_postprod:04d}.nc" if term_postprod>=0 else f"{field}_{filling}_{src}_{term_postprod:05d}.nc",
                                        "ncin_back": f"{field}_{filling}_{src}_{term_before:04d}.nc" if term_before>=0 else f"{field}_{filling}_{src}_{term_before:05d}.nc"
                                        }
                                }
                            if filter=="mean":
                                opts["ncins"].pop("ncin_back")
                            ncin_filter = [val for key, val in opts["ncins"].items() if key.startswith("ncin")]
                            if set(ncin_filter).issubset(ncins):
                                opts_list.append(opts)
    return opts_list

def get_interpolation_clargs(ncins, offset, specs, config, verticalgrid):
    """Get command opts for vertical interpolation executable"""
   
    srcs = []
                fields = specs[area]["variables"][dim]
                terms_postprod = get_postprod_term(specs, filtering=filter)
                for term_postprod in terms_postprod:
                    prec = "04d" if term_postprod >= 0 else "05d"
                    for field in fields:
                        if not field=="none":
                            opts = {"offset": offset,
                                    "ncins": {"ncin": f"{field}_{dim}_{src}_{term_postprod:{prec}}_{filter}.nc",
                                               "h": f"h_3D_{src}_{term_postprod:{prec}}_{filter}.nc"
                                               },
                                    "ncout": f"{field}_{dim}_{src}_{term_postprod:{prec}}_{filter}_interp.nc",
                                    "config": os.path.basename(config),
                                    "zgrid": verticalgrid
                                    }
                            if dim=="2D":
                                opts.pop("zgrid")
                            if set(opts["ncins"].values()).issubset(ncins):
                            else:
                                LOGGER.warning(f"The following NetCDF files are missing: \n {opts['ncins']}")
                
def get_tempconversion_clargs(ncins, specs):
    """Get command line opts for the conversion of temperature"""
    
    srcs = []
                if "temp" in fields:
                    temp = "temp"
                    saln = "saln"
                elif "sst" in fields:
                    temp = "sst"
                    saln = "sss"
                    opts = {"ncins": {"nctemp": f"{temp}_{dim}_{src}_{term_postprod:{prec}}_{filter}.nc",
                                      "ncsaln": f"{saln}_{dim}_{src}_{term_postprod:{prec}}_{filter}.nc",
                                      },
                            "ncout": f"tempis_{dim}_{src}_{term_postprod:{prec}}_{filter}.nc"
                            }
                    if set(opts["ncins"].values()).issubset(ncins):
                        opts_list.append(opts)
                    else:
                        LOGGER.warning(f"The following NetCDF files do not exist: \n {opts['ncins']}")

def get_concat_clargs(ncins, rundate, postprod, 
                      vapp, vconf, specs):
    """Concatenation of daily files args"""
    
    srcs = []
                    cmdargs = {}
                    cmdargs["ncout"] = f"HYCOM{dim}_{src}_{term_postprod:{prec}}_{filter}.nc"
                    cmdargs["ncins"] = []
                        if (dim=="3D") and (field=="temp") and ("datashom" in postprod):
                            cmdargs["ncins"].append(f"tempis_3D_{src}_{term_postprod:{prec}}_{filter}.nc")
                        else:
                            cmdargs["ncins"].append(f"{field}_{dim}_{src}_{term_postprod:{prec}}_{filter}.nc")
                    if filter=="none":
                        field_type = "instantaneous"
                    else:
                        field_type = filter
                    date = pd.to_datetime(rundate)+pd.to_timedelta(int(term_postprod), "h")
                    if "forecast" in postprod:
                        if date==(pd.to_datetime(rundate)-pd.to_timedelta(1, "D")).strftime("%Y%m%d"):
                            forecast = "analysis"
                            forecast_range = f"1-days {forecast}"
                            forecast_type = f"{forecast}"
                        else:
                            forecast = "forecast"
                            forecast_range = f"5-days {forecast}"
                    elif "spnudge" in postprod:
                        day0 = convert_to_julian_day(np.datetime64(date))
                        day1 = convert_to_julian_day(np.datetime64(rundate)-np.timedelta64(7, 'D'))
                        if day0<day1:
                            forecast = "best analysis"
                        else:
                        forecast_range = f"7-days {forecast}"
                        forecast_type = f"{forecast}"
                    elif "hindcast" in postprod:
                        forecast = "best analysis"
                    if filter == "mean":
                        glob_attrs_title = f"{vapp} - {dim} parameters - daily mean data"
                    elif (filter == "godin") or (filter == "demerliac"):
                        glob_attrs_title = f"{vapp} - {dim} parameters - {filter} filtered data"
                    else:
                        glob_attrs_title = f"{vapp} - {dim} parameters"
                    cmdargs["global_attrs"] = {"title": glob_attrs_title,
                                               "source": f"Sloop-{vapp} {vconf}",
                                               "bulletin_date": pd.to_datetime(rundate).strftime("%Y-%m-%d %H:%M:%S"),
                                               "field_date": date.strftime("%Y%m%d") if ("datashom" in postprod) else pd.to_datetime(date).strftime("%Y-%m-%d %H:%M:%S"),
                                               "field_type": field_type,
                                               "forecast_range": forecast_range,
                                               "forecast_type": forecast_type,
                                               }
                    if (set(cmdargs["ncins"]).issubset(ncins)) and (len(cmdargs["ncins"])>0):
                        cmdargs_list.append([cmdargs])
                    else:
                        LOGGER.warning(f"The following input NetCDF files do not exist: {cmdargs['ncins']}") 

def concat_field_ncfiles(args):
    list_ds = []
    for ncfile in args["ncins"]:
        list_ds.append(ds)
    ds = xr.merge(list_ds)
    ds.encoding.update(list_ds[0].encoding)
    ds.attrs.update(list_ds[0].attrs)
    ds.attrs.update(args["global_attrs"])
    ds.to_netcdf(args["ncout"])
    
def get_extract_cmdargs(ncins, specs, postprod):
    """Area extractions args"""    
                    cmdargs["name"] = area
                    cmdargs["ncin"] = f"HYCOM{dim}_{src}_{term_postprod:{prec}}_{filter}.nc"
                    cmdargs["specs"] = specs[area]
                    cmdargs["postprod"] = postprod
                    if cmdargs["ncin"] in ncins:

def regrid_soap(ds):
    """Regrid dataset with constant latitude step"""
   
    dsin = ds.assign_coords({"X":ds.lon[0,:], 
                               "Y":ds.lat[:,0]}) 
    dsin = dsin.rename({"lon":"lon_2d", "lat":"lat_2d",
                      "X":"lon", "Y":"lat"})
    dlat = np.mean(dsin.lat.diff("lat"))
    lat = np.arange(dsin.lat.min(), dsin.lat.max(), dlat, 
                    type(ds.lat.values[0,0]))
    soap_grid = xr.DataArray(dims=["lon", "lat"],
                             name="grid",
                             coords=dict(lat=lat,
                                         lon=dsin.lon),
                             )
    regridder = xesmf.Regridder(dsin, soap_grid, 
                                reuse_weights=True,
                                method="bilinear")
    dsout = []
    for dain in dsin:    
        if not "time_bnds" in dain:
            daout = regridder(dsin[dain])
        else:
            daout = dsin[dain]
        daout.attrs.update(dsin[dain].attrs)
        daout.encoding.update(dsin[dain].encoding)
        dsout.append(daout)
    dsout = xr.merge(dsout)
    dsout = dsout.rename({"lon":"X", "lat":"Y"})
    lon, lat = np.meshgrid(dsout.X, dsout.Y)
    da_lon = xr.DataArray(lon, 
                          dims=["Y", "X"],
                          name="lon",
                          ) 
    da_lat = xr.DataArray(lat,
                          name="lat",
                          )
    dsout = xr.merge([dsout, da_lon, da_lat])
    for coord in dsout.coords:
        dsout[coord].attrs.update(ds[coord].attrs)
    dsout.encoding.update(ds.encoding)
    dsout.attrs.update(ds.attrs)
    return dsout

def extract_metazone(args):
    """Metazone extractions"""

    cwd = os.getcwd()
    worker_dir = args["name"]+"_"+args["ncin"][:-3]
    os.makedirs(worker_dir)
    shutil.copy(args["ncin"], worker_dir+"/")
    os.chdir(worker_dir)
    ncin_re = "HYCOM(?P<dim>.*)_(?P<src>.*)_(?P<term>.*)_(?P<filter>.*).nc"
    matches = re.search(ncin_re, args["ncin"], re.DOTALL)
    dsin = xr.open_dataset(args["ncin"])
    ds = dsin.assign_coords({"X":dsin.lon[0,:], "Y":dsin.lat[:,0]})
    if "soap" in args["postprod"]:
        ds = regrid_soap(ds)
        if "mean" in args["ncin"]:
            ds["time"] = ds["time"].dt.floor("1d")
        for var in ["u", "v"]:
            if("u" or "v") in ds:
                ds[var].attrs["units"] = "m s-1"
    else:
                ds[var].attrs["units"] = "m/s"
    encoding = {}
    for coord in ["depth", "time", "lon", "lat"]:
        if coord in ds:
            encoding.update({coord: {"_FillValue": None}})
            ds[coord].attrs.update(dsin[coord].attrs)
    ds["time"] = [(time - pd.to_datetime("1950-01-01")) / pd.to_timedelta("1d") for time in ds["time"].values]
    ds["time"].attrs.update(dsin["time"].attrs)
    ds["time"].attrs.update({"units": "days since 1950-01-01 00:00:00"})
    if "time_bnds" in ds:
        ds["time_bnds"] = ds["time_bnds"].astype(np.double)
        encoding.update({"time_bnds": {"_FillValue": None}})
    if not ((args["specs"]["bounds_lon"]["min"]==9999.9) or (args["specs"]["bounds_lon"]["max"]==9999.9)
            or (args["specs"]["bounds_lat"]["min"]==9999.9) or (args["specs"]["bounds_lat"]["max"]==9999.9)):
        dso = ds.sel(X=slice(args["specs"]["bounds_lon"]["min"], args["specs"]["bounds_lon"]["max"]),
                     Y=slice(args["specs"]["bounds_lat"]["min"], args["specs"]["bounds_lat"]["max"]))
    else:
        dso = ds
    dso = dso.drop(labels=['X', 'Y'])
    dso.attrs.update({
            "title": dso.attrs["title"]+" zone {}".format(args["name"]),
            "geospatial_lat_min": float(dso.lat.min().round(decimals=2)),
            "geospatial_lat_max": float(dso.lat.max().round(decimals=2)),
            "geospatial_lon_min": float(dso.lon.min().round(decimals=2)),
            "geospatial_lon_max": float(dso.lon.max().round(decimals=2)),
        })
    ncout = f"HYCOM{matches['dim']}_{matches['term']}_{matches['filter']}_{args['name']}.nc"
    dso.to_netcdf(ncout, encoding=encoding)
    shutil.move(ncout, "../")
    os.chdir(cwd)

def get_vtx_geometries(cmdargs):
    """Get the post-production geometries."""

    geometries = dict()
    for cmdarg in cmdargs:
        name = cmdarg[0]["name"]
        geometries[name] = dict(
            kind = "lonlat",
            info = "Shom Hycom3d postproduction area",
            area = name,
        )
    return geometries

"""
Rivers processing
"""
import tarfile
import re

from ...interp import interp_time
from .__init__ import (HYCOM3D_RIVERS_CFG, 
                       HYCOM3D_RIVERS_INI,
                       HYCOM3D_RIVERS_FLOW_CLIM,
                       HYCOM3D_RIVERS_HYDRO_CLIM)
from ...io import read_inicfg_files

class Rivers(object):
    """ Rivers prepocessings """

            self, cfgfile=HYCOM3D_RIVERS_CFG.format(rank=0)):
        self.rivers = read_inicfg_files(
            fcfg=cfgfile, fini=HYCOM3D_RIVERS_INI
        )
        self.platforms = dict()
        for river in self.rivers.keys():
            self.rivers[river].update({"dataset": None})
            if not self.rivers[river]["clim"]:
                LOGGER.info(f"{river}: real time data used")
                for platform in self.rivers[river]["platform"]["Id"]:
                    self.platforms[platform] = {"path": [], "dataset": None}
            
    def platforms_selection(self, tarname):
        nc_pattern = r"GL_TS_RF_(?P<platform>.*)_(?P<date>.*).nc"
        tmp_tar = tarfile.open(tarname, mode="r:gz")
        for member in tmp_tar.getnames():
            if any(platform in member for platform in self.platforms.keys()):
                matches = re.search(nc_pattern, member, re.DOTALL)
                tmp_tar.extract(member)
                LOGGER.info(f"{matches['platform']}: {member}")
                self.platforms[matches["platform"]]["path"].append(member)
        tmp_tar.close()
        return self.platforms

    def platforms_nc2dataset_preprocess(self, ds):
        ds["LONGITUDE"] = ds.LONGITUDE.swap_dims({"LONGITUDE": "TIME"})
        ds["LATITUDE"] = ds.LATITUDE.swap_dims({"LATITUDE": "TIME"})
        ds = ds['RVFL'].squeeze(dim="DEPTH", drop=True)
        ds = ds.rename({"TIME": "time", "LONGITUDE": "lon", "LATITUDE": "lat"})
        ds = ds.reset_coords(["lon", "lat"])
        return ds

    def platforms_nc2datasets(self, tarname):
        LOGGER.info("NetCDF to dataset for: ")
        self.platforms = self.platforms_selection(tarname)
        for platform in self.platforms.keys():
            LOGGER.info(platform)
            dataset = xr.open_mfdataset(
                self.platforms[platform]["path"][:],
                preprocess=self.platforms_nc2dataset_preprocess,
                concat_dim=["time"],
                combine="by_coords",
            )
            self.platforms[platform]["dataset"] = dataset
        return self.rivers, self.platforms

    def write_flowrate(self, tarname, time, nc_out):
        self.rivers, self.platforms = self.platforms_nc2datasets(tarname)
            LOGGER.info(river)
                LOGGER.info("Real time data")
                    LOGGER.info(f"{platform} time interpolation")
                    ds = interp_time(
                        self.platforms[platform]["dataset"].load(), time, 
                        keep_flag=True, persistency="mean_25h"
                    )
                    ds.update(
                        {
                            "RVFL": ds["RVFL"]
                            * self.rivers[river]["platform"][platform]["debcoef"]
                        }
                    )
                        {
                            "flag": ds["flag"]
                            / len(self.rivers[river]["platform"]["Id"])
                        }
                    )
                    if self.rivers[river]["dataset"]:
                        for var in ["RVFL", "flag"]:
                            self.rivers[river]["dataset"][var] += ds[var]
                    else:
                        self.rivers[river]["dataset"] = ds
            else:
                LOGGER.info("Climatology data")
                LOGGER.info(f"Read {HYCOM3D_RIVERS_FLOW_CLIM} ascii file")
                f = open(HYCOM3D_RIVERS_FLOW_CLIM, mode="r", encoding="ascii")
                while True:
                    line = f.readline()
                    if river in line:
                        LOGGER.debug(line)
                        break
                rvfl = np.tile(np.float_(f.readline().split()), 3)
                f.close()
                LOGGER.info("Get climatology time")
                time_fr = [pd.to_datetime(t) for t in self.get_time_clim(time, len(rvfl))]
                LOGGER.debug(time_fr)
                da = xr.DataArray(rvfl,
                          name="RVFL",
                          coords={
                              "time":time_fr},
                          dims='time').to_dataset()
                LOGGER.info("Time interpolation")
                ds = interp_time(
                        da, time, keep_flag=True
                        )
                self.rivers[river]["dataset"] = ds
            LOGGER.info("{river}: save in NetCDF file")
            self.rivers[river]["dataset"][["RVFL", "flag"]].to_netcdf(
                nc_out.format(**locals())
            )

    def write_tempsaln(self, nc_in, nc_out):
        const = 2.0 * np.pi / 365.0
            LOGGER.info(f"{river}:")
            ds_river = xr.open_dataset(nc_in.format(**locals()))
                for var in ["temperature", "salinity"]:
                    tmax = xr.DataArray(
                        pd.to_datetime(
                            str(ds_river.time.values[0])[:4]
                            + self.rivers[river][var]["datemax"]
                        ),
                        name="tmax",
                    )
                    DeltaT = (ds_river["time"]-tmax)/pd.to_timedelta(1, "D")
                    DeltaT = DeltaT*const
                    VAR = self.rivers[river][var]["avg"]+self.rivers[river][var][
                        "amp"
                    ] * np.cos(DeltaT)
                    xa = xr.DataArray(
                        VAR, coords=[ds_river.time], dims=["time"], name=var
                    )
                    ds_river = xr.merge([ds_river, xa])
            else:
                LOGGER.info(f"Open {HYCOM3D_RIVERS_HYDRO_CLIM} ascii file") 
                f = open(HYCOM3D_RIVERS_HYDRO_CLIM, mode="r", encoding="ascii")
                        break
                temp = np.tile(np.float_(line.split()[2:]), 3) 
                saln = np.tile(np.float_(f.readline().split()[2:]), 3)
                f.close()
                time_fr = [pd.to_datetime(t) for t in self.get_time_clim(ds_river.time.values, len(temp))]
                da_temp = xr.DataArray(temp,
                          name="temperature",
                          coords={
                da_saln = xr.DataArray(saln,
                          name="salinity",
                          coords={
                LOGGER.debug("Time interpolation")
                        xr.merge([da_temp, da_saln]), ds_river.time.values, keep_flag=True
                        )
                ds_river = xr.merge([ds_river, ds])
            LOGGER.info(f"{river} : save in NetCDF file")
            ds_river.to_netcdf(nc_out.format(**locals()))

    def write_rfiles(self, nc_in="{river}.flx.ts.nc", r_out="{river}.r"):

            "Debits horaires du fleuve {river} \n"
            "origine : Donnees CMEMS \n\n"
            '"Date"                  '
            '"Debit"  "temp."  "sal."  "flag" \n'
            '"jour depuis 1/1/1950"  "m3/s"   "DegC"   '
            '"1e-3"   "1=data, 2=interp, 3=persistance"  \n'
            "--------------------------------------------"
            "--------------------------------------- \n"
        )
            "  {julianday}    {debit}    {temp}    {saln}    {flag}\n"
        )
        valformat = "{:.10f}"

            with open(r_out.format(**locals()), "w") as f:
                f.writelines(header_pattern.format(**locals()))
                for itime in range(len(ds_river.time)):
                    dsi = ds_river.isel({"time": itime})
                    julianday = valformat.format(
                        convert_to_julian_day(dsi.time.values)
                    )
                    debit = valformat.format(dsi.RVFL.values)
                    temp = valformat.format(dsi.temperature.values)
                    saln = valformat.format(dsi.salinity.values)
                    flag = int(dsi.flag.values)
                    f.write(line_pattern.format(**locals()))

    def write_ncfiles(self, nc_in="{river}.flx.ts.nc", nc_out="{river}.nc"):

    def get_time_clim(self, time, ntime):
        time_fr = ["{}0115".format(
                    (pd.to_datetime(time[0])-pd.DateOffset(years=1)).strftime("%Y"))]
        itime = 0
        while itime<(ntime-1):
            t = "{}15".format(
                (pd.to_datetime(time_fr[itime])+pd.DateOffset(months=1)).strftime("%Y%m"))
            time_fr.append(t)
            itime += 1
        return time_fr

"""
Spectral Nudging preprocessings
"""

from ...log import LoggingLevelContext 

from xoa.coords import get_time
from .io import read_blkdat_input, hycom2nc, read_regional_grid
from . import HYCOM3D_GRID_AFILE, HYCOM3D_BLKDAT_INPUT_FILE
from sloop.times import convert_to_julian_day
from sloop.interp import interp_time

def get_demerliac_date(ncin):
    ds = xr.open_dataset(ncin)
    return get_time(ds)

def demerliac_sel(ncins, outfile_pattern="{field}_timesel.nc", dt=None):
    LOGGER.info("Time selection for demerliac filtering") 
    LOGGER.info("time delta (dt) in day is mandatory")
    if dt==None:
        dt = pd.to_timedelta(1.46, 'D')
    else:
        dt = pd.to_timedelta(dt, 'D')
    LOGGER.debug("time delta: {dt}")
    ncouts = {}
    for ncin in ncins:
        LOGGER.info(f"Open {ncin}")
        ds = xr.open_dataset(ncin)
        ds = ds.sel(time=slice(ds.time.mean()-dt, ds.time.mean()+dt, 1))
        field = ncin.split('.')[0]
        ncout = outfile_pattern.format(**locals())
        ncouts[field] = ncout
        LOGGER.info(f"Store in {ncout}")
        ds.to_netcdf(ncout)
    date_demerliac = convert_to_julian_day(ds.time.mean().values)
    return date_demerliac, ncouts

def spectral_sel(ncins, time, outfile_pattern="{field}-mercator.nc", rank="0"):
    LOGGER.info("Time selection for spectral filtering on nesting input files")
    LOGGER.info("Read blkdat.input")
    dsb = read_blkdat_input(HYCOM3D_BLKDAT_INPUT_FILE.format(rank=rank))
    nx = int(dsb.idm)
    ny = int(dsb.jdm)
    nz = int(dsb.kdm)    
    ncouts = {}
        if ncin.endswith('.b'):
            LOGGER.info(f"Processing on {ncin}")
            nc = hycom2nc(nx, ny, nz,
                          HYCOM3D_GRID_AFILE.format(rank=rank),
                          ncin)
            ds = xr.open_dataset(nc)
            ds = interp_time(ds, [str(time.values[0])], preserve_encoding=True)
            field = (ncin.split('.')[0]).split('-')[0]
            ncout = outfile_pattern.format(**locals())
            ncouts[field] = {"mercator": ncout}
            ds = ds.fillna(0)
            LOGGER.info(f"Store results in {ncout}")
            ds.to_netcdf(ncout)
    return ncouts

def ncdiff(varname, nc1, nc2, outfile_pattern='{varname}-spn.nc', rank='0'):
    LOGGER.info("Compute difference between filtered mercator and filtered hycom fields")
    ds = read_regional_grid(HYCOM3D_GRID_AFILE.format(rank=rank), grid_loc='p')
    ds = ds.rename({"idm":"x", "jdm":"y"})
    ds1, ds2 = xr.open_dataset(nc1), xr.open_dataset(nc2)
    if varname=="t":
        longname = "temp"
    elif varname=='s':
        longname = "saln"
    else:
        longname = varname
    field = ds2[longname].load().values-ds1[longname].load().values
    field[field==0] = np.nan
    da = xr.DataArray(field, name=varname,
                      dims=('time', 'lev', 'y', 'x')
                      )
    da.encoding["_FillValue"]= np.nan
    ds = ds.assign_coords(time=ds2.time)
    ncout = outfile_pattern.format(**locals())
    LOGGER.info(f"Store results in {ncout}")

__version__ = "0.1.0"

import os
import re

#: Directory of the python library
SLOOP_PYTHON_DIR = os.path.dirname(__file__)

#: Main sloop dir
SLOOP_DIR = os.path.realpath(os.path.join(SLOOP_PYTHON_DIR, '..'))

#: Root directory of the fortran sources
SLOOP_SRC_DIR = os.path.join(SLOOP_DIR, 'src')

#: Root directory of the templates
SLOOP_TPL_DIR = os.path.join(SLOOP_DIR, 'templates')

_REQUIREMENTS_FILE = os.path.join(
    os.path.dirname(__file__), "..", "requirements.txt"
)

WORKFLOW_INI = os.path.join(
    SLOOP_PYTHON_DIR, "workflow.ini"
)

class SloopError(Exception):
    pass

def _parse_requirements_(reqfile):
    re_match_specs_match = re.compile(r"^(\w+)(\W+.+)?$").match
    reqs = {}
    with open(reqfile) as f:
            line = line.strip().strip("\n")
            if line and not line.startswith("#"):
                m = re_match_specs_match(line)
                if m:
                    reqs[m.group(1)] = m.group(2)
    return reqs

def show_versions():
    """Print the versions of sloop and of some dependencies
    """
    from pkg_resources import get_distribution, DistributionNotFound
    print("- sloop:", __version__)
    for package in _parse_requirements_(_REQUIREMENTS_FILE):
        try:
            version = get_distribution(package).version
        except DistributionNotFound:
            version = "NOT INSTALLED or UKNOWN"
        print(f"- {package}: {version}")

def show_sloop_paths():
    """Print some internal sloop paths
    """
    for label, path in [("main dir", SLOOP_DIR),
                        ("python lib dir", SLOOP_PYTHON_DIR),
                        ("model src dir", SLOOP_SRC_DIR)]:
        print("-", label+":", path)

def show_info(opt_specs=True, exp_dir=None):
    """Print sloop related info
    """
    print("# VERSIONS")
    show_versions()
    print("\n# SLOOP PATHS")
    show_sloop_paths()
    if exp_dir:
        from .env import show_exp_info
        print("\n# EXPERIMENT DIRECTORY")
        show_exp_info(exp_dir)

"""
CF utilities
"""
import re
import warnings

CF_SPECS = {
    'lon': {
        'name': ['lon', 'longitude', 'nav_lon','LON'],
        'units': ['degrees_east', 'degree_east', 'degree_E', 'degrees_E',
                  'degreeE', 'degreesE'],
        'standard_name': ['longitude'],
        'axis': 'X',
        },
    'lat': {
        'name': ['lat', 'latitude', 'nav_lat','LAT'],
        'units': ['degrees_north', 'degree_north', 'degree_N', 'degrees_N',
                  'degreeN', 'degreesN'],
        'standard_name': ['latitude'],
        'axis': 'Y'
        },
    'depth': {
        'name': ['depth', 'deptht'],
        'units': ['m'],
        'standard_name': ['depth'],
        'axis': 'Z'
        },
    'temp': {
        'name': ['temp', 'votemper'],
        'units': ['degrees_celsius', 'degrees_C'],
        'standard_name': ['sea_water_temperature',
                          "sea_water_potential_temperature"],
        },
    'saln': {
        'name': ['saln', 'votemper'],
        'units': ['PSU'],
        'standard_name': ['sea_water_salinity'],
        },
    'ssh': {
        'name': ['ssh', 'xe', 'sossheig'],
        'standard_name': ['sea_surface_height_above_sea_level',
                          'sea_surface_height_above_geoid'],
        },
    }

RE_SN_AT_LOC = re.compile(r'^(\w+)(?:_at_(\w)_location)$').match

def get_dim_name_from_axis(ds, axis):
    """
    Get a dimension name from its axis attribute in a xarray.Dataset

    ds : xarray.Dataset
    axis : str
        X, Y, Z or T

    Returns
    -------
    str

    """
    axis = axis.upper()
    assert axis in "XYZT", "Invalid value for axis attribute"
    for name, coord in ds.coords.items():
        if hasattr(coord, 'axis') and coord.axis == axis:
            return name
    raise AttributeError('Dimension from axis not found: '+axis)

def get_x_dim_name(ds):
    '''
    Get the X dimension name in a `xarray.Dataset` or a `xarray.DataArray`

    ds: xarray.DataArray or xarray.Dataset

    Returns
    -------
    str: Name of X dimension
    '''
    return get_dim_name_from_axis(ds, 'X')

def get_y_dim_name(ds):
    '''
    Get the Y dimension name in a `xarray.Dataset` or a `xarray.DataArray`

    Returns
    -------
    str: Name of Y dimension
    '''
    return get_dim_name_from_axis(ds, 'Y')

def get_t_dim_name(ds):
    '''
    Get the time dimension name in a `xarray.Dataset` or a `xarray.DataArray`

    Returns
    -------
    str: Name of time dimension
    '''
    return get_dim_name_from_axis(ds, 'T')

def get_name_from_specs(ds, specs, check=None, first=True):
    '''
    Get the name of a `xarray.DataArray` according to search specifications

    specs: dict
    check: str, list of str
        Attributes to check
    first: bool

    Returns
    -------
    str: Name of the xarray.DataArray
    '''
    # Targets
    ds_names = []
    if hasattr(ds, 'coords'):
        ds_names.extend(list(ds.coords))
    if hasattr(ds, 'variables'):
        ds_names.extend(list(ds.variables))
    if isinstance(check, str):
        check = [check]

    # Search in names first
    if 'name' in specs and (check is None or 'name' in check):
        for specs_name in specs['name']:
            if specs_name in ds_names:
                return specs_name

    # Search in units attributes
    if ('units' in specs and 'axis' in specs and
            (check is None or 'units' in check)):
        for ds_name in ds_names:
            if hasattr(ds[ds_name], 'units'):
                for specs_units in specs['units']:
                    if ds[ds_name].units == specs_units:
                        return ds_name

    # Search in standard_name attributes
    if ('standard_name' in specs and
            (check is None or 'standard_name' in check)):
            if hasattr(ds[ds_name], 'standard_name'):
                for specs_sn in specs['standard_name']:
                    if ds[ds_name].standard_name == specs_sn:
                    m = RE_SN_AT_LOC(ds[ds_name].standard_name)
                    if m and m.group(1) == specs_sn:

    # Search in axis attributes
    if 'axis' in specs and (check is None or 'axis' in check):
            if (hasattr(ds[ds_name], 'axis') and
                    ds[ds_name].axis.lower() == specs['axis'].lower()):
                return ds_name

def get_lon_coord_name(ds):
    '''
    Get the longitude name in a `xarray.Dataset` or a `xarray.DataArray`

    Returns
    -------
    '''
    lon_name = get_name_from_specs(ds, CF_SPECS['lon'])
    if lon_name is not None:
        return lon_name
    warnings.warn('longitude not found in dataset')

def get_lat_coord_name(ds):
    '''
    Get the latitude name in a `xarray.Dataset` or a `xarray.DataArray`

    Returns
    -------
    '''
    lat_name = get_name_from_specs(ds, CF_SPECS['lat'])
    if lat_name is not None:
        return lat_name
    warnings.warn('latitude not found in dataset')

def get_time_coord_name(ds):
    '''
    Get the time coordinate name from its dtype

    Returns
    -------
    '''
    for coord_name in ds.coords:
        if ds.coords[coord_name].dtype.char == 'M':
            return coord_name
    warnings.warn('time not found in dataset')

def rename_from_standard_name(ds, standard_name_to_name, mapping=None):
    """
    Rename coords and variables according to their standard_name attribute

    standard_name_to_name : dict
        Key are standard_names and values are names.

    Returns
    -------

    """
    if isinstance(ds, xr.DataArray):
        keys = ds.coords.keys()
    else:
        keys = ds.variables.keys()
    keys = list(keys)

    for standard_name, name in standard_name_to_name.items():

        names = []
        matcher = re.compile(r'^{}(_.+)?$'.format(standard_name)).match
        for vname in keys:
            if (hasattr(ds[vname], 'standard_name') and
                    matcher(ds[vname].standard_name)):
                names.append(vname)
        if not names:
            continue
        assert len(names) <= 1, ('No more than one variable must have this'
                                 ' standard_name: '+standard_name)
        vname = names[0]
        if vname != name:
            ds = ds.rename({vname: name})
            keys.remove(vname)
            keys.append(name)
            if isinstance(mapping, dict):
                mapping[name] = vname

    return ds

# with warnings.catch_warnings():
#     warnings.simplefilter(
#         "ignore", xr.core.extensions.AccessorRegistrationWarning)

#     @xr.register_dataset_accessor('cf')
#     @xr.register_dataarray_accessor('cf')
#     class CFAccessor(object):

#         def __init__(self, ds):
#             self.ds = ds

#         def __getitem__(self, cf_name):
#             return self.ds[get_name_from_specs(self.ds, CF_SPECS[cf_name])]

#         def get_lon(self):
#             return self['lon']

#         lon = property(fget=get_lon, doc='Longitude')

#         def get_lat(self):

#         lat = property(fget=get_lat, doc='Latitude')

#         def get_depth(self):
#             return self['depth']

#         depth = property(fget=get_depth, doc='Depth')

#         def get_time(self):
#             return self.ds[get_time_coord_name(self.ds)]

#         time = property(fget=get_time, doc='Time')

import sys
import fnmatch
import socket
import os
import importlib
import pathlib
import re
import subprocess
from .__init__ import SloopError
from .log import LoggingLevelContext
from configparser import ConfigParser

#: Database for guessing HPC target from hostname
HOSTNAME_TO_HPC_TARGET = {
    "datarmor": ["*.ice.ifremer.fr"],
    "aquarium": ["aqlogin*"],
    "meteofrance": ["*hpc.meteo.fr"],
}

RE_CYCLE_MATCH = re.compile(
    r"^(?P<scheme>[^:]*):(?P<cfgname>[^@]*)(?:@(?P<location>.+))?$"
).match

def load_config_from_dict(conf, sec):
    """Load environment variables from dict"""
    env_dict = conf[sec]
    try:
        for var in env_dict.keys():
            if var.startswith("env_"):
                varname = var.split("env_")[-1]
                if isinstance(env_dict[var], list):
                    value = " ".join(comp[var])
                else:
                    value = env_dict[var]
                os.environ[varname.upper()] = value
    except SloopEnvError as e:
        logger.error(e.msg)

def get_section_from_config(conf, section):
    """Load section in config.ini file"""
    config = ConfigParser(default_section=None)
    config.read(conf)
    try:
        sec = config[section]
    except Exception:
        raise SloopEnvError(f"Section {section} not in {conf} file")
    return dict(sec.items())

def get_field_from_config(conf, section, key):
    """Load config from a section in config.ini file"""
    config = ConfigParser()
    try:
        field = config[section][key]
        raise SloopEnvError(f"Section {section} or {key} not in {conf} file")

    return field

def load_env_config(conf, sec):
    # start filling env variables
    if sec in config.sections():
        comp = config[sec]
        for var in comp.keys():
                if isinstance(comp[var], list):
                else:
                    value = comp[var]
    else:
        raise SloopEnvError(f"Section {sec} not in config file")

def stripout_conda_env(env_dict=None):
    """Remove conda env from PATH, LD_LIBRARY_PATH and C_INCLUDE_PATH"""
    if env_dict is None:
        env_dict = os.environ
    for name, subdir in (
        ("PATH", "bin"),
        ("LD_LIBRARY_PATH", "lib"),
        ("LIBRARY_PATH", "lib"),
        ("INCLUDE_PATH", "include"),
        ("C_INCLUDE_PATH", "include"),
    ):
        if name in env_dict:
            paths = env_dict[name].split(":")
            path = os.path.join(sys.prefix, subdir)
            if path in paths:
                paths.remove(path)
            env_dict[name] = ":".join(paths)

class ContextEnv(object):
    """Context for changing env vars"""

    def __init__(self, stripout_conda=False, **envs):
        self.stripout_conda = stripout_conda
        self.envs = envs

    def __enter__(self):
        self.old_env = os.environ.copy()
        if self.stripout_conda:
            stripout_conda_env()
        os.environ.update(**self.envs)
        return self

    def __exit__(self, exc_type, exc_value, traceback):
        os.environ = self.old_env

def guess_hpc_target():
    """Get the current HPC target name

    It first look for the :envvar:`HPC_TARGET` env variable, then use the
    :data:`HOSTNAME_TO_HPC_TARGET` variable to guess it from hostname.
    If not a well known HPC target, "generic" is returned.
    The :envvar:`HPC_TARGET` is then set to this value.
    """
    if "HPC_TARGET" in os.environ:
        return os.environ["HPC_TARGET"]
    hostname = socket.getfqdn()
    for hpc_target, patterns in list(HOSTNAME_TO_HPC_TARGET.items()):
        for pattern in patterns:
            if fnmatch.fnmatch(hostname, pattern):
                os.environ["HPC_TARGET"] = hpc_target
                return hpc_target
    os.environ["HPC_TARGET"] = "generic"
    return "generic"

def is_os_cmd_avail(cmd):
    """Check if a command is available

    cmd: str
        Command can be executed in a terminal

    Return
    ------
    Boolean
    """
    return subprocess.run(["which", cmd], stdout=subprocess.DEVNULL).returncode == 0

class SloopEnvError(SloopError):
    def __init__(self, msg):
        self.msg = msg

def check_env_vars(env_vars):
    """Check that the given list of env variables are set

    env_vars: str, list
        List of environment variables

    Return
    ------
    SloopEnvError: if a variable is not set
    """
    if isinstance(env_vars, str):
        env_vars = [env_vars]
    for env_var in env_vars:
        if env_var not in os.environ:
            raise SloopEnvError(
                f"Environment variable not set: {env_var}. "
                + "The following variales must be set: "
                + " ".join(env_vars)
            )

def check_vortex_dir(vortex_dir):
    """Check that the vortex directory is valid

    Return
    ------
    SloopEnvError
    """
    for relpaths in [
        ("",),
        ("bin", "mkjob.py"),
        ("src", "common"),
        ("site", "footprints"),
    ]:
        if not os.path.exists(os.path.join(vortex_dir, *relpaths)):
            raise SloopEnvError(f"Invalid vortex dir: {vortex_dir}")

def check_vapp(vapp):
    """Check that vapp is valid

    vapp must be a module importable from :mod:`sloop.models`.

    Return
    ------
    """
    vapp = str(vapp)
    try:
        with LoggingLevelContext(logging.WARNING):
            vapp_module = importlib.import_module("." + vapp, "sloop.models")
    except ImportError:
        raise SloopEnvError("Invalid vapp: " + vapp)
    return vapp_module

def check_exp_dir(exp_dir, errors="raise"):
    """Check that an experiment directory is valid

    This checks parents and children.

    Return
    ------
    dict

    See also
    --------
    check_vapp
    """

    assert errors in ("raise", "ignore")
    assert "MTOOLDIR" in os.environ, "MTOOLDIR has to be set"

    # The exp dir itself
    specs = {}
    if not os.path.exists(exp_dir):
        if errors == "raise":
            raise SloopEnvError(f"Experiment directory not found: {exp_dir}")
        else:
            return specs
    else:
        exp_dir = pathlib.Path(os.path.abspath(exp_dir))
    specs["exp_dir"] = exp_dir

    # xpid/vapp/vconf
    vconf = exp_dir.name
    vapp_dir = exp_dir.parent
    vapp = vapp_dir.name
    vapp_module = check_vapp(vapp)
    xpid = vapp_dir.parent.name
    user = os.path.basename(os.environ["HOME"])  # getpass.getuser()
    location = user if "@" not in xpid else xpid.split("@")[1]
    specs.update(
        vapp=vapp,
        vconf=vconf,
        xpid=xpid,
        vapp_module=vapp_module,
        location=location,
        exp_dir=str(exp_dir),
    )

    # Links
    for tool in ["vortex", "sloop"]:
        if not (exp_dir / tool).is_symlink():
            if errors == "raise":
                raise SloopEnvError(f"Link to tool not found: {tool}")
            else:
                specs[tool] = "no local link found!"
        else:
            specs[tool] = os.path.abspath(str(exp_dir / tool))

    # Conf
    # - dir
    conf_dir = exp_dir / "conf"
    if not conf_dir.exists():
            raise SloopEnvError(f"Conf dir not found: {conf_dir}")
        else:
            specs["conf_dir"] = "not found!"
    else:
        specs["conf_dir"] = str(conf_dir)
    # - ini file
    ini_file = conf_dir / f"{vapp}_{vconf}.ini"
    if not ini_file.exists():
            raise SloopEnvError(f"Ini file not found: {ini_file}")
        else:
            specs["ini_file"] = "not found !"
    else:

        specs["ini_file"] = str(ini_file)

        # uenv file
        with ContextEnv(VORTEX_QUIET="1"):
            from vortex.util.config import GenericConfigParser

            cfg = GenericConfigParser(specs["ini_file"])
            if cfg.has_option(None, "cycle"):
                m = RE_CYCLE_MATCH(cfg.defaults()["cycle"])
                if m.group("scheme").startswith("u"):
                    from vortex import ticket

                    t = ticket()
                    location = m.group("location") or location
                    specs["uget_hack_file"] = os.path.join(
                        t.glove.configrc,
                        "hack",
                        "uget",
                        location,  # xpid location or user?
                        "env",
                        m.group("cfgname"),
                    )

    # Tasks
    # - dir
    tasks_dir = exp_dir / "tasks"
    if not tasks_dir.exists():
            raise SloopEnvError(f"Tasks dir not found: {tasks_dir}")
        else:
            specs["tasks_dir"] = "not found!"
    else:
        specs["tasks_dir"] = str(tasks_dir)
    # - init file
    init_task = tasks_dir / "__init__.py"
    if not init_task.exists():
            raise SloopEnvError(f"Tasks init file not found: {init_task}")
        else:
            specs["tasks_init_file"] = "not found!"
    else:
        specs["tasks_init_file"] = str(init_task)
    # - task files
    taskfiles = list(tasks_dir.glob("*.py"))
    specs["task_files"] = [str(t) for t in taskfiles if not t.name == "__init__.py"]

    return specs

def show_exp_info(exp_dir):
    """Print some info about an experiment

    exp_dir: str
        A valid experiment directory

    See also
    --------
    check_exp_dir
    """
    specs = check_exp_dir(exp_dir, errors="ignore")
    for label, content in specs.items():
        label = label.replace("_", " ")
        if isinstance(content, list):
            print("- " + label + ":")
            for value in content:
                print("  - " + value)
        else:
            if hasattr(content, "__file__"):
                content = content.__file__
            print("-", label + ":", content)

"""
Filters
"""

from .cf import get_x_dim_name, get_y_dim_name

def erode_coast(da, kernel=3, niter=1):
    """
    Erode the coast using a local average

    da : xarray.DataArray
        Data array to erode.
    kernel : int, dict, tuple, optional
        Size of kernel in X and Y directions.

    Returns
    -------
    xarray.DataArray.

    """
    # Masked data?
    mask = da.isnull().data
    if not mask.any():
        return da
    valid = ~mask

    # Local import to always import scipy
    from scipy.signal import convolve

    # kernel
    if not isinstance(kernel, dict):
        x_dim = get_x_dim_name(da)
        y_dim = get_y_dim_name(da)
        if isinstance(kernel, int):
            kernel = (kernel, kernel)
        kernel = {x_dim: kernel[0], y_dim: kernel[1]}
    kk = ()
    for dim in da.dims:
        kk += kernel.get(dim, 1),
    kernel = np.ones(kk, dtype='l')

    # Convolve
    vnum = da.fillna(0.)
    for i in range(niter):
        vnum = convolve(vnum, kernel, mode='same')
        vdenom = convolve((~mask).astype('l'), kernel, mode='same')
        mask = vdenom == 0
        vdenom[mask] = 1
        vnum /= vdenom
        del vdenom

    # Output
    return da.where(valid, np.where(mask, np.nan, vnum))

def erode_mask_vec(ds, kernel, param=None, until=5):
    """
    Erode the mask of a vectorized parameter

    param : str, optional
        Type of vectorized parameter
    
    Returns
    -------

    """
    LOGGER.info(f"Erode mask on {param} vectorized variables")
    if param == 'wind':
        tau = (ds.ewss**2 + ds.nsss**2)**0.5
        si10 = ds.si10**2
        si10 = si10.where(si10 > 1e-6, 1e-6)
        rhocd = tau / si10
        rhocd1 = rhocd.where(ds.mask == 1)
        rhocd2 = erode_mask(rhocd1, kernel=kernel, until=until)
        maskd = rhocd2.where(xr.ufuncs.isfinite(rhocd2), np.nan)
        maskd = maskd.where(xr.ufuncs.isnan(maskd), 1)
        w = rhocd2 / rhocd
        w = w.where(w < 1, 1)
        for var in ['ewss', 'nsss']:
            ds[var] = ds[var] * w * maskd

    return ds

"""
Grid utilities
"""

from xoa.cf import get_cf_specs

def geosel(ds, lon, lat, pad=0):
    """Geographic selection"""
    # Cf specs

    # Input grid coordinates
    glon = cfs.coords.get(ds, 'lon')
    glat = cfs.coords.get(ds, 'lat')

    # Selection specs
    lon_min = float(min(lon)) - pad
    lon_max = float(max(lon)) + pad
    lat_min = float(min(lat)) - pad
    lat_max = float(max(lat)) + pad

    # Rectangular selection
    if glon.ndim == 1:
        return ds.sel({glon.name: slice(lon_min, lon_max),
                       glat.name: slice(lat_min, lat_max)})

    # Curvilinear selection
    return ds.where((glon >= lon_min) & (glon <= lon_max) &
                    (glat >= lat_min) & (glat <= lat_max), drop=True)

class GeoSelector(object):
    """Geographic selector"""

    def __init__(self, lon, lat, pad=0):
        self.lon = lon
        self.lat = lat
        self.pad = pad

    def select(self, ds):
        """Apply slection with :func:`geosel`"""
        return geosel(ds, self.lon, self.lat, self.pad)

    __call__ = select

"""
Interpolation stuff
"""
import os
import re
from enum import IntEnum
import xoa.cf
from xoa.misc import XEnumMeta
from .cf import (rename_from_standard_name,
                 get_lon_coord_name, get_lat_coord_name)
from .filters import erode_coast
from .io import xr_open_concat
from .times import period_vortex_to_pandas

class InterpFlag(IntEnum, metaclass=XEnumMeta):
    PICKED = 1
    LINEAR = 2
    PERSISTENT = 3

def interp_time(ds, dates, preserve_encoding=False, keep_flag=False, diurnal_vars=None, persistency=None):
    """Interpolate a dataset in time

    A single interpolation per time step, with persistency.

        Dataset to interpolate.
    dates : single or list of dates
        Interpolation dates.

    Returns
    -------

    """
    # Time dim
    cfs = xoa.cf.get_cf_specs()
    time_name = cfs.search_coord(ds, cf_name="time", get="cf_name")
    time = ds.indexes[time_name]
    begindate_nc = time[0]
    enddate_nc = time[-1]
    dates = pd.to_datetime(dates)

    # For persistency
    ds_first = ds.isel({time_name: slice(0, 1)})
    ds_last = ds.isel({time_name: slice(-1, None)})
    flag = xr.DataArray(np.zeros(len(dates), dtype="i"),
                        coords=[dates], dims=[time_name], name='flag')

    # Loop on target times
    dss = []
    for date in dates:
        # Persistency
        if date < begindate_nc or date > enddate_nc:
            # New dataset
            this_ds = (ds_first if date < begindate_nc else ds_last)
        
            # New time
            this_ds = this_ds.assign_coords({time_name: [date]})
            this_ds[time_name].attrs.update(ds[time_name].attrs)

            if persistency:
                m = re.search(r"mean_(?P<term>.*)(?P<units>.)", persistency, re.DOTALL)
                if m:
                    ds_tr = ds.copy()
                    timedelta = pd.to_timedelta(int(m["term"]), m["units"])
                    if date > enddate_nc:
                        date_target = enddate_nc-timedelta
                        ds_tr = ds_tr.where(ds_tr[time_name]>=date_target).mean(dim=time_name, keepdims=True)
                    elif date < begindate_nc:
                        date_target = begindate_nc+timedelta
                        ds_tr = ds_tr.where(ds_tr[time_name]<=date_target).mean(dim=time_name, keepdims=True)
                    for vname in ds_tr:
                        this_ds[vname] = this_ds[vname].copy(data=ds_tr[vname].data)

            if diurnal_vars:
                if date < begindate_nc:
                    time_target = date+pd.Timedelta("24h")
                elif date > enddate_nc:
                    time_target = date-pd.Timedelta("24h")
                for dss_target in dss:
                    if dss_target[time_name].values[0]==time_target:
                        for vname in this_ds:
                            if vname in diurnal_vars:
                                this_ds[vname] = this_ds[vname].copy(data=dss_target[vname].data)
           
            flag.loc[date] = InterpFlag.PERSISTENT

        # Interpolation
        else:

            # Pick this date directly if available
            if date in time:
                this_ds = ds.sel({time_name: [date]})
                flag.loc[date] = InterpFlag.PICKED

            # Interpolate because not available
            else:
                # Load a two-time steps dataset
                it = ds.indexes[time_name].get_loc(date, method='ffill')
                tmp_ds = ds.isel({time_name: slice(it, it+2)}).chunk(
                    {time_name: (2,)})
                # Interpolate
                this_ds = tmp_ds.interp({time_name: [date]}).load()
                flag.loc[date] = InterpFlag.LINEAR

                # Preserve netcdf compression
                if preserve_encoding:
                    for vname, var in tmp_ds.items():
                        this_ds[vname].encoding.update(var.encoding)
       
            # New time
                
        # Append to list
        dss.append(this_ds)
        del this_ds

    # Concatenate
    dso = xr.concat(
        dss, time_name, compat='no_conflicts', coords='different',
        data_vars='all')
    flag.attrs['flags'] = str(InterpFlag)
    if keep_flag:
        dso["flag"] = flag

    # Attributes and encoding
    dso[time_name].attrs.update(ds[time_name].attrs)
    dso[time_name].encoding.update(ds[time_name].encoding)
    if preserve_encoding:
        for vname, var in ds.items():
            dso[vname].encoding.update(var.encoding)
    dso = dso.transpose(time_name,...)

    return dso

def nc_interp_time(ncfiles, dates, ncout, preproc=None, postproc=None):
    """Interpolate netcdf files to a single netcdf file at a given frequency

    ncfiles: str, list
        Single or list of netcdf files
    dates: list
        Output dates
    ncout: str
        Netcdf file with optional date patterns
    postproc: callable
        Callable that operates on the output dataset before saving it

    Return
    ------
    str
        Output netcdf file
    """
    # Input
    dsi = xr_open_concat(ncfiles)

    # Preproc
    if preproc:
        dsi = preproc(dsi)

    # Output dates
    if not isinstance(dates, list):
        dates = [dates]
    dates = [pd.to_datetime(date) for date in dates]
    date0 = dates[0]
    date1 = dates[-1]
    LOGGER.debug("Interpolating to a single file: "
                 f"from {date0:%Y-%m-%dT%H:%M} to {date1:%Y-%m-%dT%H:%M}")
    
    # Interpolate
    dso = interp_time(dsi, dates)

    # Callback postprocessing
    if postproc is not None:
        dso = postproc(dso)
    
    # Write netcdf
    ncfile = ncout.format(**locals())
    
    dso.to_netcdf(ncfile)
    LOGGER.info(f"Created interpolation file: {ncfile}")
    return ncfile

def nc_interp_at_freq_to_daily_nc(
        ncfiles, freq, ncfmt, preproc=None, postproc=None):
    """Interpolate netcdf files to daily netcdf files at a given frequency

    freq: str
        Output frequency
    ncfmt: str
        Netcdf file path format with date patterns

    Return
    ------
    list
        Effective list of output netcdf files
    """
    # Input

    # Preproc

    time = cfs.coords.get(dsi, "time")
    freq = period_vortex_to_pandas(freq)
    date0 = pd.to_datetime(time[0].values).ceil(freq)
    date1 = pd.to_datetime(time[-1].values).ceil(freq)
    day0 = date0.floor("1D")
    day1 = date1.floor("1D")
    
    LOGGER.info("Interpolating to daily files: "
                 f"from {day0:%Y-%m-%dT%H:%M} to {day1:%Y-%m-%dT%H:%M} "
                 f"at frequency {freq}")

    days = pd.date_range(day0, day1, freq='D')

    # Loop on dates
    ncfiles = []
    for iday in range(len(days)):

        # Date
        date = days[iday]
        LOGGER.debug(f"Date: {date}")
        t = time.sel(time=slice(date, None))
        # Hours of the day
        if len(t) == 1:
            times = [date]
        else:
            times = pd.date_range(
                date, date + pd.to_timedelta("24 hours"), 
                freq=freq, closed='left')
        
        # Interpolate
        dso = interp_time(dsi, times, preserve_encoding=True)

        # Callback postprocessing
        if postproc is not None:
            dso = postproc(dso)

        date = pd.to_datetime(date).strftime("%Y%m%d")
        # Save
        ncfile = ncfmt.format(**locals())
        ncfiles.append(ncfile)
        dso.to_netcdf(ncfile)
        LOGGER.info(f"Created daily interpolation file: {ncfile}")
        del dso

    return ncfiles

def nc_interp_at_freq_to_nc(
        ncfiles, freq, ncout, begindate=None, preproc=None, postproc=None):

    freq: str
    ncout: str
    begindate: date
        Start from this date

    Return
    ------
    str
    """
    # Input

    # Preproc

    time = cfs.coords.get(dsi, "time").values
    date0 = pd.to_datetime(begindate or time[0]).ceil(freq)
    date1 = pd.to_datetime(time[-1]).floor(freq)
                 f"from {date0:%Y-%m-%dT%H:%M} to {date1:%Y-%m-%dT%H:%M} "
    dates = pd.date_range(date0, date1, freq=freq)

    dso = dso.astype("float32")

class Regridder(object):
    """Regridder based on the `xesmf` library"""
    def __init__(self, dsi, dso, regridder=None, filename=None,
                 reuse_weights=True, **kwregrid):

        kwregrid.setdefault("method", "bilinear")

        self._dsi = rename_from_standard_name(
            dsi, {'longitude': get_lon_coord_name(dso),
                  'latitude': get_lat_coord_name(dso)})
        self._dso_mapping = {}
        self._dso = rename_from_standard_name(
            dso, {'longitude': get_lon_coord_name(dso),
                  'latitude': get_lat_coord_name(dso)},
            mapping=self._dso_mapping)

        if regridder is not None:
            self.regridder = regridder
        else:
            if not filename or not os.path.exists(filename):
                reuse_weights = False
            self.regridder = xesmf.Regridder(
                self._dsi, self._dso, filename=filename,
                reuse_weights=reuse_weights, **kwregrid)

    def regrid(self, dai, erode=None):
        """
        Regrid a variable

        dai : xarray.DataArray
            Single variable to regrid.

        Returns
        -------
        xarray.DataArray
            Regridded variable.

        """
        # Erode coast
        if erode:
            dai = erode_coast(dai)

        # Appropriate coord names
        dai = rename_from_standard_name(
            dai, {'longitude': 'lon', 'latitude': 'lat'})

        dao = self.regridder(dai).rename(self._dso_mapping)

        return xr.merge([self._dso, dao])

    __call__ = regrid

    def clean(self):
        """
        Remove the weights file
        """
        if hasattr(self, 'regridder'):
            self.regridder.clean_weight_file()

"""
Input/output
"""
import os
from argparse import (
    ArgumentDefaultsHelpFormatter,
    RawDescriptionHelpFormatter,
)

from appdirs import AppDirs
from configobj import ConfigObj

from xoa.cfgm import ConfigManager

from .env import SloopEnvError

class CacheManager(object):

    def __init__(self):
        self.logger = logging.getLogger(__name__)

    def _get_cache_dir_(self):
        try:
            return os.path.join(os.environ['MTOOLDIR'], 'cache',
                                'sloop', 'sessions')
        except SloopEnvError as e:
            self.logger.info(e)
        else:
            return AppDirs(os.path.join("sloop", "sessions"),
                           "$LOGNAME").user_cache_dir

    def get_session_dir(self, session):
        if not os.path.isdir(session):
            session_dir = os.path.join(self._get_cache_dir_(), session)
            self.logger.debug(session_dir)
            if not os.path.isdir(session_dir):
                os.makedirs(session_dir)
            return session_dir
        else:
            return session

    def get_files(self, session, file_format):
        return glob.glob(os.path.join(self._get_cache_dir_(), session,
                         "*." + file_format))

    def _get_session_list_(self, session):
        if session:
            return glob.glob(os.path.join(self._get_cache_dir_(),
                                          session+"*"))
        else:
            return glob.glob(os.path.join(self._get_cache_dir_(), "*"))

    def find_session(self, session=None):
        sess_list = self._get_session_list_(session)
        if len(sess_list) > 1:
            print("Choose one session:" + "\n")
            [print(sess.split('/')[-1] + "\n") for sess in sess_list]
            session_dir = input()
        elif len(sess_list) == 1:
            session_dir = sess_list[0].split('/')[-1]
        else:
            session_dir = None
        return session_dir

    def clear(self, session=None):

        session = self._get_session_list_(session)
            [print(dr + "\n") for dr in session]
            cond = input("Do you really want to erase the session"
                         "s caches listed hereabove ?(yes/no)")
            while not (cond == 'yes' or cond == 'no'):
                cond = input("Do you really want to erase the session"
                             "s caches listed hereabove ?(yes/no)")
            if cond == "yes":
                [
                    subprocess.run(
                     ["rm", "-rf", dr], capture_output=True, check=True,
                    )
                    for dr in session
                ]
        else:
            self.logger.info("The sloop cache is empty")

def tarpat(outfile, directory, patterns=['*']):
    """Create a .tar file given a list of file patterns

        File name of the output .tar file
    directory: str
        Files are considered relative this directory
    patterns: list of str
        List of glob patterns relative to `directory`

    Return
    ------
    list
        List of files that were added to the archive as absolute paths
    """
    if isinstance(patterns, str):
        patterns = [patterns]
    effective = []
    with tarfile.open(outfile, "w") as tar:
            paths = glob.glob(os.path.join(directory, pattern))
            for path in paths:
                relpath = os.path.relpath(path, directory)
                tar.add(path, relpath)
                effective.append(os.path.abspath(path))
    return effective

def xr_open_concat(ncfiles):
    """Helper to concatenate several netcdf file to a single dataset

    This basically a call to :func:`xarray.open_dataset` and
    to :func:`xarray.concat`.

        Single or list of netcdf files.
        In the case of a string (single file), no concatenation a perform.

    Return
    ------
    """
    if isinstance(ncfiles, str):
        return xr.open_dataset(ncfiles)#, chunks={"time": 24})
    #dss = [xr.open_dataset(ncfile, chunks={"time": 24}) for ncfile in ncfiles]
    dss = [xr.open_dataset(ncfile) for ncfile in ncfiles]
    dsout = xr.concat(
        dss, "time", compat="no_conflicts", coords="different",
        data_vars="all"
    )
    dsout = dsout.sortby(dsout["time"])
    for vname, var in dsout.items():
        dsout[vname].encoding.update(dss[0][vname].encoding)
        dsout[vname].attrs.update(dss[0][vname].attrs)    
        dsout[coord].encoding.update(dss[0][coord].encoding)
        dsout[coord].attrs.update(dss[0][coord].attrs)

def nc_get_time(ncfile):
    """Simply return the time axis found in a netcdf file

    ncfile

    Return
    ------
    date
    """
    ds = xr.open_dataset(ncfile)
    time = get_time(ds)
    ds.close()
    return time.dt.floor("D")

class SloopArgumentHelpFormatter(RawDescriptionHelpFormatter,
                                 ArgumentDefaultsHelpFormatter):
    pass

def read_inicfg_files(fini=None, fcfg=None):
    """Read and validate a configobj-like config file"""
    results = ConfigManager(fini).load(fcfg)
    if not results:
        raise SloopError(
                "PLEASE, CHECK [fcfg]"
            )
    return results

def write_fini(dico, fname):
    """Write specification file from dict"""

    config = ConfigObj(dico)
    config.filename = fname
    config.write()

def concat_ascii_files(fin, fout):
    """Concatenation of ascii files"""
    LOGGER.info(f"concatenate ascii files : {fin} + {fout}")
    if not isinstance(fin, list):
        fin = [fin]
    with open(fout, 'w') as outfile:
        for infile in fin:
            with open(infile) as inf:
                outfile.write(inf.read())

"""
Job management
"""
from enum import Enum
import datetime
from .env import is_os_cmd_avail
from .io import CacheManager
import json
import secrets
import configparser

ALLOWED_JOB_MANAGER_TYPES = ["basic", "slurm", "pbspro"]

logger = logging.getLogger(__name__)

cm = CacheManager()

class Job(object):
    """Single job"""

    overview_format = dict(
        name="20",
        jobid="8",
        session="22",
        queue="10",
        realqueue="10",
        time="5",
        status="10",
        submission_date="20",
    )

        self,
        manager,
        name,
        args,
        queue=None,
        jobid=None,
        session=None,
        submission_date=None,
    ):

        self.manager = manager
        self.name = name
        self.queue = queue
        self.jobid = jobid
        self.args = args
        self.status = JobStatus.UNKOWN
        self.realqueue = None
        self.time = None
        self.memory=None
        self.session = session
        self.submission_date = submission_date

    def __str__(self):
        return "Job(name={}, status={}, jobid={}, session={})".format(
            self.name, self.status.name, self.jobid, self.session
        )

    def __repr__(self):
        return "<{}>".format(self)

    def get_status(self):
        return self.manager.get_status(self)  # or JobStatus.FINISHED

    def is_running(self):
        return self.get_status() is JobStatus.RUNNING

    def update_status(self, status=None):
        if status is None:
            status = self.get_status()
            return
        if isinstance(status, JobStatus):
            self.status = status
        else:  # dict
            self.jobid = status["jobid"]
            self.realqueue = status["queue"]
            self.status = status["status"]
            self.time = status["time"]

    @classmethod
    def get_overview_header(cls):
        heads = []
        tails = []
        for name, fmt in cls.overview_format.items():
            name = name.upper()
            heads.append(f"{name:{fmt}}")
            tails.append("-" * len(heads[-1]))
        return "      ".join(heads) + "\n" + "    ".join(tails)

    def get_overview(self, update=True):
        if update:
            self.update_status()
        name = self.name
        jobid = self.jobid
        queue = self.queue
        realqueue = self.realqueue
        status = self.status.name
        session = self.session
        submission_date = self.submission_date
        if self.time is not None:
            hours = self.time.seconds // 3600
            minutes = (self.time.seconds - hours * 3600) // 60
            time = f"{hours:02}h{minutes:02}"
        else:
            time = "--h--"
        fmt = "    ".join(
            [
                ("{" + key + f"!s:{ff}" + "}")
                for key, ff in self.overview_format.items()
            ]
        )
        return fmt.format(**locals())
    
class JobStatus(Enum):
    UNKOWN = -1
    FINISHED = 0
    PENDING = 1
    RUNNING = 2
    INQUEUE = 3
    EXITING = 4
    COMPLETING = 5

class BasicJobManager(object):
    """Basic job manager"""
    commands = {
        "submit": {
            "command": "python",
            "options": {
                "script": "{}",
                "log_out": "-o {}",
            },
        },
    }

    status_names = {
        "F": JobStatus.FINISHED,
    }
    
        self.jobs = []
        self.session = None
        logger.debug("Started job manager: " + self.__class__.__name__)

    def from_type(jmtype):
        jmtype = jmtype.lower()
        assert jmtype in ALLOWED_JOB_MANAGER_TYPES
        cls_name = jmtype.title() + "JobManager"
        from . import job

        return getattr(job, cls_name)()

    def _get_command_args_(cls, command, **kwargs):
        args = []
        if "command" in cls.commands[command]:
            args.append(cls.commands[command]["command"])
        if "options" in cls.commands[command]:
            for oname, ovalue in kwargs.items():
                if oname in cls.commands[command]["options"].keys():
                    if not type(ovalue) == type(None):
                        fmt = cls.commands[command]["options"][oname]
                        if isinstance(ovalue, list):
                            ovalue = [val for val in ovalue if val]
                            if not isinstance(
                                cls.commands[command]["options"][oname], tuple
                            ):
                                sep = ","
                            else:
                                fmt, sep = fmt
                            for val in ovalue:
                                args.append(fmt.format(val))
                        else:
                            fmt = fmt.format(ovalue).split(" ")
                            args += fmt
        return args

    def _get_opts_(cls, command, opts_in):
        sep = ","
        opts = {}
        if cls.commands[command]['options']:
            for context in cls.commands[command]['options'].keys():
                try:
                    value = opts_in[context]
                    if sep in value:
                        value = value.split(sep)
                except:
                    continue
                else:
                    opts[context] = value
        return opts

    def _parse_submit_res_(self, res, jobargs, session):
        if res.stderr:
            logger.debug(
                "Job submit stderr: "
                + res.stderr.decode("utf-8", errors="ignore")
            )
        if res.stdout:
                "Job submit stdout: "
                + res.stdout.decode("utf-8", errors="ignore")
            )
        job = Job(
            manager=self,
            name=jobargs["name"],
            queue=jobargs["queue"],
            args=res.args,
            jobid=None,
            session=session,
            submission_date=str(datetime.datetime.now())[:-7],
        )
        return job
        
    def _get_session_id_(self):
        return secrets.token_hex(8)
        
    def to_json(self, job):
        session_dir = cm.get_session_dir(self.session)
        jobdict = self.to_dict(job)
        with open(
            session_dir + "/" + jobdict["name"] + ".json", "w"
        ) as outfile:
            json.dump(jobdict, outfile, indent=4)
        
    def from_json(self, json_files):
        for json_file in json_files:
            with open(json_file) as jsonf:
                content = json.load(jsonf)
            job = Job(
                manager=self,
                name=content["name"],
                args=content["args"],
                jobid=content["jobid"],
                queue=content["queue"],
                session=content["session"],
                submission_date=content["submission_date"],
            )
            self.jobs.append(job)
        return self.jobs

    def to_dict(self, job):
        dict_job = dict()
        for (key, value) in job.__dict__.items():
                if isinstance(value, str):
                    if value != "":
                        dict_job[key] = value
                elif key == 'manager':
                      dict_job[key] = value.__class__.__name__
                elif key == "status":
                    dict_job[key] = str(value.name)
                elif key == "time":
                    if value is not None:
                        hours = value.seconds // 3600
                        minutes = (value.seconds - hours * 3600) // 60
                        dict_job[key] = f"{hours:02}h{minutes:02}"
                    else:
                        dict_job[key] = "--h--"
                else:
                    dict_job[key] = value

        return dict_job

    def _parse_status_res_(self, res):
                "Job status stderr: "
            )
                "Job status stdout: "
            )
        return res.stdout.decode("utf-8", errors="ignore")

    def get_jobids(self, name=None, queue=None):
        if isinstance(name, Job):
            assert name in self.jobs
            jobs = [name]
        else:
            jobs = self.get_jobs(name=name, queue=queue)
        jobids = [job.jobid for job in jobs if job.jobid is not None]
        return jobs, jobids

    def update_status(self, name=None, queue=None, jobids=None):
        if self.__class__.__name__ != "BasicJobManager" :
            if jobids is None:
                jobs, jobids = self.get_jobids(name=name, queue=queue)
            args = self._extra_status_args_(self._get_command_args_(
                    "status", jobid=self._jobid_sep_().join(jobids)
                    ))
            if args:
                logger.debug("Get status: " + " ".join(args))
                res = subprocess.run(args, capture_output=True, check=True)
                logger.debug("Got status")
                status_list = self._parse_status_res_(res)
            if status_list:
                for status in status_list:
                    self.get_job(status["jobid"]).update_status(status)
        else:
            [job.update_status(JobStatus.FINISHED) for job in self.jobs]

    def get_status(self, name=None, queue=None):
        jobs = self.update_status(name=name, queue=queue)
        if jobs:
            if not isinstance(jobs, list):
                return jobs.status
            return [job.status for job in jobs]

    def get_job(self, jobid):
        if jobid is None:
            return
        for job in self.jobs:
            if job.jobid == jobid:
                return job

    def get_jobs(self, name=None, queue=None):
        jobs = []
        if name:
            if len(name)>1:
                for job in self.jobs:
                    if job.name in name:
                        jobs.append(job)
        else:    
            for job in self.jobs:
                if (name is not None and job.name != name) or (
                    queue is not None and job.queue != queue
                ):
                    continue
                jobs.append(job)
        return jobs

    def _get_jobs_session_(self, session=None):
        self.session = cm.find_session(session)
        json_files = cm.get_files(self.session, 'json')
        self.from_json(json_files)

    def __getitem__(self, name):
        return self.get_jobs(name=name)

    def get_overview(self,jobids=None):
        self.update_status(jobids=jobids)
        header = Job.get_overview_header()
        overviews = [job.get_overview(update=False) for job in self.jobs]
        return header + "\n" + "\n".join(overviews)

        return self.get_overview()    
        
    def submit(self, opts):
        self.session = opts["session"]
        script = f"{opts['job']}"
        opts = self._get_opts_('submit', opts)
        if "queue" not in opts.keys():
            opts.update({"queue":None})
        jobargs = {}
        if "depend" in opts.keys():
            opts["depend"] = ":".join(opts["depend"])
        jobargs.update(opts)
        jobargs.update(dict(script=script))
        args = self._get_command_args_("submit", **jobargs)
        logger.debug("Submit: " + " ".join(args))
        res = subprocess.run(args, capture_output=True, check=True)
        logger.debug("Submitted")
        job = self._parse_submit_res_(res, jobargs, self.session)
        self.jobs.append(job)
        self.to_json(job)
        self.check_status(session=self.session, show=False)
        return job.jobid
         
    def check_status(self, session, show=True):
        self._get_jobs_session_(session=session)
        self.get_overview(jobids=[job.jobid for job in self.jobs])
        [self.to_json(job) for job in self.jobs]
        if show == True:
            print(self)

    def delete(self, session=None):
        self.check_status(session=session)
        cond = input(
            "Do you really want to delete the jobs listed hereabove ?(yes/no)"
        )
        if cond == "yes":
            args = self._get_command_args_(
                "delete", jobid=self._jobid_sep_().join([job.jobid for job in self.jobs])
            )
            subprocess.run(args, capture_output=True, check=True)

    def delete_force(self, session=None):
        args = self._get_command_args_(
            "delete", jobid=self._jobid_sep_().join([job.jobid for job in self.jobs])
            )
        subprocess.run(args, capture_output=True, check=True)
    
class PbsJobManager(BasicJobManager):
    """Pbs Job Manager"""
            "command": "qsub",
                "name": "-N {}",
                "queue": "-V -q {}",
                "time": "-l walltime={}",
                "memory":"-l mem={}",
                "log_out": "-koed -o {}",
                "depend": ("-W depend=afterok:{}"),
                "mail": "-M {}",
            },
        },
        "status": {
            "command": "qstat",
            "options": {"jobid": "{}", "logname": "-u $LOGNAME",},
        },
        "delete": {"command": "qdel", "options": {"jobid": "{}",}},
    }

        "R": JobStatus.RUNNING,
        "E": JobStatus.EXITING,
        "Q": JobStatus.INQUEUE,
        "H": JobStatus.PENDING,
    }

        jobid = BasicJobManager.submit(self, opts)
        return jobid

        job = BasicJobManager._parse_submit_res_(self, res, jobargs, session)
        job.jobid = res.stdout.decode("utf-8", errors="ignore").split(".")[0]
        return job

    def _jobid_sep_(self):
        return ' '

    def _extra_status_args_(self,args):
        args.append('-x')
        args.append('-u $LOGNAME')

        """JOBID PARTITION NAME USER ST TIME NODES NODELIST(REASON)"""
        res = BasicJobManager._parse_status_res_(self, res)
        lines = res.splitlines()[5:]
        out = []
        for line in lines:
            (
                jobid,
                user,
                queue,
                name,
                session,
                nodes,
                task,
                mem,
                time,
                status,
                elaptime,
            ) = line.split()
            if (elaptime == "--:--") or elaptime == "--":
                elaptime = None
            else:
                hms = elaptime.split(":")
                hh = 0
                mm = 0
                ss = 0
                if len(hms) == 1:
                    ss = hms[0]
                elif len(hms) == 2:
                    hh, mm = hms
                else:
                    hh, mm, ss = hms
                elaptime = datetime.timedelta(
                    seconds=int(ss), minutes=int(mm), hours=int(hh)
                )
            if status in self.status_names:
                status = self.status_names[status]
            else:
                status = JobStatus.UNKOWN
            out.append(
                {
                    "jobid": jobid.split(".")[0],
                    "queue": queue,
                    "name": name,
                    "time": elaptime,
                    "status": status,
                }
            )
        return out

class SlurmJobManager(BasicJobManager):
    """Slurm Job Manager"""
            "command": "sbatch",
                "name": "--exclusive -J {}",
                "queue": "-p {}",
                "nnodes": "-N {}",
                "ncpus": "-c {}",
                "mem": "--mem={}",
                "time": "--time={}",
                "depend": "--dependency=afterok:{}",
                "mail": "--mail-type=ALL --mail-user={}",
            },
        },
            "command": "squeue",
                "jobid": "--jobs={}",
                "queue": "--partition={}",
                "name": "--name={}",
                "users": "--users={}",
                "noheader": "--noheader",
            },
        },
        "delete": {"command": "scancel", "options": {"jobid": "{}",}},
    }

        "CD": JobStatus.FINISHED,
        "PD": JobStatus.PENDING,
        "CG": JobStatus.COMPLETING,
    }

        args.append('--noheader')

        return ','

        job.jobid = res.stdout.decode("utf-8", errors="ignore").split()[-1]
        return job

        res=BasicJobManager._parse_status_res_(self, res)
        out = []
        lines = res.splitlines()
        if lines:
            for line in lines:
                info = line.split()
                jobid = info[0]
                queue = info[1]
                name = info[2]
                user = info[3]
                status = info[4]
                time = info[5]
                nodes = info[6]
                nodelist = " ".join(info[7:])
 
                hms = time.split(":")
                hh = 0
                mm = 0
                    ss = hms[0]
                    mm, ss = hms
                else:
                time = datetime.timedelta(seconds=int(mm), minutes=int(hh))
                if status in self.status_names:
                    status = self.status_names[status]
                else:
                    status = JobStatus.UNKOWN
                out.append(
                    {"jobid": jobid, "queue": queue, "name": name, "time": time, "status": status}
                )
        return out

if( is_os_cmd_avail("squeue") ):  
    jobmanager = SlurmJobManager()  
elif( is_os_cmd_avail("qsub") ):  
    jobmanager = PbsJobManager()  
else:  
    jobmanager = BasicJobManager()

"""
if __name__ == '__main__':
    jm = BasicJobManager()
    #job = Job(jm, name='toto', args='mkjob.py rundate=10', queue='seq',jobid='007')
    job = Job(jm, name='fake_task', args='fake_task.py rundate=10', queue='None')
    jm.jobs.append(job)
    print(job.get_status())
    #print(job)
    print(Job.get_overview_header())
    #print(job.get_overview())
    #print(jm)
"""

import logging.config

# import colorlog

DEFAULT_LOGGING_CONFIG = {
    "version": 1,
    #    'incremental': True,
    "disable_existing_loggers": False,
    "formatters": {
        "brief": {
            "format": ("%(log_color)s%(name)-12s: %(bold)s%(levelname)-8s"
                       "%(bold)s%(reset)s%(log_color)s %(message)s%(reset)s"),
            "()": "colorlog.ColoredFormatter",
            "log_colors": {
                "DEBUG": "thin",
                "INFO": "green",
                "WARNING": "yellow",
                "ERROR": "red",
                "CRITICAL": "red,bg_white",
            },
        },
        "brief_no_color": {
            "format": "%(name)-12s: %(levelname)-8s %(message)s",
        },
        "precise": {
            "format": "%(asctime)s %(name)-12s: %(levelname)-8s %(message)s",
        },
    },
    "handlers": {
        "console": {
            "class": "logging.StreamHandler",
            "formatter": "brief",
            "level": "INFO",
        },
        "file": {
            "class": "logging.handlers.RotatingFileHandler",
            "formatter": "precise",
            "filename": "sloop.log",
            "maxBytes": 1024 * 1000,
            "backupCount": 3,
            "level": "DEBUG",
        },
    },
    "root": {"handlers": ["console", "file"], "level": "DEBUG"},
}

def setup_logging(
    console_level=None, file_level=None, to_file=True, no_color=False,
    show_init_msg=False, **kwargs
):
    """Setup the logging"""

    #    for handler in logging.root.handlers:
    #        logging.root.handlers.remove(handler)
    #    del logging.root.handlers[:]

    # Alter the config
    logging_config = DEFAULT_LOGGING_CONFIG.copy()
    if console_level is not None:
        logging_config["handlers"]["console"]["level"] = console_level.upper()
    if to_file is False:
        logging_config["root"]["handlers"].remove("file")
    else:
        fconfig = logging_config["handlers"]["file"]
        if isinstance(to_file, str):
            fconfig["filename"] = to_file
        if file_level is not None:
            fconfig["level"] = file_level.upper()
    if no_color:
        logging_config["handlers"]["console"]["formatter"] = "brief_no_color"

    # Update config
    logging_config.update(kwargs)

    # Load it
    logging.config.dictConfig(logging_config)
    if show_init_msg:
        logging.getLogger(__name__).debug(
            "*** INITIALISED LOGGING SESSION ***")

def add_logging_parser_arguments(parser, default_level="info"):
    add_log_level_parser_arguments(parser, default_level=default_level)
    parser.add_argument(
        "--log-file", help="logging file name", default="sloop.log"
    )
        "--log-no-color",
        help="suppress colors in console",
        action="store_true",
    )

def add_log_level_parser_arguments(parser, default_level="info"):
        "--log-level",
        help="logging level",
        choices=["debug", "info", "warning", "error", "critical"],
        default=default_level,
    )

def main_setup_logging(args):
    setup_logging(
        console_level=args.log_level,
        to_file=args.log_file,
        no_color=args.log_no_color,
    )

class LoggingLevelContext(object):
    """Temporarily change the logging level within a context

    Example
    -------
    >>> with LoggingLevelContext(logging.WARNING) as logger:
    >>> ... logger.warning("warning level here")
    """

    def __init__(self, level, logger=None):
        self.logger = logger if logger else logging.root
        self.level = level

        self.old_level = self.logger.level
        self.logger.level = self.level
        return self.logger

        self.logger.level = self.old_level

def celsius2kelvin(ds):
    """Convert Celsius degree to Kelvin

    Returns
    -------

    """
    LOGGER.info("Convert temperature from Celsius to Kelvin degrees")
    ds['t2m'] -= 273.15
    return ds

def windstress(ds, method=None):
    """Compute wind stress

        ds : xarray.Dataset
        method : str         
        
    Returns
    -------
        xarray.Dataset

    """
    LOGGER.info("Compute the wind stress following the "
        "{method} method")
    if method == 'Speich':
        rhocd = 1.3*(0.61+0.063*ds.si10)*1e-3
        taux = rhocd * ds.si10 * ds.u10
        tauy = rhocd * ds.si10 * ds.v10
        taux.name = 'ewss'
        tauy.name = 'nsss'
        return xr.merge([ds, taux, tauy], compat='override')

def radiativeflux(ds):
    """Compute radiative flux

    Returns
    -------

    """
    LOGGER.info("Compute the radiative flux")
    ds['radflx'] = ds.ssr + ds.str
    return ds

def watervapormixingratio(ds):
    """Compute water vapor mixing ratio

    Returns
    -------

    """
    LOGGER.info("Compute the water vapor mixing ratio")
    t = ds.t2m
    qsatur =.622e-3*(6.107799961e+00+t*(4.436518521e-01
                    +t*(1.428945805e-02+t*(2.650648471e-04
                    +t*(3.031240396e-06+t*(2.034080948e-08
                    +t* 6.136820929e-11))))))
    ds['vapmix'] = qsatur * ds.r2 * 0.01
    return ds

"""
Time utilities
"""

def convert_to_julian_day(time):
    """Conversion of date in Julian day"

    time

    Output
    ------
    xr.Dataset
    """
    LOGGER.info("Convert time to julian day")
    return (time - np.datetime64("1950-01-01")) / np.timedelta64(1, "D")

def convert_from_julian_day(time, unit='D'):
    """Conversion of Julian day in date"

    time

    Output
    ------
    """
    LOGGER.info("Convert julian day to time")
    res = np.timedelta64(int(time), unit) + np.datetime64("1950-01-01")
    return res

def find_monday(date):
    """Find the lastest monday

    date: string format

    Output
    ------
    date of the latest monday in string format
    """
    date = pd.to_datetime(date)
    monday = date - pd.to_timedelta(date.dayofweek, "D")
    return monday.strftime("%Y%m%d")

def running_time(start=None, maxterm=None, step=None):
    """Build the list of running dates

    start: starting date
    maxterm: maximum term
    step: time step

    Output
    ------
    list of running dates
    """
    step = period_vortex_to_pandas(step)
    end = pd.to_datetime(start)+pd.to_timedelta(maxterm)
    time = pd.date_range(start=start, end=end.strftime('%Y%m%d%H%M'),
                         freq=step)

    return list(time)

def period_vortex_to_pandas(period):
    """Convert a vortex period string to a pandas period string

    For pandas, see : https://pandas.pydata.org/pandas-docs/stable/user_guide/timeseries.html#timeseries-offset-aliases
    """
    if period.startswith("PT"):
        if period.endswith("M"):
            period = period.replace('M', 'min')
        period = period[2:]
    elif period.startswith("P"):
        period = period[1:]
    return period