Changeset 6190 for TOOLS


Ignore:
Timestamp:
07/06/22 11:06:07 (22 months ago)
Author:
omamce
Message:

O.M. : Evolution on MOSAIX

  • Add licensing information
  • Update and beautify README.md
  • Add few model version to known list
Location:
TOOLS/MOSAIX
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/MOSAIX/Build_coordinates_mask.py

    r6124 r6190  
    66# creates file coordinates_mask.nc used by mosaix from NEMO netcdf files  
    77# coordinates.nc and bathymetry.nc (on the same grid as coordinates.nc) 
    8  
     8##  
     9##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     10##  file for an english version of the licence and 
     11##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     12## 
     13##  Permission is hereby granted, free of charge, to any person or 
     14##  organization obtaining a copy of the software and accompanying 
     15##  documentation covered by this license (the "Software") to use, 
     16##  reproduce, display, distribute, execute, and transmit the 
     17##  Software, and to prepare derivative works of the Software, and to 
     18##  permit third-parties to whom the Software is furnished to do so, 
     19##  all subject to the following: 
     20## 
     21##  Warning, to install, configure, run, use any of MOSAIX software or 
     22##  to read the associated documentation you'll need at least one (1) 
     23##  brain in a reasonably working order. Lack of this implement will 
     24##  void any warranties (either express or implied).  Authors assumes 
     25##  no responsability for errors, omissions, data loss, or any other 
     26##  consequences caused directly or indirectly by the usage of his 
     27##  software by incorrectly or partially configured 
     28## 
     29## 
    930import numpy  as np 
    1031import xarray as xr 

  • TOOLS/MOSAIX/CalvingWeights.py

    r6105 r6190  
    55### =========================================================================== 
    66## 
    7 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    8 ##  software or to read the associated documentation you'll need at least 
    9 ##  one (1) brain in a reasonably working order. Lack of this implement 
    10 ##  will void any warranties (either express or implied). 
    11 ##  O. Marti assumes no responsability for errors, omissions, 
    12 ##  data loss, or any other consequences caused directly or indirectly by 
    13 ##  the usage of his software by incorrectly or partially configured 
    14 ##  personal. 
     7##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     8##  file for an english version of the licence and 
     9##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     10## 
     11##  Permission is hereby granted, free of charge, to any person or 
     12##  organization obtaining a copy of the software and accompanying 
     13##  documentation covered by this license (the "Software") to use, 
     14##  reproduce, display, distribute, execute, and transmit the 
     15##  Software, and to prepare derivative works of the Software, and to 
     16##  permit third-parties to whom the Software is furnished to do so, 
     17##  all subject to the following: 
     18## 
     19##  Warning, to install, configure, run, use any of MOSAIX software or 
     20##  to read the associated documentation you'll need at least one (1) 
     21##  brain in a reasonably working order. Lack of this implement will 
     22##  void any warranties (either express or implied).  Authors assumes 
     23##  no responsability for errors, omissions, data loss, or any other 
     24##  consequences caused directly or indirectly by the usage of his 
     25##  software by incorrectly or partially configured 
     26## 
    1527## 
    1628import numpy as np, xarray as xr 
    … …  
    3749# Adding arguments 
    3850parser.add_argument ('--oce'     , help='oce model name', type=str, default='eORCA1.2', 
    39                          choices=['ORCA2.3', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 
     51                         choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 
    4052parser.add_argument ('--atm'     , help='atm model name (ICO* or LMD*)', type=str, default='ICO40'    ) 
    4153parser.add_argument ('--type'    , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full'  ) 

  • TOOLS/MOSAIX/ComputeNemoCoast.py

    r3901 r6190  
    55### =========================================================================== 
    66## 
    7 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    8 ##  software or to read the associated documentation you'll need at least 
    9 ##  one (1) brain in a reasonably working order. Lack of this implement 
    10 ##  will void any warranties (either express or implied). 
    11 ##  O. Marti assumes no responsability for errors, omissions, 
    12 ##  data loss, or any other consequences caused directly or indirectly by 
    13 ##  the usage of his software by incorrectly or partially configured 
    14 ##  personal. 
     7##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     8##  file for an english version of the licence and 
     9##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     10## 
     11##  Permission is hereby granted, free of charge, to any person or 
     12##  organization obtaining a copy of the software and accompanying 
     13##  documentation covered by this license (the "Software") to use, 
     14##  reproduce, display, distribute, execute, and transmit the 
     15##  Software, and to prepare derivative works of the Software, and to 
     16##  permit third-parties to whom the Software is furnished to do so, 
     17##  all subject to the following: 
     18## 
     19##  Warning, to install, configure, run, use any of MOSAIX software or 
     20##  to read the associated documentation you'll need at least one (1) 
     21##  brain in a reasonably working order. Lack of this implement will 
     22##  void any warranties (either express or implied).  Authors assumes 
     23##  no responsability for errors, omissions, data loss, or any other 
     24##  consequences caused directly or indirectly by the usage of his 
     25##  software by incorrectly or partially configured 
     26## 
    1527## 
    1628## SVN information 

  • TOOLS/MOSAIX/CreateOasisGrids.bash

    r6107 r6190  
    66### =========================================================================== 
    77## 
    8 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    9 ##  software or to read the associated documentation you'll need at least 
    10 ##  one (1) brain in a reasonably working order. Lack of this implement 
    11 ##  will void any warranties (either express or implied). 
    12 ##  O. Marti assumes no responsability for errors, omissions, 
    13 ##  data loss, or any other consequences caused directly or indirectly by 
    14 ##  the usage of his software by incorrectly or partially configured 
    15 ##  personal. 
     8##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     9##  file for an english version of the licence and 
     10##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     11## 
     12##  Permission is hereby granted, free of charge, to any person or 
     13##  organization obtaining a copy of the software and accompanying 
     14##  documentation covered by this license (the "Software") to use, 
     15##  reproduce, display, distribute, execute, and transmit the 
     16##  Software, and to prepare derivative works of the Software, and to 
     17##  permit third-parties to whom the Software is furnished to do so, 
     18##  all subject to the following: 
     19## 
     20##  Warning, to install, configure, run, use any of MOSAIX software or 
     21##  to read the associated documentation you'll need at least one (1) 
     22##  brain in a reasonably working order. Lack of this implement will 
     23##  void any warranties (either express or implied).  Authors assumes 
     24##  no responsability for errors, omissions, data loss, or any other 
     25##  consequences caused directly or indirectly by the usage of his 
     26##  software by incorrectly or partially configured 
    1627## 
    1728### 

  • TOOLS/MOSAIX/CreateWeightsMask.bash

    r6107 r6190  
    2121### =========================================================================== 
    2222## 
    23 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    24 ##  software or to read the associated documentation you'll need at least 
    25 ##  one (1) brain in a reasonably working order. Lack of this implement 
    26 ##  will void any warranties (either express or implied). 
    27 ##  O. Marti assumes no responsability for errors, omissions, 
    28 ##  data loss, or any other consequences caused directly or indirectly by 
    29 ##  the usage of his software by incorrectly or partially configured 
    30 ##  personal. 
     23##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     24##  file for an english version of the licence and 
     25##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     26## 
     27##  Permission is hereby granted, free of charge, to any person or 
     28##  organization obtaining a copy of the software and accompanying 
     29##  documentation covered by this license (the "Software") to use, 
     30##  reproduce, display, distribute, execute, and transmit the 
     31##  Software, and to prepare derivative works of the Software, and to 
     32##  permit third-parties to whom the Software is furnished to do so, 
     33##  all subject to the following: 
     34## 
     35##  Warning, to install, configure, run, use any of MOSAIX software or 
     36##  to read the associated documentation you'll need at least one (1) 
     37##  brain in a reasonably working order. Lack of this implement will 
     38##  void any warranties (either express or implied).  Authors assumes 
     39##  no responsability for errors, omissions, data loss, or any other 
     40##  consequences caused directly or indirectly by the usage of his 
     41##  software by incorrectly or partially configured 
     42## 
    3143## 
    3244## SVN information 
    … …  
    7082echo ${Titre}"Defines model"${Norm} 
    7183# ================================= 
    72 #CplModel=ORCA2.3xLMD9695 
     84CplModel=ORCA2.3xLMD9695 
    7385#CplModel=ORCA2.3xICO30 
    7486#CplModel=ORCA2.3xICO40 
    … …  
    7890#CplModel=eORCA1.2xLMD256256 
    7991#CplModel=eORCA1.2xICO40 
    80 CplModel=eORCA1.4.2xICO40 
     92#CplModel=eORCA1.4.2xICO40 
    8193#CplModel=eORCA1.2xICO450 
    8294#CplModel=eORCA025.1xLMD256256 
    … …  
    898910echo ${Titre}"Save results"${Norm} 
    899911## =========================================================================== 
    900 for File in $(ls dia_*.nc rmp_*.nc ${ATM}_grid_maskFrom_${OCE}*.nc areas_${OCE}x${ATM}*.nc grids_${OCE}x${ATM}*.nc masks_${OCE}x${ATM}*.nc 2> /dev/null ) ; do 
    901     cp ${File} ${SUBMIT_DIR} 
    902 done 
     912#for File in $(ls dia_*.nc rmp_*.nc ${ATM}_grid_maskFrom_${OCE}*.nc areas_${OCE}x${ATM}*.nc grids_${OCE}x${ATM}*.nc masks_${OCE}x${ATM}*.nc 2> /dev/null ) ; do 
     913#    cp ${File} ${SUBMIT_DIR} 
     914#done 
    903915 
    904916 

  • TOOLS/MOSAIX/README.md

    r6044 r6190  
    1616or `svn co svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/MOSAIX MOSAIX` 
    1717 
    18 ### [XIOS](http://forge.ipsl.jussieu.fr/ioserver/) 
     18## [XIOS](http://forge.ipsl.jussieu.fr/ioserver/) 
    1919MOSAIC is known to work with XIOS revision 2286 
    2020 
    … …  
    2323Extraction : svn co http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/trunk XIOS 
    2424 
     25## Licensing 
     26 
     27MOSAIX is under CeCILL\_V2 licence. See `Licence_CeCILL_V2-en.txt`file 
     28 for an english version of the licence and `Licence_CeCILL_V2-fr.txt` 
     29 for a french version. 
     30  
     31Permission is hereby granted, free of charge, to any person or 
     32organization obtaining a copy of the software and accompanying 
     33documentation covered by this license (the "Software") to use, 
     34reproduce, display, distribute, execute, and transmit the Software, 
     35and to prepare derivative works of the Software, and to permit 
     36third-parties to whom the Software is furnished to do so, all subject 
     37to the following warning : 
     38 
     39### Warning 
     40 
     41Warning, to install, configure, run, use any of MOSAIX software or to 
     42read the associated documentation you'll need at least one (1) brain 
     43in a reasonably working order. Lack of this implement will void any 
     44warranties (either express or implied). Authors assumes no 
     45responsability for errors, omissions, data loss, or any other 
     46consequences caused directly or indirectly by the usage of his 
     47software by incorrectly or partially configured 
     48 
    2549## MOSAIX 
    2650 
    27 MOSAIX is a small set of fortran routines, bash scripts, python scripts and xml files allowing to generate coupling weights for IPSL Earth System Model. It is under development, and presently not used by the production runs. It aims to replace the old MOSAIC software. It uses XIOS as a library to compute weights. In the present state, it can handle ORCA, LMDZ lon/lat and DYNAMICO grids. 
    28  
    29 MOSAIX can generate 2<sup>nd</sup> order conservative interpolation weights. They are used in place of bilinear interpolation by [OASIS-MCT](https://portal.enes.org/oasis) 
     51MOSAIX is a small set of fortran routines, bash scripts, python 
     52scripts and xml files allowing to generate coupling weights for the 
     53IPSL Earth System Model. It is under development, and presently not 
     54used by the production runs. It aims to replace the old MOSAIC 
     55software. It uses XIOS as a library to compute weights. In the present 
     56state, it can handle ORCA, LMDZ lon/lat and DYNAMICO grids. 
     57 
     58MOSAIX can generate second order conservative interpolation 
     59weights.  They are used in place of bilinear interpolation by 
     60[OASIS-MCT](https://portal.enes.org/oasis) 
    3061 
    3162As MOSAIX uses XIOS, it works in parallel, using MPI distributed memory. 
    … …  
    3364## Compatibility with MOSAIC 
    3465 
    35 -   Weights are not compatible with the weights generated by MOSAIC. The precision of the interpolation (conservation) is improved. The ocean mask interpolated on the atmosphere grid will be slightly different, and not compatible. 
    36  
    37 -   **Run-off** : runoff weights are computed with a different algorithm, the result can be quite different. 
    38  
    39 -   **Calving** : In IPSLCM6, the calving was integrated for each latitude of the grid on the LMDZ grid. Weights were use to aggregate several latitudes to compute an integral over a latitudinal band.. This is not possible with DYNAMICO, and was abandon for rectilinear LMD grid. The atmosphere model must the integral over a region (presently latitudinal bands). The final result should be the same as with MOSAIC. 
     66-   Weights are not compatible with the weights generated by 
     67    MOSAIC. The precision of the interpolation (conservation) is 
     68    improved. The ocean mask interpolated on the atmosphere grid will 
     69    be slightly different, and not compatible. 
     70 
     71-   **Run-off** : runoff weights are computed with a different 
     72    algorithm, the result can be quite different. 
     73 
     74-   **Calving** : In IPSLCM6, the calving was integrated for each 
     75    latitude of the grid on the LMDZ grid. Weights were use to 
     76    aggregate several latitudes to compute an integral over a 
     77    latitudinal band. This is not possible with DYNAMICO, and was 
     78    abandon for rectilinear LMD grid. The atmosphere model must the 
     79    integral over a region (presently latitudinal bands). The final 
     80    result should be the same as with MOSAIC. 
     81         
     82## Computer 
     83 
     84MOSAIX been developped and extensively used on Irene Skylake computer, 
     85using the Intel compilation tools. MOSAIX is supposed to work on any 
     86computer where XIOS works. 
     87 
     88It has been tested on Mac OS X, using `gcc` as a compiler. It is 
     89fragile, as regurlarly XIOS updates break the compatibility. 
    4090 
    4191## TODO 
    4292 
    4393-   Creates a non-masked grid for atmosphere, to compute integrals of run-off. 
    44 -   Creates 2<sup>nd</sup> order interpolation weights that conserve posivity. They won't be conservative. 
     94-   Creates 2<sup>nd</sup> order interpolation weights that conserve 
     95    posivity. They won't be conservative. 
     96-   Test on other computers : Irene AMD, Jean-Zay, etc ... 
    4597 
    4698## Known problems 
    … …  
    65117 
    66118-   Fortran compiler, C++ compiler, MPI. 
    67 -   A working [XIOS](http://forge.ipsl.jussieu.fr/ioserver) library. Revision 1615 is known to work. 
     119-   A working [XIOS](http://forge.ipsl.jussieu.fr/ioserver) 
     120    library. Revision 1615 is known to work. 
    68121-   fcm version 1.2 
    69 -   NetCDF fortran librairies. See [XIOS](http://forge.ipsl.jussieu.fr/ioserver) for a suitable version. 
    70 -   Python, version>=3.6 with [netCDF4](http://unidata.github.io/netcdf4-python/) and [numpy](http://www.numpy.org) modules. Python version>=2.7 mights work, but has not been tested recently. 
     122-   NetCDF fortran librairies. See 
     123    [XIOS](http://forge.ipsl.jussieu.fr/ioserver) for a suitable version. 
     124-   Python, version>=3.6 with 
     125    [netCDF4](http://unidata.github.io/netcdf4-python/) and 
     126    [numpy](http://www.numpy.org) modules. Python version>=2.7 mights work, but has not been tested recently. 
    71127-   Bash version>=4. 
    72128-   [​NCO](http://nco.sourceforge.net). Version 4.6.0 is known to work. 
    … …  
    74130## Input files 
    75131 
    76 MOSAIX needs files describing the grid of the models. See `CreateWeightsMask.bash` and the bottom of this page. 
     132MOSAIX needs files describing the grid of the models. See 
     133`CreateWeightsMask.bash` and the bottom of this page. 
    77134 
    78135## Output files 
    … …  
    80137All netCDF files have a detailed header. 
    81138 
    82 -   Weight files `rmp_*`. Should have a rather explicit name. The name is configurable at the beginning of `CreateWeightsMask.bash`. 
    83 -   Diagnostic file used for debug `dia_*.nc` for each rmp file, except for runoff and calving. Debug information are directly in weight files in this last cases. 
    84 -   Grid description file for OASIS-MCT `grids.nc`, `areas.nc` and `masks.nc`. 
    85 -   A file containing the fraction of ocean in each atmosphere grid box. For instance `ICO40_grid_maskFrom_eORCA1.2.nc`. 
     139-   Weight files `rmp_*`. Should have a rather explicit name. The name 
     140    is configurable at the beginning of `CreateWeightsMask.bash`. 
     141-   Diagnostic file used for debug `dia_*.nc` for each rmp file, 
     142    except for runoff and calving. Debug information are directly in 
     143    weight files in this last cases. 
     144-   Grid description file for OASIS-MCT `grids.nc`, `areas.nc` and 
     145    `masks.nc`. 
     146-   A file containing the fraction of ocean in each atmosphere grid 
     147    box. For instance `ICO40_grid_maskFrom_eORCA1.2.nc`. 
    86148-   A `README` file. Includes a checksum for each generated files. 
    87149 
    88150## Scripts and codes 
    89151 
    90 *   **CreateWeightsMask.bash**. Creates the weights to interpolate between atmospheric and ocean grid. Weight files are suitable for OASIS-MCT. Creates the fraction of ocean in atmospheric grid cells. Uses parallel processing with MPI. Has a configuration section at the beginning. 
    91  
    92 -   **CreateOasisGrids.bash**. Creates files `grids.nc`, `areas.nc` and `masks.nc` with information from both models, suitable for OASIS-MCT. Mono CPU. Launched by `CreateWeights.bash`. 
    93  
    94 -   **update_xml.py**. Python script used to perform on the fly editing of xml files. Used by `CreateWeights.bash`. More information with `python update_xml.py -h`. 
    95  
    96 -   **iodef\_atm\_to_oce.xml** and **iodef\_oce\_to_atm.xml**. xml files read by `interpol.f90`. These files are edited by `update_xml.py`. 
    97  
    98 -   **MOSAIX/src/MOSAIX/interpol.f90**. Fortran source using the XIOS library. 
    99  
    100 -   **make_mosaix**. Small script for compiling **interpol.f90**. You need to compile XIOS first. Uses compiler, compiler options and librairies defined by XIOS. 
     152*   **CreateWeightsMask.bash**. Creates the weights to interpolate 
     153    between atmospheric and ocean grid. Weight files are suitable for 
     154    OASIS-MCT. Creates the fraction of ocean in atmospheric grid 
     155    cells. Uses parallel processing with MPI. Has a configuration 
     156    section at the beginning. 
     157 
     158-   **CreateOasisGrids.bash**. Creates files `grids.nc`, `areas.nc` 
     159    and `masks.nc` with information from both models, suitable for 
     160    OASIS-MCT. Mono CPU. Launched by `CreateWeights.bash`. 
     161 
     162-   **update_xml.py**. Python script used to perform on the fly 
     163    editing of xml files. Used by `CreateWeights.bash`. More 
     164    information with `python update_xml.py -h`. 
     165 
     166-   **iodef\_atm\_to_oce.xml** and **iodef\_oce\_to_atm.xml**. xml 
     167    files read by `interpol.f90`. These files are edited by 
     168    `update_xml.py`. 
     169 
     170-   **MOSAIX/src/MOSAIX/interpol.f90**. Fortran source using the XIOS 
     171    library. 
     172 
     173-   **make_mosaix**. Small script for compiling **interpol.f90**. You 
     174    need to compile XIOS first. Uses compiler, compiler options and 
     175    librairies defined by XIOS. 
     176 
    101177```bash  
    102178        make_xios \[options\] 
    … …  
    107183``` 
    108184                         
    109 - **RunoffWeights.py**. Python code to generates weights for run-off. Launched by `CreateWeights.bash`. More information with python `RunOffWeights.py -h`. 
    110  
    111   - Defines a coastal band on land in the atmosphere model. For LMDZ rectilinear grids, the width of this band is parametrized, in grid points. For ico grids, the band contains only point with a fraction of ocean in \]0,1\[. 
    112  
    113   - Defines a coastal band on ocean in the ocean model. The width of this band is parametrized. 
    114  
    115   - An atmosphere point of the coastal band is linked to an ocean point in the coastal band if the distance between the two is less than searchRadius. 
    116  
    117 - **CalvingWeights.py**. Python code to generates weights for calving. Launched by `CreateWeights.bash`. More information with python `CalvingWeights.py -h`. 
    118  
    119   - The atmosphere model integrate the calving over on specific regions. Presently the regions are three latitudinal bands with limits 90°S/40°S/50°N/90°N. 
    120  
    121   - The weights distribute the calving in the ocean in the corresponding latitudinal bands. By default, the distribution is uniform. 
    122  
    123   - For the southernmost band, it is possible to use a geographical distribution read in a file. These files are currently available for eORCA1 and eORCA025. 
     185- **RunoffWeights.py**. Python code to generates weights for 
     186  run-off. Launched by `CreateWeights.bash`. More information with 
     187  python `RunOffWeights.py -h`. 
     188 
     189  - Defines a coastal band on land in the atmosphere model. For LMDZ 
     190    rectilinear grids, the width of this band is parametrized, in grid 
     191    points. For ico grids, the band contains only point with a 
     192    fraction of ocean in \]0,1\[. 
     193 
     194  - Defines a coastal band on ocean in the ocean model. The width of 
     195    this band is parametrized. 
     196 
     197  - An atmosphere point of the coastal band is linked to an ocean 
     198    point in the coastal band if the distance between the two is less 
     199    than searchRadius. 
     200 
     201- **CalvingWeights.py**. Python code to generates weights for 
     202  calving. Launched by `CreateWeights.bash`. More information with 
     203  python `CalvingWeights.py -h`. 
     204 
     205  - The atmosphere model integrate the calving over on specific 
     206    regions. Presently the regions are three latitudinal bands with 
     207    limits 90°S/40°S/50°N/90°N. 
     208 
     209  - The weights distribute the calving in the ocean in the 
     210    corresponding latitudinal bands. By default, the distribution is 
     211    uniform. 
     212 
     213  - For the southernmost band, it is possible to use a geographical 
     214    distribution read in a file. These files are currently available 
     215    for eORCA1 and eORCA025. 
    124216 
    125217## Grids 
    126218 
    127 -   _\[tuv\]orc_. ORCA grid at _\[TUV\]_ points. Masked on land, with area equal to grid box surface. 
    128 -   _tlmd_. LMD rectilinear grid at scalar (physics) point. Masked on land, with area equal to grid box surface. 
    129 -   _tico_. LMD icosahedron grid at scalar (physics) point. Masked on land, with area equal to grid box surface. 
    130 -   _olmd_. LMD rectilinear grid at scalar (physics) point. Not masked, with areas equal to 1. 
    131 -   _oico_. LMD icosahedron grid at scalar (physics) point. Not masked, with areas equal to 1. 
    132 -   _corc_. ORCA grid at _T_ point. Masked on land, duplicated (from periodicity) point masked, with area equal to grid box surface. 
    133 -   _oorc_. ORCA grid at _T_ point. Masked on land, duplicated (from periodicity) point masked, with area equal 1. 
     219  -  _\[tuv\]orc_. ORCA grid at _\[TUV\]_ points. Masked on land, with 
     220   area equal to grid box surface. 
     221  -  _tlmd_. LMD rectilinear grid at scalar (physics) point. Masked on 
     222   land, with area equal to grid box surface. 
     223  -  _tico_. LMD icosahedron grid at scalar (physics) point. Masked on 
     224   land, with area equal to grid box surface. 
     225  -  _olmd_. LMD rectilinear grid at scalar (physics) point. Not 
     226   masked, with areas equal to 1. 
     227  -  _oico_. LMD icosahedron grid at scalar (physics) point. Not 
     228   masked, with areas equal to 1. 
     229  -  _corc_. ORCA grid at _T_ point. Masked on land, duplicated (from 
     230   periodicity) point masked, with area equal to grid box surface. 
     231  -  _oorc_. ORCA grid at _T_ point. Masked on land, duplicated (from 
     232   periodicity) point masked, with area equal 1. 
    134233 
    135234## Generating input grid files 

  • TOOLS/MOSAIX/README.txt

    r3902 r6190  
    2020 
    2121 
    22 XIOS : http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/trunk   revision 1257 
     22Obsolete file : please refer to README.md 

  • TOOLS/MOSAIX/RunoffWeights.py

    r6105 r6190  
    88### =========================================================================== 
    99## 
    10 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    11 ##  software or to read the associated documentation you'll need at least 
    12 ##  one (1) brain in a reasonably working order. Lack of this implement 
    13 ##  will void any warranties (either express or implied). 
    14 ##  O. Marti assumes no responsability for errors, omissions, 
    15 ##  data loss, or any other consequences caused directly or indirectly by 
    16 ##  the usage of his software by incorrectly or partially configured 
    17 ##  personal. 
     10##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     11##  file for an english version of the licence and 
     12##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     13## 
     14##  Permission is hereby granted, free of charge, to any person or 
     15##  organization obtaining a copy of the software and accompanying 
     16##  documentation covered by this license (the "Software") to use, 
     17##  reproduce, display, distribute, execute, and transmit the 
     18##  Software, and to prepare derivative works of the Software, and to 
     19##  permit third-parties to whom the Software is furnished to do so, 
     20##  all subject to the following: 
     21## 
     22##  Warning, to install, configure, run, use any of MOSAIX software or 
     23##  to read the associated documentation you'll need at least one (1) 
     24##  brain in a reasonably working order. Lack of this implement will 
     25##  void any warranties (either express or implied).  Authors assumes 
     26##  no responsability for errors, omissions, data loss, or any other 
     27##  consequences caused directly or indirectly by the usage of his 
     28##  software by incorrectly or partially configured 
     29## 
    1830## 
    1931# SVN information 
    … …  
    5668# Adding arguments 
    5769parser.add_argument ('--oce'          , help='oce model name', type=str, default='eORCA1.2', 
    58                          choices=['ORCA2.3', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 
     70                         choices=['ORCA2.3',  'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 
    5971parser.add_argument ('--atm'          , help='atm model name', type=str, default='LMD9695'    ) 
    6072parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', type=int, default=1 ) 

  • TOOLS/MOSAIX/check_conserv.py

    r4153 r6190  
    55### =========================================================================== 
    66## 
    7 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    8 ##  software or to read the associated documentation you'll need at least 
    9 ##  one (1) brain in a reasonably working order. Lack of this implement 
    10 ##  will void any warranties (either express or implied). 
    11 ##  O. Marti assumes no responsability for errors, omissions, 
    12 ##  data loss, or any other consequences caused directly or indirectly by 
    13 ##  the usage of his software by incorrectly or partially configured 
    14 ##  personal. 
     7##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     8##  file for an english version of the licence and 
     9##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     10## 
     11##  Permission is hereby granted, free of charge, to any person or 
     12##  organization obtaining a copy of the software and accompanying 
     13##  documentation covered by this license (the "Software") to use, 
     14##  reproduce, display, distribute, execute, and transmit the 
     15##  Software, and to prepare derivative works of the Software, and to 
     16##  permit third-parties to whom the Software is furnished to do so, 
     17##  all subject to the following: 
     18## 
     19##  Warning, to install, configure, run, use any of MOSAIX software or 
     20##  to read the associated documentation you'll need at least one (1) 
     21##  brain in a reasonably working order. Lack of this implement will 
     22##  void any warranties (either express or implied).  Authors assumes 
     23##  no responsability for errors, omissions, data loss, or any other 
     24##  consequences caused directly or indirectly by the usage of his 
     25##  software by incorrectly or partially configured 
    1526## 
    1627# 

  • TOOLS/MOSAIX/cotes_etal.py

    r4186 r6190  
    77### =========================================================================== 
    88## 
    9 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    10 ##  software or to read the associated documentation you'll need at least 
    11 ##  one (1) brain in a reasonably working order. Lack of this implement 
    12 ##  will void any warranties (either express or implied). 
    13 ##  O. Marti assumes no responsability for errors, omissions, 
    14 ##  data loss, or any other consequences caused directly or indirectly by 
    15 ##  the usage of his software by incorrectly or partially configured 
    16 ##  personal. 
     9##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     10##  file for an english version of the licence and 
     11##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     12## 
     13##  Permission is hereby granted, free of charge, to any person or 
     14##  organization obtaining a copy of the software and accompanying 
     15##  documentation covered by this license (the "Software") to use, 
     16##  reproduce, display, distribute, execute, and transmit the 
     17##  Software, and to prepare derivative works of the Software, and to 
     18##  permit third-parties to whom the Software is furnished to do so, 
     19##  all subject to the following: 
     20## 
     21##  Warning, to install, configure, run, use any of MOSAIX software or 
     22##  to read the associated documentation you'll need at least one (1) 
     23##  brain in a reasonably working order. Lack of this implement will 
     24##  void any warranties (either express or implied).  Authors assumes 
     25##  no responsability for errors, omissions, data loss, or any other 
     26##  consequences caused directly or indirectly by the usage of his 
     27##  software by incorrectly or partially configured 
    1728## 
    1829# SVN information 

  • TOOLS/MOSAIX/make_mosaix

    r6065 r6190  
    88#  $HeadURL$ 
    99# 
     10##  MOSAIX s under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     11##  file for an english version of the licence and 
     12##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     13## 
     14##  Permission is hereby granted, free of charge, to any person or 
     15##  organization obtaining a copy of the software and accompanying 
     16##  documentation covered by this license (the "Software") to use, 
     17##  reproduce, display, distribute, execute, and transmit the 
     18##  Software, and to prepare derivative works of the Software, and to 
     19##  permit third-parties to whom the Software is furnished to do so, 
     20##  all subject to the following: 
     21## 
     22##  Warning, to install, configure, run, use any of MOSAIX software or 
     23##  to read the associated documentation you'll need at least one (1) 
     24##  brain in a reasonably working order. Lack of this implement will 
     25##  void any warranties (either express or implied).  Authors assumes 
     26##  no responsability for errors, omissions, data loss, or any other 
     27##  consequences caused directly or indirectly by the usage of his 
     28##  software by incorrectly or partially configured 
     29## 
    1030 
    1131# Default values for command line parameters 

  • TOOLS/MOSAIX/nemo.py

    r6105 r6190  
    22## =========================================================================== 
    33## 
    4 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    5 ##  software or to read the associated documentation you'll need at least 
    6 ##  one (1) brain in a reasonably working order. Lack of this implement 
    7 ##  will void any warranties (either express or implied). 
    8 ##  O. Marti assumes no responsability for errors, omissions, 
    9 ##  data loss, or any other consequences caused directly or indirectly by 
    10 ##  the usage of his software by incorrectly or partially configured 
    11 ##  personal. 
     4##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     5##  file for an english version of the licence and 
     6##  "Licence_CeCILL_V2-fr.txt" for a french version. 
    127## 
     8##  Permission is hereby granted, free of charge, to any person or 
     9##  organization obtaining a copy of the software and accompanying 
     10##  documentation covered by this license (the "Software") to use, 
     11##  reproduce, display, distribute, execute, and transmit the 
     12##  Software, and to prepare derivative works of the Software, and to 
     13##  permit third-parties to whom the Software is furnished to do so, 
     14##  all subject to the following: 
     15## 
     16##  Warning, to install, configure, run, use any of MOSAIX software or 
     17##  to read the associated documentation you'll need at least one (1) 
     18##  brain in a reasonably working order. Lack of this implement will 
     19##  void any warranties (either express or implied).  Authors assumes 
     20##  no responsability for errors, omissions, data loss, or any other 
     21##  consequences caused directly or indirectly by the usage of his 
     22##  software by incorrectly or partially configured 
     23## 
     24## =========================================================================== 
    1325''' 
    1426Utilities to plot NEMO ORCA fields 
    … …  
    1729olivier.marti@lsce.ipsl.fr 
    1830''' 
    19  
    20 import sys, numpy as np 
    21 try    : import xarray as xr 
    22 except : pass 
    23  
    24 rpi = np.pi ; rad = rpi / 180.0 
    25  
    26 nperio_valid_range = [0, 1, 4, 4.2, 5, 6, 6.2] 
    2731 
    2832## SVN information 
    … …  
    3337__HeadURL    = "$HeadURL$" 
    3438 
     39import sys, numpy as np 
     40try    : import xarray as xr 
     41except : pass 
     42 
     43rpi = np.pi ; rad = np.deg2rad (1.0) ; dar = np.rad2deg (1.0) 
     44 
     45nperio_valid_range = [0, 1, 4, 4.2, 5, 6, 6.2] 
     46 
     47rday   = 24.*60.*60.     # Day length [s] 
     48rsiyea = 365.25 * rday * 2. * rpi / 6.283076 # Sideral year length [s] 
     49rsiday = rday / (1. + rday / rsiyea) 
     50raamo  =  12.        # Number of months in one year 
     51rjjhh  =  24.        # Number of hours in one day 
     52rhhmm  =  60.        # Number of minutes in one hour 
     53rmmss  =  60.        # Number of seconds in one minute 
     54omega  = 2. * rpi / rsiday # Earth rotation parameter [s-1] 
     55ra     = 6371229.    # Earth radius [m] 
     56grav   = 9.80665     # Gravity [m/s2] 
     57repsi  = np.finfo (1.0).eps 
     58 
     59xList = [ 'x', 'X', 'lon'   , 'longitude' ] 
     60yList = [ 'y', 'Y', 'lat'   , 'latitude'  ] 
     61zList = [ 'z', 'Z', 'depth' , ] 
     62tList = [ 't', 'T', 'time'  , ] 
     63 
     64## SVN information 
     65__Author__   = "$Author$" 
     66__Date__     = "$Date$" 
     67__Revision__ = "$Revision$" 
     68__Id__       = "$Id$" 
     69__HeadURL    = "$HeadURL$" 
     70 
     71def __math__ (tab, default=None) : 
     72    math = default 
     73    try    : 
     74        if type (tab) == xr.core.dataarray.DataArray : math = xr 
     75    except : 
     76        pass 
     77 
     78    try    : 
     79        if type (tab) == np.ndarray : math = np 
     80    except : 
     81        pass 
     82             
     83    return math 
     84 
    3585def __guessNperio__ (jpj, jpi, nperio=None) : 
    3686    ''' 
    … …  
    4494    if nperio == None : 
    4595        ## Values for NEMO version < 4.2 
    46         if jpi ==  182 : nperio = 4   #   ORCA2. We choose legacy orca2. 
    47         if jpi ==  362 : nperio = 6   #   eORCA1. 
    48         if jpi == 1442 : nperio = 6   #   ORCA025. 
    49  
    50         if jpj ==  149 : nperio = 4   #   ORCA2. We choose legacy orca2. 
    51         if jpj ==  332 : nperio = 6   #   eORCA1. 
    52  
     96        if jpi ==  182 : 
     97            nperio = 4   # ORCA2. We choose legacy orca2. 
     98            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'T' 
     99        if jpi ==  362 : # eORCA1. 
     100            nperio = 6   
     101            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     102        if jpi == 1442 :  # ORCA025. 
     103            nperio = 6  
     104            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     105        if jpj ==  149 : # ORCA2. We choose legacy orca2. 
     106            nperio = 4 
     107            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     108        if jpj ==  294 : # ORCA1 
     109            nperio = 6 
     110            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     111        if jpj ==  332 : # eORCA1. 
     112            nperio = 6   
     113            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     114             
    53115        ## Values for NEMO version >= 4.2. No more halo points 
    54         if jpi ==  180 : nperio = 4.2 #   ORCA2. We choose legacy orca2. 
    55         if jpi ==  360 : nperio = 0   #   eORCA1. 
    56         if jpi == 1440 : nperio = 0   #   ORCA025. 
    57  
    58         if jpj ==  148 : nperio = 4.2 #   ORCA2. We choose legacy orca2. 
    59         if jpj ==  330 : nperio = 0   #   eORCA1. 
    60  
    61         # 
     116        if jpi ==  180 : 
     117            nperio = 4.2 # ORCA2. We choose legacy orca2. 
     118            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     119        if jpi ==  360 : # eORCA1. 
     120            nperio = 6.2 
     121            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     122        if jpi == 1440 : # ORCA025. 
     123            nperio = 6.2 
     124            Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 
     125 
     126        if jpj ==  148 : nperio = 4.2 # ORCA2. We choose legacy orca2. 
     127        if jpj ==  330 : nperio = 6.2 # eORCA1. 
     128 
    62129        if nperio == None : 
    63             sys.exit ('in nemo module : nperio not found, and cannot by guessed' ) 
     130            sys.exit ('in nemo module : nperio not found, and cannot by guessed') 
    64131        else : 
    65             sys.ext  ('nperio set as {:} (deduced from jpi={:d}) : nemo.py is not ready for this value'.format (nperio, jpi) ) 
     132            if nperio in nperio_valid_range : 
     133                print ('nperio set as {:} (deduced from jpi={:d})'.format (nperio, jpi)) 
     134            else :  
     135                sys.exit  ('nperio set as {:} (deduced from jpi={:d}) : nemo.py is not ready for this value'.format (nperio, jpi)) 
    66136             
    67             print ('nperio set as {:} (deduced from jpi={:d})'.format (nperio, jpi)) 
    68  
    69137    return nperio 
    70138 
    … …  
    76144 
    77145    Inputs 
    78     ptab : xarray array 
     146         ptab : xarray array 
    79147 
    80148    Credits : who is the original author ? 
    81149    ''' 
    82150    gP = None 
    83     if isinstance (ptab, xr.core.dataarray.DataArray) : 
     151    math = __math__ (ptab) 
     152    if math == xr : 
    84153        if 'x_c' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'T' 
    85154        if 'x_f' in ptab.dims and 'y_c' in ptab.dims                        : gP = 'U' 
    … …  
    93162              
    94163        if gP == None : 
    95             sys.exit ('in nemo module : cd_type not found, and cannot by guessed' ) 
     164            sys.exit ('in nemo module : cd_type not found, and cannot by guessed') 
    96165        else : 
    97166            print ('Grid set as', gP, 'deduced from dims ', ptab.dims) 
    98167            return gP 
    99168    else : 
    100          sys.exit ('in nemo module : cd_type not found, input is not an xarray data' ) 
     169         sys.exit ('in nemo module : cd_type not found, input is not an xarray data') 
    101170         #return gP 
    102171 
    103 def fixed_lon (lon) : 
     172def __findAxis__ (tab, axis='z') : 
     173    ''' 
     174    Find number and name of the requested axis 
     175    ''' 
     176    math = __math__ (tab) 
     177    ix = None ; ax = None 
     178 
     179    if axis in xList : 
     180        axList = [ 'x', 'X', 
     181                   'lon', 'nav_lon', 'nav_lon_T', 'nav_lon_U', 'nav_lon_V', 'nav_lon_F', 'nav_lon_W', 
     182                   'x_grid_T', 'x_grid_U', 'x_grid_V', 'x_grid_F', 'x_grid_W', 
     183                   'glam', 'glamt', 'glamu', 'glamv', 'glamf', 'glamw' ] 
     184        unList = [ 'degrees_east' ] 
     185    if axis in yList : 
     186        axList = [ 'y', 'Y', 'lat', 
     187                   'nav_lat', 'nav_lat_T', 'nav_lat_U', 'nav_lat_V', 'nav_lat_F', 'nav_lat_W', 
     188                   'y_grid_T', 'y_grid_U', 'y_grid_V', 'y_grid_F', 'y_grid_W', 
     189                   'gphi', 'gphi', 'gphiu', 'gphiv', 'gphif', 'gphiw'] 
     190        unList = [ 'degrees_north' ] 
     191    if axis in zList : 
     192        axList = [ 'z', 'Z', 
     193                   'depth', 'deptht', 'depthu', 'depthv', 'depthf', 'depthw', 
     194                   'olevel' ] 
     195        unList = [ 'm', 'meter' ] 
     196    if axis in tList : 
     197        axList = [ 't', 'T', 'time', 'time_counter' ] 
     198        unList = [ 'second', 'minute', 'hour', 'day', 'month' ] 
     199     
     200    if math == xr : 
     201        for Name in axList : 
     202            try    : 
     203                ix = tab.dims.index (Name) 
     204                ax = Name 
     205            except : pass 
     206 
     207        for i, dim in enumerate (tab.dims) : 
     208            if 'units' in tab.coords[dim].attrs.keys() : 
     209                for name in unList : 
     210                    if name in tab.coords[dim].attrs['units'] : 
     211                        ix = i 
     212                        ax = dim 
     213    else : 
     214        if axis in xList : ix=-1 
     215        if axis in yList : 
     216            if len(tab.shape) >= 2 : ix=-2 
     217        if axis in zList : 
     218            if len(tab.shape) >= 3 : ix=-3 
     219        if axis in tList : 
     220            if len(tab.shape) >=3  : ix=-3 
     221            if len(tab.shape) >=4  : ix=-4 
     222        
     223    return ix, ax 
     224                 
     225def fixed_lon (lon, center_lon=0.0) : 
    104226    ''' 
    105227    Returns corrected longitudes for nicer plots 
    106228 
    107     fixed_lon (lon) 
    108     lon : longitudes of the grid. At least 2D. 
     229    lon        : longitudes of the grid. At least 2D. 
     230    center_lon : center longitude. Default=0. 
    109231 
    110232    Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7 
    111233    ''' 
     234    math = __math__ (lon) 
     235     
    112236    fixed_lon = lon.copy () 
    113     for i, start in enumerate (np.argmax (np.abs (np.diff (lon, axis=-1)) > 180, axis=-1)) : 
    114         fixed_lon[..., i, start+1:] += 360 
     237         
     238    fixed_lon = math.where (fixed_lon > center_lon+180., fixed_lon-360.0, fixed_lon) 
     239    fixed_lon = math.where (fixed_lon < center_lon-180., fixed_lon+360.0, fixed_lon) 
     240     
     241    for i, start in enumerate (np.argmax (np.abs (np.diff (fixed_lon, axis=-1)) > 180., axis=-1)) : 
     242        fixed_lon [..., i, start+1:] += 360 
    115243 
    116244    # Special case for eORCA025 
    117     if lon.shape[-1] == 1442 : lon [..., -2, :] =  lon [..., -3, :] 
    118     if lon.shape[-1] == 1440 : lon [..., -1, :] =  lon [..., -2, :] 
     245    if fixed_lon.shape [-1] == 1442 : fixed_lon [..., -2, :] = fixed_lon [..., -3, :] 
     246    if fixed_lon.shape [-1] == 1440 : fixed_lon [..., -1, :] = fixed_lon [..., -2, :] 
     247 
     248    if fixed_lon.min () > center_lon : fixed_lon = fixed_lon-360.0 
    119249                 
    120250    return fixed_lon 
    … …  
    126256    lat : latitudes of the grid. At least 2D. 
    127257    ''' 
    128     jeq = np.argmin (np.abs (np.float64(lat[:,0]))) 
     258    math = __math__ (lat) 
     259    ix, ax = __findAxis__ (lat, 'x') 
     260    iy, ay = __findAxis__ (lat, 'y') 
     261 
     262    if math == xr : 
     263        jeq = int ( np.mean ( np.argmin (np.abs (np.float64 (lat)), axis=iy) ) ) 
     264    else :  
     265        jeq = np.argmin (np.abs (np.float64 (lat[...,:, 0]))) 
    129266    return jeq 
    130267 
    … …  
    137274    ''' 
    138275    if np.max (lat) != None : 
    139         je     = jeq (lat) 
    140         lon1D = lon[..., je, :] 
     276        je    = jeq (lat) 
     277        lon1D = lon [..., je, :] 
    141278    else : 
    142         jpj, jpi = lon.shape[-2:] 
    143         lon1D = lon[..., jpj//3, :] 
     279        jpj, jpi = lon.shape [-2:] 
     280        lon1D    = lon [..., jpj//3, :] 
    144281 
    145282    start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1) 
    146     lon1D[..., start+1:] += 360 
     283    lon1D [..., start+1:] += 360 
    147284         
    148285    return lon1D 
    … …  
    155292    diff [optional] : tolerance 
    156293    ''' 
     294    math = __math__ (lat) 
    157295    if diff == None : 
    158         dy = np.float64 (np.mean (np.abs (lat - np.roll(lat,shift=1,axis=-2)))) 
     296        dy   = np.float64 (np.mean (np.abs (lat - np.roll (lat,shift=1,axis=-2, roll_coords=False)))) 
    159297        diff = dy/100. 
    160298     
    … …  
    172310    lat : latitudes of the grid (2D) 
    173311    ''' 
     312    math = __math__ (lat) 
    174313    jpj, jpi = lat.shape[-2:] 
    175     try : 
    176         if type (lat) == xr.core.dataarray.DataArray : math = xr 
    177         else : math = np 
    178     except : math = np 
    179314 
    180315    dy     = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2)))) 
    … …  
    185320    lat_ave = np.mean (lat, axis=-1) 
    186321 
    187     if ( np.abs (lat_eq) < dy/100. ) : # T, U or W grid 
     322    if (np.abs (lat_eq) < dy/100.) : # T, U or W grid 
    188323        dys    = (90.-lat_reg) / (jpj-jreg-1)*0.5 
    189324        yrange = 90.-dys-lat_reg 
    190     else                             :  # V or F grid 
     325    else                           :  # V or F grid 
    191326        yrange = 90.    -lat_reg 
    192327         
    193     lat1D = math.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) )    
     328    lat1D = math.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1))    
    194329         
    195330    if math == xr : lat1D.attrs = lat.attrs 
    … …  
    206341 
    207342def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) : 
     343    math = __math__ (ptab) 
    208344    try : 
    209         lon = lon.to_masked_array() 
    210         lat = lat.to_masked_array() 
     345        lon = lon.copy().to_masked_array() 
     346        lat = lat.copy().to_masked_array() 
    211347    except : pass 
    212348             
    213     mask = np.logical_and ( np.logical_and(lat>y0, lat<y1),  
    214         np.logical_or (np.logical_or (np.logical_and(lon>x0, lon<x1), np.logical_and(lon+360>x0, lon+360<x1)), 
    215                                          np.logical_and(lon-360>x0, lon-360<x1)) ) 
    216     try    : tab = xr.where (mask, ptab, np.nan) 
    217     except : tab = np.where (mask, ptab, np.nan) 
     349    mask = np.logical_and (np.logical_and(lat>y0, lat<y1),  
     350            np.logical_or (np.logical_or (np.logical_and(lon>x0, lon<x1), np.logical_and(lon+360>x0, lon+360<x1)), 
     351                                      np.logical_and(lon-360>x0, lon-360<x1))) 
     352    tab = math.where (mask, ptab, np.nan) 
    218353     
    219354    return tab 
    … …  
    231366     
    232367    ''' 
    233     if tab.shape[-1] == 1 : 
    234         extend = tab 
     368    math = __math__ (tab) 
     369     
     370    if tab.shape[-1] == 1 : extend = tab 
    235371 
    236372    else : 
    237373        if jpi == None : jpi = tab.shape[-1] 
    238      
    239         try :  
    240             if type (tab) ==  xr.core.dataarray.DataArray : math = xr 
    241             else : math = np 
    242         except : math = np 
    243374 
    244375        if Lon : xplus = -360.0 
    … …  
    248379            extend = tab 
    249380        else : 
     381            if nperio == 0 or nperio == 4.2 : 
     382                istart = 0 ; le=jpi+1 ; la=0 
    250383            if nperio == 1 : 
    251384                istart = 0 ; le=jpi+1 ; la=0 
    252             if nperio > 1 : # OPA case with to halo points for periodicity 
     385            if nperio == 4 or nperio == 6 : # OPA case with two halo points for periodicity 
    253386                istart = 1 ; le=jpi-2 ; la=1  # Perfect, except at the pole that should be masked by lbc_plot 
    254                  
    255             if math == xr :  
    256                 extend = xr.concat      ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), dim=tab.dims[-1] ) 
    257             if math == np : 
    258                 extend = np.concatenate ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), axis=-1 ) 
    259              
     387            
     388            if math == xr : 
     389                extend = np.concatenate ((tab.values[..., istart   :istart+le+1    ] + xplus, 
     390                                          tab.values[..., istart+la:istart+la+jplus]         ), axis=-1) 
     391                lon    = tab.dims[-1] 
     392                new_coords = [] 
     393                for coord in tab.dims : 
     394                    if coord == lon : new_coords.append ( np.arange( extend.shape[-1])) 
     395                    else            : new_coords.append ( tab.coords[coord].values) 
     396                extend = xr.DataArray ( extend, dims=tab.dims, coords=new_coords ) 
     397            else :  
     398                extend = np.concatenate ((tab [..., istart   :istart+le+1    ] + xplus, 
     399                                          tab [..., istart+la:istart+la+jplus]          ), axis=-1) 
    260400    return extend 
    261401 
    … …  
    286426                 'time_counter', 'time', 'tbnds',  
    287427                 'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] : 
    288         if dim in ff.dims : coord_order.insert (0, dim ) 
     428        if dim in ff.dims : coord_order.insert (0, dim) 
    289429         
    290430    ff = ff.transpose (*coord_order) 
    291431    return ff 
    292432 
    293 def lbc_init (ptab, nperio=None, cd_type='T') : 
    294     """ 
     433def lbc_init (ptab, nperio=None) : 
     434    ''' 
    295435    Prepare for all lbc calls 
    296436     
    … …  
    305445 
    306446    See NEMO documentation for further details 
    307     """ 
     447    ''' 
    308448    jpj, jpi = ptab.shape[-2:] 
    309449    if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio) 
    310450     
    311451    if nperio not in nperio_valid_range : 
    312         print ( 'nperio=', nperio, ' is not in the valid range', nperio_valid_range ) 
    313         sys.exit() 
    314  
    315     if cd_type == None and isinstance (ptab, xr.core.dataarray.DataArray) : 
    316         cd_type = __guessPoint__ (ptab) 
    317  
    318     if cd_type not in ['T', 'U', 'V', 'F', 'W', 'F' ]: 
    319         print ( 'cd_type=', cd_type, ' is not in the list T, U, V, F, W, F' ) 
     452        print ('nperio=', nperio, ' is not in the valid range', nperio_valid_range) 
    320453        sys.exit () 
    321454 
    322     return jpj, jpi, nperio, cd_type 
     455    return jpj, jpi, nperio 
    323456         
    324457         
    325458def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) : 
    326     """ 
     459    ''' 
    327460    Set periodicity on input field 
    328     ptab      : Input array (works  
    329       rank 2 at least : ptab[...., lat, lon] 
     461    ptab      : Input array (works for rank 2 at least : ptab[...., lat, lon]) 
    330462    nperio    : Type of periodicity 
    331463    cd_type   : Grid specification : T, U, V or F 
    … …  
    333465     
    334466    See NEMO documentation for further details 
    335     """ 
    336     jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) 
     467    ''' 
     468    jpj, jpi, nperio = lbc_init (ptab, nperio) 
    337469    psgn   = ptab.dtype.type (psgn) 
    338  
     470    math = __math__ (ptab) 
     471     
     472    if math == xr : ztab = ptab.values.copy () 
     473    else          : ztab = ptab.copy () 
    339474    # 
    340475    #> East-West boundary conditions 
    … …  
    342477    if nperio in [1, 4, 6] : 
    343478        # ... cyclic 
    344         ptab [..., :,  0] = ptab [..., :, -2] 
    345         ptab [..., :, -1] = ptab [..., :,  1] 
     479        ztab [..., :,  0] = ztab [..., :, -2] 
     480        ztab [..., :, -1] = ztab [..., :,  1] 
    346481    # 
    347482    #> North-South boundary conditions 
    … …  
    349484    if nperio in [3, 4] :  # North fold T-point pivot 
    350485        if cd_type in [ 'T', 'W' ] : # T-, W-point 
    351             ptab [..., -1, 1:       ] = psgn * ptab [..., -3, -1:0:-1] 
    352             ptab [..., -1, 0        ] = psgn * ptab [..., -3, 2            ] 
    353             ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2:0:-1  ] 
     486            ztab [..., -1, 1:       ] = psgn * ztab [..., -3, -1:0:-1      ] 
     487            ztab [..., -1, 0        ] = psgn * ztab [..., -3, 2            ] 
     488            ztab [..., -2, jpi//2:  ] = psgn * ztab [..., -2, jpi//2:0:-1  ] 
    354489                 
    355490        if cd_type == 'U' : 
    356             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ]        
    357             ptab [..., -1,  0       ] = psgn * ptab [..., -3,  1           ] 
    358             ptab [..., -1, -1       ] = psgn * ptab [..., -3, -2           ] 
     491            ztab [..., -1, 0:-1     ] = psgn * ztab [..., -3, -1:0:-1      ]        
     492            ztab [..., -1,  0       ] = psgn * ztab [..., -3,  1           ] 
     493            ztab [..., -1, -1       ] = psgn * ztab [..., -3, -2           ] 
    359494                 
    360495            if nemo_4U_bug : 
    361                 ptab [..., -2, jpi//2+1:-1] = psgn * ptab [..., -2, jpi//2-2:0:-1] 
    362                 ptab [..., -2, jpi//2-1   ] = psgn * ptab [..., -2, jpi//2     ] 
     496                ztab [..., -2, jpi//2+1:-1] = psgn * ztab [..., -2, jpi//2-2:0:-1] 
     497                ztab [..., -2, jpi//2-1   ] = psgn * ztab [..., -2, jpi//2       ] 
    363498            else : 
    364                 ptab [..., -2, jpi//2-1:-1] = psgn * ptab [..., -2, jpi//2:0:-1] 
     499                ztab [..., -2, jpi//2-1:-1] = psgn * ztab [..., -2, jpi//2:0:-1] 
    365500                 
    366501        if cd_type == 'V' :  
    367             ptab [..., -2, 1:       ] = psgn * ptab [..., -3, jpi-1:0:-1   ] 
    368             ptab [..., -1, 1:       ] = psgn * ptab [..., -4, -1:0:-1      ]    
    369             ptab [..., -1, 0        ] = psgn * ptab [..., -4, 2            ] 
     502            ztab [..., -2, 1:       ] = psgn * ztab [..., -3, jpi-1:0:-1   ] 
     503            ztab [..., -1, 1:       ] = psgn * ztab [..., -4, -1:0:-1      ]    
     504            ztab [..., -1, 0        ] = psgn * ztab [..., -4, 2            ] 
    370505                 
    371506        if cd_type == 'F' : 
    372             ptab [..., -2, 0:-1     ] = psgn * ptab [..., -3, -1:0:-1      ] 
    373             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -4, -1:0:-1      ] 
    374             ptab [..., -1,  0       ] = psgn * ptab [..., -4,  1           ] 
    375             ptab [..., -1, -1       ] = psgn * ptab [..., -4, -2           ] 
     507            ztab [..., -2, 0:-1     ] = psgn * ztab [..., -3, -1:0:-1      ] 
     508            ztab [..., -1, 0:-1     ] = psgn * ztab [..., -4, -1:0:-1      ] 
     509            ztab [..., -1,  0       ] = psgn * ztab [..., -4,  1           ] 
     510            ztab [..., -1, -1       ] = psgn * ztab [..., -4, -2           ] 
    376511 
    377512    if nperio in [4.2] :  # North fold T-point pivot 
    378513        if cd_type in [ 'T', 'W' ] : # T-, W-point 
    379             ptab [..., -1, jpi//2:  ] = psgn * ptab [..., -1, jpi//2:0:-1  ] 
     514            ztab [..., -1, jpi//2:  ] = psgn * ztab [..., -1, jpi//2:0:-1  ] 
    380515                 
    381516        if cd_type == 'U' : 
    382             ptab [..., -1, jpi//2-1:-1] = psgn * ptab [..., -1, jpi//2:0:-1] 
     517            ztab [..., -1, jpi//2-1:-1] = psgn * ztab [..., -1, jpi//2:0:-1] 
    383518                 
    384519        if cd_type == 'V' :  
    385             ptab [..., -1, 1:       ] = psgn * ptab [..., -2, jpi-1:0:-1   ] 
     520            ztab [..., -1, 1:       ] = psgn * ztab [..., -2, jpi-1:0:-1   ] 
    386521                 
    387522        if cd_type == 'F' : 
    388             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -2, -1:0:-1      ] 
    389           
    390          
     523            ztab [..., -1, 0:-1     ] = psgn * ztab [..., -2, -1:0:-1      ] 
     524 
    391525    if nperio in [5, 6] :            #  North fold F-point pivot   
    392526        if cd_type in ['T', 'W']  : 
    393             ptab [..., -1, 0:       ] = psgn * ptab [..., -2, -1::-1       ] 
     527            ztab [..., -1, 0:       ] = psgn * ztab [..., -2, -1::-1       ] 
    394528                 
    395529        if cd_type == 'U' : 
    396             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -2, -2::-1       ]        
    397             ptab [..., -1, -1       ] = psgn * ptab [..., -2, 0            ] # Bug ? 
     530            ztab [..., -1, 0:-1     ] = psgn * ztab [..., -2, -2::-1       ]        
     531            ztab [..., -1, -1       ] = psgn * ztab [..., -2, 0            ] # Bug ? 
    398532                 
    399533        if cd_type == 'V' : 
    400             ptab [..., -1, 0:       ] = psgn * ptab [..., -3, -1::-1       ] 
    401             ptab [..., -2, jpi//2:  ] = psgn * ptab [..., -2, jpi//2-1::-1 ] 
     534            ztab [..., -1, 0:       ] = psgn * ztab [..., -3, -1::-1       ] 
     535            ztab [..., -2, jpi//2:  ] = psgn * ztab [..., -2, jpi//2-1::-1 ] 
    402536                 
    403537        if cd_type == 'F' : 
    404             ptab [..., -1, 0:-1     ] = psgn * ptab [..., -3, -2::-1       ] 
    405             ptab [..., -1, -1       ] = psgn * ptab [..., -3, 0            ] 
    406             ptab [..., -2, jpi//2:-1] = psgn * ptab [..., -2, jpi//2-2::-1 ] 
    407                  
    408     return ptab 
    409  
    410 def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) : 
    411     # 
    412     """ 
    413     Mask fields on duplicated points 
    414     ptab      : Input array 
    415       rank 2 at least : ptab[...., lat, lon] 
    416     nperio    : Type of periodicity 
    417     cd_type   : Grid specification : T, U, V or F 
    418      
    419     See NEMO documentation for further details 
    420     """ 
    421  
    422     jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) 
    423           
     538            ztab [..., -1, 0:-1     ] = psgn * ztab [..., -3, -2::-1       ] 
     539            ztab [..., -1, -1       ] = psgn * ztab [..., -3, 0            ] 
     540            ztab [..., -2, jpi//2:-1] = psgn * ztab [..., -2, jpi//2-2::-1 ] 
     541 
    424542    # 
    425543    #> East-West boundary conditions 
    … …  
    427545    if nperio in [1, 4, 6] : 
    428546        # ... cyclic 
    429         ptab [...,:,  0] = sval 
    430         ptab [...,:, -1] = sval 
     547        ztab [..., :,  0] = ztab [..., :, -2] 
     548        ztab [..., :, -1] = ztab [..., :,  1] 
     549 
     550    if math == xr : 
     551        ztab = xr.DataArray ( ztab, dims=ptab.dims, coords=ptab.coords) 
     552         
     553    return ztab 
     554 
     555def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) : 
     556    # 
     557    ''' 
     558    Mask fields on duplicated points 
     559    ptab      : Input array. Rank 2 at least : ptab [...., lat, lon] 
     560    nperio    : Type of periodicity 
     561    cd_type   : Grid specification : T, U, V or F 
     562     
     563    See NEMO documentation for further details 
     564    ''' 
     565    jpj, jpi, nperio = lbc_init (ptab, nperio) 
     566    ztab = ptab.copy () 
     567 
     568    # 
     569    #> East-West boundary conditions 
     570    # ------------------------------ 
     571    if nperio in [1, 4, 6] : 
     572        # ... cyclic 
     573        ztab [..., :,  0] = sval 
     574        ztab [..., :, -1] = sval 
    431575 
    432576    # 
    … …  
    434578    # -------------------------------- 
    435579    if nperio in [1, 3, 4, 5, 6] : 
    436         ptab [...,0,:] = sval 
     580        ztab [..., 0, :] = sval 
    437581         
    438582    # 
    … …  
    441585    if nperio in [3, 4] :  # North fold T-point pivot 
    442586        if cd_type in [ 'T', 'W' ] : # T-, W-point 
    443             ptab [..., -1, 1:       ] = sval 
    444             ptab [..., -1, 0        ] = sval  
    445             ptab [..., -2, jpi//2:  ] = sval 
     587            ztab [..., -1,  :         ] = sval 
     588            ztab [..., -2, :jpi//2  ] = sval 
    446589                 
    447590        if cd_type == 'U' : 
    448             ptab [..., -1,  :       ] = sval   
    449             ptab [..., -2, jpi//2:  ] = sval 
     591            ztab [..., -1,  :         ] = sval   
     592            ztab [..., -2, jpi//2+1:  ] = sval 
    450593                 
    451594        if cd_type == 'V' : 
    452             ptab [..., -2, :       ] = sval 
    453             ptab [..., -1, :       ] = sval    
     595            ztab [..., -2, :       ] = sval 
     596            ztab [..., -1, :       ] = sval    
    454597                 
    455598        if cd_type == 'F' : 
    456             ptab [..., -2, :     ] = sval 
    457             ptab [..., -1, :     ] = sval 
     599            ztab [..., -2, :       ] = sval 
     600            ztab [..., -1, :       ] = sval 
    458601 
    459602    if nperio in [4.2] :  # North fold T-point pivot 
    460603        if cd_type in [ 'T', 'W' ] : # T-, W-point 
    461             ptab [..., -1, jpi//2:  ] = sval 
     604            ztab [..., -1, jpi//2  :  ] = sval 
    462605                 
    463606        if cd_type == 'U' : 
    464             ptab [..., -1, jpi//2-1:-1] = sval 
     607            ztab [..., -1, jpi//2-1:-1] = sval 
    465608                 
    466609        if cd_type == 'V' :  
    467             ptab [..., -1, 1:       ] = sval 
     610            ztab [..., -1, 1:       ] = sval 
    468611                 
    469612        if cd_type == 'F' : 
    470             ptab [..., -1, 0:-1     ] = sval 
     613            ztab [..., -1, 0:-1     ] = sval 
    471614     
    472615    if nperio in [5, 6] :            #  North fold F-point pivot 
    473616        if cd_type in ['T', 'W']  : 
    474             ptab [..., -1, 0:       ] = sval 
     617            ztab [..., -1, 0:       ] = sval 
    475618                 
    476619        if cd_type == 'U' : 
    477             ptab [..., -1, 0:-1     ] = sval        
    478             ptab [..., -1, -1       ] = sval 
     620            ztab [..., -1, 0:-1     ] = sval        
     621            ztab [..., -1, -1       ] = sval 
    479622              
    480623        if cd_type == 'V' : 
    481             ptab [..., -1, 0:       ] = sval 
    482             ptab [..., -2, jpi//2:  ] = sval 
     624            ztab [..., -1, 0:       ] = sval 
     625            ztab [..., -2, jpi//2:  ] = sval 
    483626                              
    484627        if cd_type == 'F' : 
    485             ptab [..., -1, 0:-1     ] = sval 
    486             ptab [..., -1, -1       ] = sval 
    487             ptab [..., -2, jpi//2:-1] = sval 
    488  
    489     return ptab 
    490  
     628            ztab [..., -1, 0:-1       ] = sval 
     629            ztab [..., -1, -1         ] = sval 
     630            ztab [..., -2, jpi//2+1:-1] = sval 
     631 
     632    return ztab 
    491633 
    492634def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) : 
    493     """ 
     635    ''' 
    494636    Set periodicity on input field, adapted for plotting for any cartopy projection 
    495     ptab      : Input array (works  
    496       rank 2 at least : ptab[...., lat, lon] 
     637    ptab      : Input array. Rank 2 at least : ptab[...., lat, lon] 
    497638    nperio    : Type of periodicity 
    498639    cd_type   : Grid specification : T, U, V or F 
    … …  
    500641     
    501642    See NEMO documentation for further details 
    502     """ 
    503  
    504     jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) 
     643    ''' 
     644 
     645    jpj, jpi, nperio = lbc_init (ptab, nperio) 
    505646    psgn   = ptab.dtype.type (psgn) 
    506      
     647    ztab   = ptab.copy () 
    507648    # 
    508649    #> East-West boundary conditions 
    … …  
    510651    if nperio in [1, 4, 6] : 
    511652        # ... cyclic 
    512         ptab [..., :,  0] = ptab [..., :, -2] 
    513         ptab [..., :, -1] = ptab [..., :,  1] 
     653        ztab [..., :,  0] = ztab [..., :, -2] 
     654        ztab [..., :, -1] = ztab [..., :,  1] 
     655 
     656    #> Masks south 
     657    # ------------ 
     658    if nperio in [4, 6] : ztab [..., 0, : ] = sval 
    514659         
    515660    # 
    … …  
    518663    if nperio in [3, 4] :  # North fold T-point pivot 
    519664        if cd_type in [ 'T', 'W' ] : # T-, W-point 
    520             ptab [..., -1, :] = sval 
    521             ptab [..., -2, :jpi//2] = sval 
    522             #ptab [..., -2, :] = sval 
    523                        
     665            ztab [..., -1,  :      ] = sval 
     666            #ztab [..., -2, jpi//2: ] = sval 
     667            ztab [..., -2, :jpi//2 ] = sval # Give better plots than above 
    524668        if cd_type == 'U' : 
    525             ptab [..., -1, : ] = sval 
     669            ztab [..., -1, : ] = sval 
    526670 
    527671        if cd_type == 'V' :  
    528             ptab [..., -2, : ] = sval 
    529             ptab [..., -1, : ] = sval 
     672            ztab [..., -2, : ] = sval 
     673            ztab [..., -1, : ] = sval 
    530674             
    531675        if cd_type == 'F' : 
    532             ptab [..., -2, : ] = sval 
    533             ptab [..., -1, : ] = sval 
     676            ztab [..., -2, : ] = sval 
     677            ztab [..., -1, : ] = sval 
    534678 
    535679    if nperio in [4.2] :  # North fold T-point pivot 
    536680        if cd_type in [ 'T', 'W' ] : # T-, W-point 
    537             ptab [..., -1, jpi//2:  ] = sval 
     681            ztab [..., -1, jpi//2:  ] = sval 
    538682                 
    539683        if cd_type == 'U' : 
    540             ptab [..., -1, jpi//2-1:-1] = sval 
     684            ztab [..., -1, jpi//2-1:-1] = sval 
    541685                 
    542686        if cd_type == 'V' :  
    543             ptab [..., -1, 1:       ] = sval 
     687            ztab [..., -1, 1:       ] = sval 
    544688                 
    545689        if cd_type == 'F' : 
    546             ptab [..., -1, 0:-1     ] = sval 
     690            ztab [..., -1, 0:-1     ] = sval 
    547691       
    548692    if nperio in [5, 6] :            #  North fold F-point pivot   
    549693        if cd_type in ['T', 'W']  : 
    550             ptab [..., -1, : ] = sval 
     694            ztab [..., -1, : ] = sval 
    551695                 
    552696        if cd_type == 'U' : 
    553             ptab [..., -1, : ] = sval       
     697            ztab [..., -1, : ] = sval       
    554698              
    555699        if cd_type == 'V' : 
    556             ptab [..., -1, : ] = sval 
    557             ptab [..., -2, jpi//2:  ] = sval 
     700            ztab [..., -1, :        ] = sval 
     701            ztab [..., -2, jpi//2:  ] = sval 
    558702                              
    559703        if cd_type == 'F' : 
    560             ptab [..., -1, : ] = sval 
    561             ptab [..., -2, jpi//2:-1] = sval 
    562  
     704            ztab [..., -1, :          ] = sval 
     705            ztab [..., -2, jpi//2+1:-1] = sval 
     706 
     707    return ztab 
     708 
     709def lbc_add (ptab, nperio=None, cd_type =None) : 
     710    ''' 
     711    Handle NEMO domain changes between NEMO 4.0 to NEMO 4.2 
     712      Peridodicity halo has been removes 
     713    This routine adds the halos if needed 
     714 
     715    ptab      : Input array (works  
     716      rank 2 at least : ptab[...., lat, lon] 
     717    nperio    : Type of periodicity 
     718  
     719    See NEMO documentation for further details 
     720    ''' 
     721 
     722    jpj, jpi, nperio = lbc_init (ptab, nperio) 
     723 
     724    t_shape = np.array (ptab.shape) 
     725 
     726    if nperio == 4.2 or nperio == 6.2 : 
     727        math = __math__ (ptab)  
     728       
     729        ext_shape = t_shape 
     730        ext_shape[-1] = ext_shape[-1] + 2 
     731        ext_shape[-2] = ext_shape[-2] + 1 
     732 
     733        if math == xr : ptab_ext = xr.DataArray (np.empty (ext_shape), dims=ptab.dims)  
     734        else          : ptab_ext =               np.empty (ext_shape) 
     735             
     736        ptab_ext[..., :-1, 1:-1] = ptab 
     737         
     738        if nperio == 4.2 : ptab_ext = lbc (ptab_ext, 4, cd_type) 
     739        if nperio == 6.2 : ptab_ext = lbc (ptab_ext, 6, cd_type) 
     740 
     741        if math == xr : 
     742            for attr in ptab.attrs : 
     743                ptab_ext.attrs [attr] = ptab.attrs [attr] 
     744 
     745    else : ptab_ext = ptab 
     746         
    563747    return ptab 
    564748 
    565 def geo2oce (pxx, pyy, pzz, glam, gphi) :  
     749def lbc_del (ptab, nperio=None) : 
     750    ''' 
     751    Handle NEMO domain changes between NEMO 4.0 to NEMO 4.2 
     752      Peridodicity halo has been removes 
     753    This routine removes the halos if needed 
     754 
     755    ptab      : Input array (works  
     756      rank 2 at least : ptab[...., lat, lon] 
     757    nperio    : Type of periodicity 
     758  
     759    See NEMO documentation for further details 
     760    ''' 
     761 
     762    jpj, jpi, nperio = lbc_init (ptab, nperio) 
     763 
     764    if nperio == 4 or nperio == 6 : 
     765        return ptab[..., :-1, 1:-1] 
     766    else : 
     767        return ptab 
     768 
     769def lbc_index (jj, ii, jpj, jpi, nperio=None, cd_type='T') : 
     770    ''' 
     771    For indexes of a NEMO point, give the corresponding point inside the util domain 
     772    jj, ii    : indexes 
     773    jpi, jpi  : size of domain 
     774    nperio    : type of periodicity 
     775    cd_type   : grid specification : T, U, V or F 
     776     
     777    See NEMO documentation for further details 
     778    ''' 
     779 
     780    if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio) 
     781     
     782    ## For the sake of simplicity, switch to the convention of original lbc Fortran routine from NEMO 
     783    ## : starts indexes at 1 
     784    jy = jj + 1 ; ix = ii + 1 
     785 
     786    math = __math__ (jj) 
     787    if math == None : math=np 
     788 
     789    # 
     790    #> East-West boundary conditions 
     791    # ------------------------------ 
     792    if nperio in [1, 4, 6] : 
     793        #... cyclic 
     794        ix = math.where (ix==jpi, 2   , ix) 
     795        ix = math.where (ix== 1 ,jpi-1, ix) 
     796 
     797    # 
     798    def modIJ (cond, jy_new, ix_new) : 
     799        jy_r = math.where (cond, jy_new, jy) 
     800        ix_r = math.where (cond, ix_new, ix) 
     801        return jy_r, ix_r 
     802    # 
     803    #> North-South boundary conditions 
     804    # -------------------------------- 
     805    if nperio in [ 3 , 4 ]  : 
     806        if cd_type in  [ 'T' , 'W' ] : 
     807            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix>=2       ), jpj-2, jpi-ix+2) 
     808            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1       ), jpj-1, 3       )    
     809            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy   , jpi-ix+2)  
     810 
     811        if cd_type in [ 'U' ] : 
     812            (jy, ix) = modIJ (np.logical_and (jy==jpj  , np.logical_and (ix>=1, ix <= jpi-1)   ), jy   , jpi-ix+1) 
     813            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1  )                               , jpj-2, 2       ) 
     814            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==jpi)                               , jpj-2, jpi-1   ) 
     815            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, np.logical_and (ix>=jpi//2, ix<=jpi-1)), jy   , jpi-ix+1) 
     816           
     817        if cd_type in [ 'V' ] : 
     818            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=2  ), jpj-2, jpi-ix+2) 
     819            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix>=2  ), jpj-3, jpi-ix+2) 
     820            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1  ), jpj-3,  3      ) 
     821             
     822        if cd_type in [ 'F' ] : 
     823            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix<=jpi-1), jpj-2, jpi-ix+1) 
     824            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix<=jpi-1), jpj-3, jpi-ix+1) 
     825            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==1    ), jpj-3, 2       ) 
     826            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==jpi  ), jpj-3, jpi-1   ) 
     827 
     828    if nperio in [ 5 , 6 ] : 
     829        if cd_type in [ 'T' , 'W' ] :                        # T-, W-point 
     830             (jy, ix) = modIJ (jy==jpj, jpj-1, jpi-ix+1) 
     831  
     832        if cd_type in [ 'U' ] :                              # U-point 
     833            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix<=jpi-1   ), jpj-1, jpi-ix  ) 
     834            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix==jpi     ), jpi-1, 1       ) 
     835             
     836        if cd_type in [ 'V' ] :    # V-point 
     837            (jy, ix) = modIJ (jy==jpj                                 , jy   , jpi-ix+1) 
     838            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy   , jpi-ix+1) 
     839             
     840        if cd_type in [ 'F' ] :                              # F-point 
     841            (jy, ix) = modIJ (np.logical_and (jy==jpj  , ix<=jpi-1   ), jpj-2, jpi-ix  ) 
     842            (jy, ix) = modIJ (np.logical_and (ix==jpj  , ix==jpi     ), jpj-2, 1       ) 
     843            (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy   , jpi-ix  ) 
     844 
     845    ## Restore convention to Python/C : indexes start at 0 
     846    jy += -1 ; ix += -1 
     847 
     848    if isinstance (jj, int) : jy = jy.item () 
     849    if isinstance (ii, int) : ix = ix.item () 
     850 
     851    return jy, ix 
     852     
     853def geo2en (pxx, pyy, pzz, glam, gphi) :  
    566854    ''' 
    567855    Change vector from geocentric to east/north 
    568856 
    569     Input 
    570         pxx, pyy, pzz : components on the geoce,tric system 
    571         glam, gphi : longitue and latitude of the points 
    572  
    573     ''' 
     857    Inputs : 
     858        pxx, pyy, pzz : components on the geocentric system 
     859        glam, gphi : longitude and latitude of the points 
     860    ''' 
     861 
    574862    gsinlon = np.sin (rad * glam) 
    575863    gcoslon = np.cos (rad * glam) 
    … …  
    582870    return pte, ptn 
    583871 
    584 def oce2geo (pte, ptn, glam, gphi) : 
    585     ##---------------------------------------------------------------------- 
    586     ##                    ***  ROUTINE oce2geo  *** 
    587     ##       
    588     ## ** Purpose : 
    589     ## 
    590     ## ** Method  :   Change vector from east/north to geocentric 
    591     ## 
    592     ## History : 
    593     ##        #         (A. Caubel)  oce2geo - Original code 
    594     ##---------------------------------------------------------------------- 
     872def en2geo (pte, ptn, glam, gphi) : 
     873    ''' 
     874    Change vector from east/north to geocentric 
     875 
     876    Inputs :  
     877        pte, ptn   : eastward/northward components 
     878        glam, gphi : longitude and latitude of the points 
     879    ''' 
     880     
    595881    gsinlon = np.sin (rad * glam) 
    596882    gcoslon = np.cos (rad * glam) 
    … …  
    603889     
    604890    return pxx, pyy, pzz 
     891 
     892def findJI (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False) : 
     893    ''' 
     894    Description: seeks J,I indices of the grid point which is the closest of a given point  
     895    Usage: go FindJI  <data latitude> <data longitude> <grid latitudes> <grid longitudes> [mask] 
     896    <longitude fields> <latitude field> are 2D fields on J/I (Y/X) dimensions 
     897    mask : if given, seek only non masked grid points (i.e with mask=1) 
     898     
     899    Example : findIJ (40, -20, nav_lat, nav_lon, mask=1.0) 
     900 
     901    Note : all longitudes and latitudes in degrees 
     902         
     903    Note : may work with 1D lon/lat (?) 
     904    ''' 
     905    # Get grid dimensions 
     906    if len (lon_grid.shape) == 2 : (jpj, jpi) = lon_grid.shape 
     907    else                         : jpj = len(lat_grid) ; jpi=len(lon_grid) 
     908 
     909    math = __math__ (lat_grid) 
     910         
     911    # Compute distance from the point to all grid points (in radian) 
     912    arg      = np.sin (rad*lat_data) * np.sin (rad*lat_grid) \ 
     913             + np.cos (rad*lat_data) * np.cos (rad*lat_grid) * np.cos(rad*(lon_data-lon_grid)) 
     914    distance = np.arccos (arg) + 4.0*rpi*(1.0-mask) # Send masked points to 'infinite'  
     915 
     916    # Truncates to alleviate some precision problem with some grids 
     917    prec = int (1E7) 
     918    distance = (distance*prec).astype(int) / prec 
     919 
     920    # Compute minimum of distance, and index of minimum 
     921    # 
     922    distance_min = distance.min    () 
     923    jimin        = int (distance.argmin ()) 
     924     
     925    # Compute 2D indices  
     926    jmin = jimin // jpi ; imin = jimin - jmin*jpi 
     927 
     928    # Compute distance achieved 
     929    mindist = distance[jmin, imin] 
     930     
     931    # Compute azimuth 
     932    dlon = lon_data-lon_grid[jmin,imin] 
     933    arg  = np.sin (rad*dlon) /  (np.cos(rad*lat_data)*np.tan(rad*lat_grid[jmin,imin]) - np.sin(rad*lat_data)*np.cos(rad*dlon)) 
     934    azimuth = dar*np.arctan (arg) 
     935     
     936    # Result 
     937    if verbose :  
     938        print ('I={:d} J={:d} - Data:{:5.1f}°N {:5.1f}°E - Grid:{:4.1f}°N {:4.1f}°E - Dist: {:6.1f}km {:5.2f}° - Azimuth: {:3.2f}rad - {:5.1f}°' 
     939            .format (imin, jmin, lat_data, lon_data, lat_grid[jmin,imin], lon_grid[jmin,imin], ra*distance[jmin,imin], dar*distance[jmin,imin], rad*azimuth, azimuth)) 
     940 
     941    return jmin, imin 
     942 
     943def clo_lon (lon, lon0) : 
     944    '''Choose closest to lon0 longitude, adding or substacting 360° if needed''' 
     945    math = __math__ (lon, np) 
     946         
     947    clo_lon = lon 
     948    clo_lon = math.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon) 
     949    clo_lon = math.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon) 
     950    clo_lon = math.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon) 
     951    clo_lon = math.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon) 
     952    if clo_lon.shape == () : clo_lon = clo_lon.item () 
     953    return clo_lon 
     954 
     955def angle_full (glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif, nperio=None) : 
     956    '''Compute sinus and cosinus of model line direction with respect to east''' 
     957    math = __math__ (glamt) 
     958     
     959    # north pole direction & modulous (at T-point) 
     960    zxnpt = 0. - 2.0 * np.cos (rad*glamt) * np.tan (rpi/4.0 - rad*gphit/2.0) 
     961    zynpt = 0. - 2.0 * np.sin (rad*glamt) * np.tan (rpi/4.0 - rad*gphit/2.0) 
     962    znnpt = zxnpt*zxnpt + zynpt*zynpt 
     963     
     964    # north pole direction & modulous (at U-point) 
     965    zxnpu = 0. - 2.0 * np.cos (rad*glamu) * np.tan (rpi/4.0 - rad*gphiu/2.0) 
     966    zynpu = 0. - 2.0 * np.sin (rad*glamu) * np.tan (rpi/4.0 - rad*gphiu/2.0) 
     967    znnpu = zxnpu*zxnpu + zynpu*zynpu 
     968     
     969    # north pole direction & modulous (at V-point) 
     970    zxnpv = 0. - 2.0 * np.cos (rad*glamv) * np.tan (rpi/4.0 - rad*gphiv/2.0) 
     971    zynpv = 0. - 2.0 * np.sin (rad*glamv) * np.tan (rpi/4.0 - rad*gphiv/2.0) 
     972    znnpv = zxnpv*zxnpv + zynpv*zynpv 
     973 
     974    # north pole direction & modulous (at F-point) 
     975    zxnpf = 0. - 2.0 * np.cos( rad*glamf ) * np.tan ( rpi/4. - rad*gphif/2. ) 
     976    zynpf = 0. - 2.0 * np.sin( rad*glamf ) * np.tan ( rpi/4. - rad*gphif/2. ) 
     977    znnpf = zxnpf*zxnpf + zynpf*zynpf 
     978 
     979    # j-direction: v-point segment direction (around T-point) 
     980    zlam = glamv   
     981    zphi = gphiv 
     982    zlan = np.roll ( glamv, axis=-2, shift=1)  # glamv (ji,jj-1) 
     983    zphh = np.roll ( gphiv, axis=-2, shift=1)  # gphiv (ji,jj-1) 
     984    zxvvt =  2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     985          -  2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     986    zyvvt =  2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     987          -  2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     988    znvvt = np.sqrt ( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt )  ) 
     989 
     990    # j-direction: f-point segment direction (around u-point) 
     991    zlam = glamf 
     992    zphi = gphif 
     993    zlan = np.roll (glamf, axis=-2, shift=1) # glamf (ji,jj-1) 
     994    zphh = np.roll (gphif, axis=-2, shift=1) # gphif (ji,jj-1) 
     995    zxffu =  2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     996          -  2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     997    zyffu =  2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     998          -  2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     999    znffu = np.sqrt ( znnpu * ( zxffu*zxffu + zyffu*zyffu )  ) 
     1000 
     1001    # i-direction: f-point segment direction (around v-point) 
     1002    zlam = glamf   
     1003    zphi = gphif 
     1004    zlan = np.roll (glamf, axis=-1, shift=1) # glamf (ji-1,jj) 
     1005    zphh = np.roll (gphif, axis=-1, shift=1) # gphif (ji-1,jj) 
     1006    zxffv =  2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     1007          -  2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     1008    zyffv =  2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     1009          -  2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     1010    znffv = np.sqrt ( znnpv * ( zxffv*zxffv + zyffv*zyffv )  ) 
     1011 
     1012    # j-direction: u-point segment direction (around f-point) 
     1013    zlam = np.roll (glamu, axis=-2, shift=-1) # glamu (ji,jj+1)  
     1014    zphi = np.roll (gphiu, axis=-2, shift=-1) # gphiu (ji,jj+1) 
     1015    zlan = glamu 
     1016    zphh = gphiu 
     1017    zxuuf =  2. * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     1018          -  2. * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     1019    zyuuf =  2. * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )   \ 
     1020          -  2. * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. ) 
     1021    znuuf = np.sqrt ( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf )  ) 
     1022 
     1023     
     1024    # cosinus and sinus using scalar and vectorial products 
     1025    gsint = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 
     1026    gcost = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 
     1027     
     1028    gsinu = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 
     1029    gcosu = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 
     1030     
     1031    gsinf = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 
     1032    gcosf = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 
     1033     
     1034    gsinv = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 
     1035    gcosv =-( zxnpv*zyffv - zynpv*zxffv ) / znffv  # (caution, rotation of 90 degres) 
     1036     
     1037    gsint = lbc (gsint, cd_type='T', nperio=nperio, psgn=-1.) 
     1038    gcost = lbc (gcost, cd_type='T', nperio=nperio, psgn=-1.) 
     1039    gsinu = lbc (gsinu, cd_type='U', nperio=nperio, psgn=-1.) 
     1040    gcosu = lbc (gcosu, cd_type='U', nperio=nperio, psgn=-1.) 
     1041    gsinv = lbc (gsinv, cd_type='V', nperio=nperio, psgn=-1.) 
     1042    gcosv = lbc (gcosv, cd_type='V', nperio=nperio, psgn=-1.) 
     1043    gsinf = lbc (gsinf, cd_type='F', nperio=nperio, psgn=-1.) 
     1044    gcosf = lbc (gcosf, cd_type='F', nperio=nperio, psgn=-1.) 
     1045 
     1046    if math == xr : 
     1047        gsint = gsint.assign_coords ( glamt.coords ) 
     1048        gcost = gcost.assign_coords ( glamt.coords ) 
     1049        gsinu = gsinu.assign_coords ( glamu.coords ) 
     1050        gcosu = gcosu.assign_coords ( glamu.coords ) 
     1051        gsinv = gsinv.assign_coords ( glamv.coords ) 
     1052        gcosv = gcosv.assign_coords ( glamv.coords ) 
     1053        gsinf = gsinf.assign_coords ( glamf.coords ) 
     1054        gcosf = gcosf.assign_coords ( glamf.coords ) 
     1055 
     1056    return gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf 
     1057 
     1058def angle (glam, gphi, nperio, cd_type='T') : 
     1059    '''Compute sinus and cosinus of model line direction with respect to east''' 
     1060    math = __math__ (glam)  
     1061    # north pole direction & modulous 
     1062    zxnp = 0. - 2.0 * np.cos (rad*glam) * np.tan (rpi/4.0 - rad*gphi/2.0) 
     1063    zynp = 0. - 2.0 * np.sin (rad*glam) * np.tan (rpi/4.0 - rad*gphi/2.0) 
     1064    znnp = zxnp*zxnp + zynp*zynp 
     1065 
     1066    # j-direction: segment direction (around point) 
     1067    zlan_n = np.roll (glam, axis=-2, shift=-1) # glam [jj+1, ji] 
     1068    zphh_n = np.roll (gphi, axis=-2, shift=-1) # gphi [jj+1, ji] 
     1069    zlan_s = np.roll (glam, axis=-2, shift= 1) # glam [jj-1, ji] 
     1070    zphh_s = np.roll (gphi, axis=-2, shift= 1) # gphi [jj-1, ji] 
     1071     
     1072    zxff = 2.0 * np.cos (rad*zlan_n) * np.tan (rpi/4.0 - rad*zphh_n/2.0) \ 
     1073        -  2.0 * np.cos (rad*zlan_s) * np.tan (rpi/4.0 - rad*zphh_s/2.0) 
     1074    zyff = 2.0 * np.sin (rad*zlan_n) * np.tan (rpi/4.0 - rad*zphh_n/2.0) \ 
     1075        -  2.0 * np.sin (rad*zlan_s) * np.tan (rpi/4.0 - rad*zphh_s/2.0) 
     1076    znff = np.sqrt (znnp * (zxff*zxff + zyff*zyff) ) 
     1077  
     1078    gsin = (zxnp*zyff - zynp*zxff) / znff 
     1079    gcos = (zxnp*zxff + zynp*zyff) / znff 
     1080 
     1081    gsin = lbc (gsin, cd_type=cd_type, nperio=nperio, psgn=-1.) 
     1082    gcos = lbc (gcos, cd_type=cd_type, nperio=nperio, psgn=-1.) 
     1083 
     1084    if math == xr : 
     1085        gsin = gsin.assign_coords ( glam.coords ) 
     1086        gcos = gcos.assign_coords ( glam.coords ) 
     1087         
     1088    return gsin, gcos 
     1089 
     1090def rot_en2ij ( u_e, v_n, gsin, gcos, nperio, cd_type ) : 
     1091    ''' 
     1092    ** Purpose :   Rotate the Repere: Change vector componantes between 
     1093    geographic grid --> stretched coordinates grid. 
     1094    All components are on the same grid (T, U, V or F)  
     1095    ''' 
     1096 
     1097    u_i = + u_e * gcos + v_n * gsin 
     1098    v_j = - u_e * gsin + v_n * gcos 
     1099     
     1100    u_i = lbc (u_i, nperio=nperio, cd_type=cd_type, psgn=-1.0) 
     1101    v_j = lbc (v_j, nperio=nperio, cd_type=cd_type, psgn=-1.0) 
     1102     
     1103    return u_i, v_j 
     1104 
     1105def rot_ij2en ( u_i, v_j, gsin, gcos, nperio, cd_type='T' ) : 
     1106    ''' 
     1107    ** Purpose :   Rotate the Repere: Change vector componantes from 
     1108    stretched coordinates grid --> geographic grid 
     1109    All components are on the same grid (T, U, V or F)  
     1110    ''' 
     1111    u_e = + u_i * gcos - v_j * gsin 
     1112    v_n = + u_i * gsin + v_j * gcos 
     1113     
     1114    u_e = lbc (u_e, nperio=nperio, cd_type=cd_type, psgn= 1.0) 
     1115    v_n = lbc (v_n, nperio=nperio, cd_type=cd_type, psgn= 1.0) 
     1116     
     1117    return u_e, v_n 
     1118 
     1119def rot_uv2en ( uo, vo, gsint, gcost, nperio, zdim='deptht' ) : 
     1120    ''' 
     1121    ** Purpose :   Rotate the Repere: Change vector componantes from 
     1122    stretched coordinates grid --> geographic grid 
     1123    uo is on the U grid point, vo is on the V grid point 
     1124    east-north components on the T grid point    
     1125    ''' 
     1126    math = __math__ (uo) 
     1127 
     1128    ut = U2T (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 
     1129    vt = V2T (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 
     1130     
     1131    u_e = + ut * gcost - vt * gsint 
     1132    v_n = + ut * gsint + vt * gcost 
     1133 
     1134    u_e = lbc (u_e, nperio=nperio, cd_type='T', psgn=1.0) 
     1135    v_n = lbc (v_n, nperio=nperio, cd_type='T', psgn=1.0) 
     1136     
     1137    return u_e, v_n 
     1138 
     1139def rot_uv2enF ( uo, vo, gsinf, gcosf, nperio, zdim='deptht' ) : 
     1140    ''' 
     1141    ** Purpose :   Rotate the Repere: Change vector componantes from 
     1142    stretched coordinates grid --> geographic grid 
     1143    uo is on the U grid point, vo is on the V grid point 
     1144    east-north components on the T grid point    
     1145    ''' 
     1146    math = __math__ (uo) 
     1147 
     1148    uf = U2F (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 
     1149    vf = V2F (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 
     1150     
     1151    u_e = + uf * gcosf - vf * gsinf 
     1152    v_n = + uf * gsinf + vf * gcosf 
     1153 
     1154    u_e = lbc (u_e, nperio=nperio, cd_type='F', psgn= 1.0) 
     1155    v_n = lbc (v_n, nperio=nperio, cd_type='F', psgn= 1.0) 
     1156     
     1157    return u_e, v_n 
     1158              
     1159def U2T (utab, nperio=None, psgn=-1.0, zdim='deptht') : 
     1160    '''Interpolate an array from U grid to T grid i-mean)''' 
     1161    math = __math__ (utab) 
     1162    utab_0 = math.where ( np.isnan(utab), 0., utab) 
     1163    utab_0 = lbc (utab_0 , nperio=nperio, cd_type='U', psgn=psgn) 
     1164    ix, ax = __findAxis__ (utab_0, 'x') 
     1165    iz, az = __findAxis__ (utab_0, 'z') 
     1166    ttab =  lbc ( 0.5 * (utab_0 + np.roll (utab_0, axis=ix, shift=1)), cd_type='T', nperio=nperio, psgn=psgn) 
     1167     
     1168    if math == xr : 
     1169        if ax != None : 
     1170            ttab = ttab.assign_coords({ax:np.arange (ttab.shape[ix])+1.}) 
     1171        if zdim != None and iz != None  and az != 'olevel' :  
     1172            ttab = ttab.rename( {az:zdim})  
     1173    return ttab 
     1174 
     1175def V2T (vtab, nperio=None, psgn=-1.0, zdim='deptht') : 
     1176    '''Interpolate an array from V grid to T grid (j-mean)''' 
     1177    math = __math__ (vtab) 
     1178    vtab_0 = math.where ( np.isnan(vtab), 0., vtab) 
     1179    vtab_0 = lbc (vtab_0 , nperio=nperio, cd_type='V', psgn=psgn) 
     1180    iy, ay = __findAxis__ (vtab_0, 'y') 
     1181    iz, az = __findAxis__ (vtab_0, 'z') 
     1182    ttab = lbc ( 0.5 * (vtab_0 + np.roll (vtab_0, axis=iy, shift=1)), cd_type='T', nperio=nperio, psgn=psgn) 
     1183                      
     1184    if math == xr : 
     1185        if ay !=None :  
     1186            ttab = ttab.assign_coords({ay:np.arange(ttab.shape[iy])+1.}) 
     1187        if zdim != None and iz != None  and az != 'olevel' : 
     1188            ttab = ttab.rename( {az:zdim})  
     1189    return ttab 
     1190 
     1191def F2T (ftab, nperio=None, psgn=1.0, zdim='depthf') : 
     1192    '''Interpolate an array from F grid to T grid (i- and j- means)''' 
     1193    math = __math__ (ftab) 
     1194    ftab_0 = math.where ( np.isnan(ftab), 0., ftab) 
     1195    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 
     1196    ttab = V2T(F2V(ftab_0, nperio=nperio, psgn=psgn, zdim=zdim), nperio=nperio, psgn=psgn, zdim=zdim) 
     1197    return ttab 
     1198 
     1199def T2U (ttab, nperio=None, psgn=1.0, zdim='depthu') : 
     1200    '''Interpolate an array from T grid to U grid (i-mean)''' 
     1201    math = __math__ (ttab) 
     1202    ttab_0 = math.where ( np.isnan(ttab), 0., ttab) 
     1203    ttab_0 = lbc (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn) 
     1204    ix, ax = __findAxis__ (ttab_0, 'x') 
     1205    iz, az = __findAxis__ (ttab_0, 'z') 
     1206    utab = lbc ( 0.5 * (ttab_0 + np.roll (ttab_0, axis=ix, shift=-1)), cd_type='U', nperio=nperio, psgn=psgn) 
     1207 
     1208    if math == xr :     
     1209        if ax != None :  
     1210            utab = ttab.assign_coords({ax:np.arange(utab.shape[ix])+1.}) 
     1211        if zdim != None  and iz != None  and az != 'olevel' : 
     1212            utab = utab.rename( {az:zdim})  
     1213    return utab 
     1214 
     1215def T2V (ttab, nperio=None, psgn=1.0, zdim='depthv') : 
     1216    '''Interpolate an array from T grid to V grid (j-mean)''' 
     1217    math = __math__ (ttab) 
     1218    ttab_0 = math.where ( np.isnan(ttab), 0., ttab) 
     1219    ttab_0 = lbc (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn) 
     1220    iy, ay = __findAxis__ (ttab_0, 'y') 
     1221    iz, az = __findAxis__ (ttab_0, 'z') 
     1222    vtab = lbc ( 0.5 * (ttab_0 + np.roll (ttab_0, axis=iy, shift=-1)), cd_type='V', nperio=nperio, psgn=psgn) 
     1223 
     1224    if math == xr : 
     1225        if ay != None :  
     1226            vtab = vtab.assign_coords({ay:np.arange(vtab.shape[iy])+1.}) 
     1227        if zdim != None  and iz != None and az != 'olevel' : 
     1228            vtab = vtab.rename( {az:zdim})  
     1229    return vtab 
     1230 
     1231def V2F (vtab, nperio=None, psgn=-1.0, zdim='depthf') : 
     1232    '''Interpolate an array from V grid to F grid (i-mean)''' 
     1233    math = __math__ (vtab) 
     1234    vtab_0 = math.where ( np.isnan(vtab), 0., vtab) 
     1235    vtab_0 = lbc (vtab_0 , nperio=nperio, cd_type='V', psgn=psgn) 
     1236    ix, ax = __findAxis__ (vtab_0, 'x') 
     1237    iz, az = __findAxis__ (vtab_0, 'z') 
     1238    ftab = lbc ( 0.5 * (vtab_0 + np.roll (vtab_0, axis=ix, shift=-1)), cd_type='F', nperio=nperio, psgn=psgn) 
     1239 
     1240    if math == xr : 
     1241        if ax != None :  
     1242            ftab = ftab.assign_coords({ax:np.arange(ftab.shape[ix])+1.}) 
     1243        if zdim != None and iz != None and az != 'olevel' : 
     1244            ftab = ftab.rename( {az:zdim})  
     1245    return ftab 
     1246 
     1247def U2F (utab, nperio=None, psgn=-1.0, zdim='depthf') : 
     1248    '''Interpolate an array from U grid to F grid i-mean)''' 
     1249    math = __math__ (utab) 
     1250    utab_0 = math.where ( np.isnan(utab), 0., utab) 
     1251    utab_0 = lbc (utab_0 , nperio=nperio, cd_type='U', psgn=psgn) 
     1252    iy, ay = __findAxis__ (utab_0, 'y') 
     1253    iz, az = __findAxis__ (utab_0, 'z') 
     1254    ftab = lbc ( 0.5 * (utab_0 + np.roll (utab_0, axis=iy, shift=-1)), cd_type='F', nperio=nperio, psgn=psgn) 
     1255 
     1256    if math == xr : 
     1257        if ay != None :  
     1258            ftab = ftab.assign_coords({'y':np.arange(ftab.shape[iy])+1.}) 
     1259        if zdim != None and iz != None and az != 'olevel' : 
     1260            ftab = ftab.rename( {az:zdim})  
     1261    return ftab 
     1262 
     1263def F2T (ftab, nperio=None, psgn=1.0, zdim='deptht') : 
     1264    '''Interpolate an array on F grid to T grid (i- and j- means)''' 
     1265    math = __math__ (ftab) 
     1266    ftab_0 = math.where ( np.isnan(ttab), 0., ttab) 
     1267    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 
     1268    ttab = U2T(F2U(ftab_0, nperio=nperio, psgn=psgn, zdim=zdim), nperio=nperio, psgn=psgn, zdim=zdim) 
     1269    return ttab 
     1270 
     1271def T2F (ttab, nperio=None, psgn=1.0, zdim='deptht') : 
     1272    '''Interpolate an array on T grid to F grid (i- and j- means)''' 
     1273    math = __math__ (ftab) 
     1274    ttab_0 = math.where ( np.isnan(ttab), 0., ttab) 
     1275    ttab_0 = lbc (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn) 
     1276    ftab = T2U(U2F(ttab, nperio=nperio, psgn=psgn, zdim=zdim), nperio=nperio, psgn=psgn, zdim=zdim) 
     1277    return ftab 
     1278 
     1279def F2U (ftab, nperio=None, psgn=1.0, zdim='depthu') : 
     1280    '''Interpolate an array on F grid to FUgrid (i-mean)''' 
     1281    math = __math__ (ftab) 
     1282    ftab_0 = math.where ( np.isnan(ftab), 0., ftab) 
     1283    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 
     1284    iy, ay = __findAxis__ (ftab_0, 'y') 
     1285    iz, az = __findAxis__ (ftab_0, 'z') 
     1286    utab = lbc ( 0.5 * (ftab_0 + np.roll (ftab_0, axis=iy, shift=-1)), cd_type='U', nperio=nperio, psgn=psgn) 
     1287     
     1288    if math == xr : 
     1289        utab = utab.assign_coords({ay:np.arange(ftab.shape[iy])+1.}) 
     1290        if zdim != None and iz != None and az != 'olevel' : 
     1291            utab = utab.rename( {az:zdim})  
     1292    return utab 
     1293 
     1294def F2V (ftab, nperio=None, psgn=1.0, zdim='depthv') : 
     1295    '''Interpolate an array from F grid to V grid (i-mean)''' 
     1296    math = __math__ (ftab) 
     1297    ftab_0 = math.where ( np.isnan(ftab), 0., ftab) 
     1298    ftab_0 = lbc (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 
     1299    ix, ax = __findAxis__ (ftab_0, 'x') 
     1300    iz, az = __findAxis__ (ftab_0, 'z') 
     1301    vtab = lbc ( 0.5 * (ftab_0 + np.roll (ftab_0, axis=ix, shift=-1)), cd_type='V', nperio=nperio) 
     1302 
     1303    if math == xr : 
     1304        vtab = vtab.assign_coords({ax:np.arange(ftab.shape[ix])+1.}) 
     1305        if zdim != None and iz != None and az != 'olevel' : 
     1306            vtab = vtab.rename( {az:zdim})  
     1307    return vtab 
     1308 
     1309def W2T (wtab, zcoord=None, zdim='deptht', sval=np.nan) : 
     1310    ''' 
     1311    Interpolate an array on W grid to T grid (k-mean) 
     1312    sval is the bottom value 
     1313    ''' 
     1314    math = __math__ (wtab) 
     1315    wtab_0 = math.where ( np.isnan(wtab), 0., wtab) 
     1316 
     1317    iz, az = __findAxis__ (wtab_0, 'z') 
     1318        
     1319    ttab = 0.5 * ( wtab_0 + np.roll (wtab_0, axis=iz, shift=-1) ) 
     1320     
     1321    if math == xr : 
     1322        ttab[{az:iz}] = sval 
     1323        if zdim != None and iz != None and az != 'olevel' : 
     1324            ttab = ttab.rename ( {az:zdim} ) 
     1325        try    : ttab = ttab.assign_coords ( {zdim:zcoord} ) 
     1326        except : pass 
     1327    else : 
     1328        ttab[..., -1, :, :] = sval 
     1329 
     1330    return ttab 
     1331 
     1332def T2W (ttab, zcoord=None, zdim='depthw', sval=np.nan, extrap_surf=False) : 
     1333    '''Interpolate an array from T grid to W grid (k-mean) 
     1334    sval is the surface value 
     1335    if extrap_surf==True, surface value is taken from 1st level value. 
     1336    ''' 
     1337    math = __math__ (ttab) 
     1338    ttab_0 = math.where ( np.isnan(ttab), 0., ttab) 
     1339    iz, az = __findAxis__ (ttab_0, 'z') 
     1340    wtab = 0.5 * ( ttab_0 + np.roll (ttab_0, axis=iz, shift=1) ) 
     1341 
     1342    if math == xr : 
     1343        if extrap_surf : wtab[{az:0}] = ttabb[{az:0}] 
     1344        else           : wtab[{az:0}] = sval 
     1345    else :  
     1346        if extrap_surf : wtab[..., 0, :, :] = ttab[..., 0, :, :] 
     1347        else           : wtab[..., 0, :, :] = sval 
     1348 
     1349    if math == xr : 
     1350        if zdim != None and iz != None and az != 'olevel' : 
     1351                wtab = wtab.rename ( {az:zdim}) 
     1352        if zcoord != None : wtab = wtab.assign_coords ( {zdim:zcoord}) 
     1353        else              : ztab = wtab.assign_coords ( {zdim:np.arange(ttab.shape[iz])+1.} ) 
     1354    return wtab 
     1355 
     1356def fill (ptab, nperio, cd_type='T', npass=1) : 
     1357    ''' 
     1358    Fill 0. or np.nan values with mean of neighbours 
     1359    
     1360    Inputs : 
     1361       ptab : input field to fill 
     1362       nperio, cd_type : periodicity characteristics 
     1363    '''        
     1364 
     1365    math = __math__ (ptab) 
     1366    ztab   = math.where (np.isnan(ptab), 0., ptab) 
     1367 
     1368    DoPerio = False 
     1369    if nperio == 4.2 or nperio == 6.2 : DoPerio = True 
     1370         
     1371    if DoPerio :  ztab = lbc_add (ztab) 
     1372    
     1373    for nn in np.arange (npass) :  
     1374        zmask  = math.where (ztab==0., 0., 1.  ) 
     1375        # Compte du nombre de voisins 
     1376        zcount = zmask \ 
     1377          + np.roll(zmask, shift=1, axis=-1) + np.roll(zmask, shift=-1, axis=-1) \ 
     1378          + np.roll(zmask, shift=1, axis=-2) + np.roll(zmask, shift=-1, axis=-2) \ 
     1379          + 0.5 * ( \ 
     1380                + np.roll(np.roll(zmask, shift= 1, axis=-2), shift= 1, axis=-1) \ 
     1381                + np.roll(np.roll(zmask, shift=-1, axis=-2), shift= 1, axis=-1) \ 
     1382                + np.roll(np.roll(zmask, shift= 1, axis=-2), shift=-1, axis=-1) \ 
     1383                + np.roll(np.roll(zmask, shift=-1, axis=-2), shift=-1, axis=-1) ) 
     1384 
     1385        znew = zmask \ 
     1386          + np.roll(ztab, shift=1, axis=-1) + np.roll(ztab, shift=-1, axis=-1) \ 
     1387          + np.roll(ztab, shift=1, axis=-2) + np.roll(ztab, shift=-1, axis=-2) \ 
     1388          + 0.5 * ( \ 
     1389            + np.roll(np.roll(ztab, shift= 1, axis=-2), shift= 1, axis=-1) \ 
     1390            + np.roll(np.roll(ztab, shift=-1, axis=-2), shift= 1, axis=-1) \ 
     1391            + np.roll(np.roll(ztab, shift= 1, axis=-2), shift=-1, axis=-1) \ 
     1392            + np.roll(np.roll(ztab, shift=-1, axis=-2), shift=-1, axis=-1) ) 
     1393 
     1394        zcount = lbc (zcount, nperio=nperio, cd_type=cd_type) 
     1395        znew   = lbc (znew  , nperio=nperio, cd_type=cd_type) 
     1396     
     1397        ztab = math.where (np.logical_and (zmask==0., zcount>0), znew/zcount, ztab) 
     1398         
     1399    ztab = math.where (ztab==0., np.nan, ztab) 
     1400 
     1401    if DoPerio : ztab = lbc_del (ztab) 
     1402 
     1403    return ztab 
     1404 
     1405def correct_uv (u, v, lat) : 
     1406    ''' 
     1407    Correct a Cartopy bug in Orthographic projection 
     1408 
     1409    See https://github.com/SciTools/cartopy/issues/1179 
     1410 
     1411    The correction is needed with cartopy <= 0.20 
     1412    It seems that version 0.21 will correct the bug (https://github.com/SciTools/cartopy/pull/1926) 
     1413 
     1414    Inputs : 
     1415       u, v : eastward/nothward components 
     1416       lat  : latitude of the point (degrees north) 
     1417 
     1418    Outputs :  
     1419       modified eastward/nothward components have correct polar projection in cartopy 
     1420    ''' 
     1421    math = __math__ (u) 
     1422    uv = np.sqrt (u*u + v*v)           # Original modulus 
     1423    zu = u 
     1424    zv = v * np.cos (rad*lat) 
     1425    zz = np.sqrt ( zu*zu + zv*zv )     # Corrected modulus 
     1426    uc = zu*uv/zz ; vc = zv*uv/zz      # Final corrected values 
     1427    return uc, vc 
    6051428 
    6061429## =========================================================================== 
    … …  
    6091432## 
    6101433## =========================================================================== 
     1434 
     1435def is_orca_north_fold ( Xtest, cname_long='T' ) : 
     1436    ''' 
     1437    Ported (pirated !!?) from Sosie 
     1438 
     1439    Tell if there is a 2/point band overlaping folding at the north pole typical of the ORCA grid 
     1440 
     1441    0 => not an orca grid (or unknown one) 
     1442    4 => North fold T-point pivot (ex: ORCA2) 
     1443    6 => North fold F-point pivot (ex: ORCA1) 
     1444 
     1445    We need all this 'cname_long' stuff because with our method, there is a 
     1446    confusion between "Grid_U with T-fold" and "Grid_V with F-fold" 
     1447    => so knowing the name of the longitude array (as in namelist, and hence as 
     1448    in netcdf file) might help taking the righ decision !!! UGLY!!! 
     1449    => not implemented yet 
     1450    ''' 
     1451     
     1452    ifld_nord =  0 ; cgrd_type = 'X' 
     1453    ny, nx = Xtest.shape[-2:] 
     1454 
     1455    if ny > 3 : # (case if called with a 1D array, ignoring...) 
     1456        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. : 
     1457          ifld_nord = 4 ; cgrd_type = 'T' # T-pivot, grid_T       
     1458 
     1459        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2  :-1] ).sum() == 0. : 
     1460            if cnlon == 'U' : ifld_nord = 4 ;  cgrd_type = 'U' # T-pivot, grid_T 
     1461                ## LOLO: PROBLEM == 6, V !!! 
     1462 
     1463        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. : 
     1464            ifld_nord = 4 ; cgrd_type = 'V' # T-pivot, grid_V 
     1465 
     1466        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-2, nx-1-1:nx-nx//2:-1] ).sum() == 0. : 
     1467            ifld_nord = 6 ; cgrd_type = 'T'# F-pivot, grid_T 
     1468 
     1469        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-1, nx-1:nx-nx//2-1:-1] ).sum() == 0. : 
     1470            ifld_nord = 6 ;  cgrd_type = 'U' # F-pivot, grid_U 
     1471 
     1472        if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2  :-1] ).sum() == 0. : 
     1473            if cnlon == 'V' : ifld_nord = 6 ; cgrd_type = 'V' # F-pivot, grid_V 
     1474                ## LOLO: PROBLEM == 4, U !!! 
     1475 
     1476    return ifld_nord, cgrd_type 

  • TOOLS/MOSAIX/update_xml.py

    r6093 r6190  
    66### =========================================================================== 
    77## 
    8 ##  Warning, to install, configure, run, use any of Olivier Marti's 
    9 ##  software or to read the associated documentation you'll need at least 
    10 ##  one (1) brain in a reasonably working order. Lack of this implement 
    11 ##  will void any warranties (either express or implied). 
    12 ##  O. Marti assumes no responsability for errors, omissions, 
    13 ##  data loss, or any other consequences caused directly or indirectly by 
    14 ##  the usage of his software by incorrectly or partially configured 
    15 ##  personal. 
     8##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" 
     9##  file for an english version of the licence and 
     10##  "Licence_CeCILL_V2-fr.txt" for a french version. 
     11## 
     12##  Permission is hereby granted, free of charge, to any person or 
     13##  organization obtaining a copy of the software and accompanying 
     14##  documentation covered by this license (the "Software") to use, 
     15##  reproduce, display, distribute, execute, and transmit the 
     16##  Software, and to prepare derivative works of the Software, and to 
     17##  permit third-parties to whom the Software is furnished to do so, 
     18##  all subject to the following: 
     19## 
     20##  Warning, to install, configure, run, use any of MOSAIX software or 
     21##  to read the associated documentation you'll need at least one (1) 
     22##  brain in a reasonably working order. Lack of this implement will 
     23##  void any warranties (either express or implied).  Authors assumes 
     24##  no responsability for errors, omissions, data loss, or any other 
     25##  consequences caused directly or indirectly by the usage of his 
     26##  software by incorrectly or partially configured 
     27## 
    1628## 
    1729## SVN information 
Note: See TracChangeset for help on using the changeset viewer.