Changeset 6190 for TOOLS/MOSAIX
- Timestamp:
- 07/06/22 11:06:07 (23 months ago)
- Location:
- TOOLS/MOSAIX
- Files:
-
- 13 edited
Legend:
- Unmodified
- Added
- Removed
-
TOOLS/MOSAIX/Build_coordinates_mask.py
r6124 r6190 6 6 # creates file coordinates_mask.nc used by mosaix from NEMO netcdf files 7 7 # 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 ## 9 30 import numpy as np 10 31 import xarray as xr -
TOOLS/MOSAIX/CalvingWeights.py
r6105 r6190 5 5 ### =========================================================================== 6 6 ## 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 ## 15 27 ## 16 28 import numpy as np, xarray as xr … … 37 49 # Adding arguments 38 50 parser.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'] ) 40 52 parser.add_argument ('--atm' , help='atm model name (ICO* or LMD*)', type=str, default='ICO40' ) 41 53 parser.add_argument ('--type' , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full' ) -
TOOLS/MOSAIX/ComputeNemoCoast.py
r3901 r6190 5 5 ### =========================================================================== 6 6 ## 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 ## 15 27 ## 16 28 ## SVN information -
TOOLS/MOSAIX/CreateOasisGrids.bash
r6107 r6190 6 6 ### =========================================================================== 7 7 ## 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 16 27 ## 17 28 ### -
TOOLS/MOSAIX/CreateWeightsMask.bash
r6107 r6190 21 21 ### =========================================================================== 22 22 ## 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 ## 31 43 ## 32 44 ## SVN information … … 70 82 echo ${Titre}"Defines model"${Norm} 71 83 # ================================= 72 #CplModel=ORCA2.3xLMD969584 CplModel=ORCA2.3xLMD9695 73 85 #CplModel=ORCA2.3xICO30 74 86 #CplModel=ORCA2.3xICO40 … … 78 90 #CplModel=eORCA1.2xLMD256256 79 91 #CplModel=eORCA1.2xICO40 80 CplModel=eORCA1.4.2xICO4092 #CplModel=eORCA1.4.2xICO40 81 93 #CplModel=eORCA1.2xICO450 82 94 #CplModel=eORCA025.1xLMD256256 … … 898 910 echo ${Titre}"Save results"${Norm} 899 911 ## =========================================================================== 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 ) ; do901 cp ${File} ${SUBMIT_DIR}902 done912 #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 903 915 904 916 -
TOOLS/MOSAIX/README.md
r6044 r6190 16 16 or `svn co svn+ssh://omamce@forge.ipsl.jussieu.fr/ipsl/forge/projets/igcmg/svn/TOOLS/MOSAIX MOSAIX` 17 17 18 ## #[XIOS](http://forge.ipsl.jussieu.fr/ioserver/)18 ## [XIOS](http://forge.ipsl.jussieu.fr/ioserver/) 19 19 MOSAIC is known to work with XIOS revision 2286 20 20 … … 23 23 Extraction : svn co http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/trunk XIOS 24 24 25 ## Licensing 26 27 MOSAIX 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 31 Permission is hereby granted, free of charge, to any person or 32 organization obtaining a copy of the software and accompanying 33 documentation covered by this license (the "Software") to use, 34 reproduce, display, distribute, execute, and transmit the Software, 35 and to prepare derivative works of the Software, and to permit 36 third-parties to whom the Software is furnished to do so, all subject 37 to the following warning : 38 39 ### Warning 40 41 Warning, to install, configure, run, use any of MOSAIX software or to 42 read the associated documentation you'll need at least one (1) brain 43 in a reasonably working order. Lack of this implement will void any 44 warranties (either express or implied). Authors assumes no 45 responsability for errors, omissions, data loss, or any other 46 consequences caused directly or indirectly by the usage of his 47 software by incorrectly or partially configured 48 25 49 ## MOSAIX 26 50 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) 51 MOSAIX is a small set of fortran routines, bash scripts, python 52 scripts and xml files allowing to generate coupling weights for the 53 IPSL Earth System Model. It is under development, and presently not 54 used by the production runs. It aims to replace the old MOSAIC 55 software. It uses XIOS as a library to compute weights. In the present 56 state, it can handle ORCA, LMDZ lon/lat and DYNAMICO grids. 57 58 MOSAIX can generate second order conservative interpolation 59 weights. They are used in place of bilinear interpolation by 60 [OASIS-MCT](https://portal.enes.org/oasis) 30 61 31 62 As MOSAIX uses XIOS, it works in parallel, using MPI distributed memory. … … 33 64 ## Compatibility with MOSAIC 34 65 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 84 MOSAIX been developped and extensively used on Irene Skylake computer, 85 using the Intel compilation tools. MOSAIX is supposed to work on any 86 computer where XIOS works. 87 88 It has been tested on Mac OS X, using `gcc` as a compiler. It is 89 fragile, as regurlarly XIOS updates break the compatibility. 40 90 41 91 ## TODO 42 92 43 93 - 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 ... 45 97 46 98 ## Known problems … … 65 117 66 118 - 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. 68 121 - 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. 71 127 - Bash version>=4. 72 128 - [âNCO](http://nco.sourceforge.net). Version 4.6.0 is known to work. … … 74 130 ## Input files 75 131 76 MOSAIX needs files describing the grid of the models. See `CreateWeightsMask.bash` and the bottom of this page. 132 MOSAIX needs files describing the grid of the models. See 133 `CreateWeightsMask.bash` and the bottom of this page. 77 134 78 135 ## Output files … … 80 137 All netCDF files have a detailed header. 81 138 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`. 86 148 - A `README` file. Includes a checksum for each generated files. 87 149 88 150 ## Scripts and codes 89 151 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 101 177 ```bash 102 178 make_xios \[options\] … … 107 183 ``` 108 184 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. 124 216 125 217 ## Grids 126 218 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. 134 233 135 234 ## Generating input grid files -
TOOLS/MOSAIX/README.txt
r3902 r6190 20 20 21 21 22 XIOS : http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/trunk revision 1257 22 Obsolete file : please refer to README.md -
TOOLS/MOSAIX/RunoffWeights.py
r6105 r6190 8 8 ### =========================================================================== 9 9 ## 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 ## 18 30 ## 19 31 # SVN information … … 56 68 # Adding arguments 57 69 parser.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'] ) 59 71 parser.add_argument ('--atm' , help='atm model name', type=str, default='LMD9695' ) 60 72 parser.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 5 5 ### =========================================================================== 6 6 ## 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 15 26 ## 16 27 # -
TOOLS/MOSAIX/cotes_etal.py
r4186 r6190 7 7 ### =========================================================================== 8 8 ## 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 17 28 ## 18 29 # SVN information -
TOOLS/MOSAIX/make_mosaix
r6065 r6190 8 8 # $HeadURL$ 9 9 # 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 ## 10 30 11 31 # Default values for command line parameters -
TOOLS/MOSAIX/nemo.py
r6105 r6190 2 2 ## =========================================================================== 3 3 ## 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. 12 7 ## 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 ## =========================================================================== 13 25 ''' 14 26 Utilities to plot NEMO ORCA fields … … 17 29 olivier.marti@lsce.ipsl.fr 18 30 ''' 19 20 import sys, numpy as np21 try : import xarray as xr22 except : pass23 24 rpi = np.pi ; rad = rpi / 180.025 26 nperio_valid_range = [0, 1, 4, 4.2, 5, 6, 6.2]27 31 28 32 ## SVN information … … 33 37 __HeadURL = "$HeadURL$" 34 38 39 import sys, numpy as np 40 try : import xarray as xr 41 except : pass 42 43 rpi = np.pi ; rad = np.deg2rad (1.0) ; dar = np.rad2deg (1.0) 44 45 nperio_valid_range = [0, 1, 4, 4.2, 5, 6, 6.2] 46 47 rday = 24.*60.*60. # Day length [s] 48 rsiyea = 365.25 * rday * 2. * rpi / 6.283076 # Sideral year length [s] 49 rsiday = rday / (1. + rday / rsiyea) 50 raamo = 12. # Number of months in one year 51 rjjhh = 24. # Number of hours in one day 52 rhhmm = 60. # Number of minutes in one hour 53 rmmss = 60. # Number of seconds in one minute 54 omega = 2. * rpi / rsiday # Earth rotation parameter [s-1] 55 ra = 6371229. # Earth radius [m] 56 grav = 9.80665 # Gravity [m/s2] 57 repsi = np.finfo (1.0).eps 58 59 xList = [ 'x', 'X', 'lon' , 'longitude' ] 60 yList = [ 'y', 'Y', 'lat' , 'latitude' ] 61 zList = [ 'z', 'Z', 'depth' , ] 62 tList = [ 't', 'T', 'time' , ] 63 64 ## SVN information 65 __Author__ = "$Author$" 66 __Date__ = "$Date$" 67 __Revision__ = "$Revision$" 68 __Id__ = "$Id$" 69 __HeadURL = "$HeadURL$" 70 71 def __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 35 85 def __guessNperio__ (jpj, jpi, nperio=None) : 36 86 ''' … … 44 94 if nperio == None : 45 95 ## 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 53 115 ## 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 62 129 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') 64 131 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)) 66 136 67 print ('nperio set as {:} (deduced from jpi={:d})'.format (nperio, jpi))68 69 137 return nperio 70 138 … … 76 144 77 145 Inputs 78 ptab : xarray array146 ptab : xarray array 79 147 80 148 Credits : who is the original author ? 81 149 ''' 82 150 gP = None 83 if isinstance (ptab, xr.core.dataarray.DataArray) : 151 math = __math__ (ptab) 152 if math == xr : 84 153 if 'x_c' in ptab.dims and 'y_c' in ptab.dims : gP = 'T' 85 154 if 'x_f' in ptab.dims and 'y_c' in ptab.dims : gP = 'U' … … 93 162 94 163 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') 96 165 else : 97 166 print ('Grid set as', gP, 'deduced from dims ', ptab.dims) 98 167 return gP 99 168 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') 101 170 #return gP 102 171 103 def fixed_lon (lon) : 172 def __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 225 def fixed_lon (lon, center_lon=0.0) : 104 226 ''' 105 227 Returns corrected longitudes for nicer plots 106 228 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. 109 231 110 232 Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7 111 233 ''' 234 math = __math__ (lon) 235 112 236 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 115 243 116 244 # 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 119 249 120 250 return fixed_lon … … 126 256 lat : latitudes of the grid. At least 2D. 127 257 ''' 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]))) 129 266 return jeq 130 267 … … 137 274 ''' 138 275 if np.max (lat) != None : 139 je 140 lon1D = lon [..., je, :]276 je = jeq (lat) 277 lon1D = lon [..., je, :] 141 278 else : 142 jpj, jpi = lon.shape [-2:]143 lon1D = lon[..., jpj//3, :]279 jpj, jpi = lon.shape [-2:] 280 lon1D = lon [..., jpj//3, :] 144 281 145 282 start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1) 146 lon1D [..., start+1:] += 360283 lon1D [..., start+1:] += 360 147 284 148 285 return lon1D … … 155 292 diff [optional] : tolerance 156 293 ''' 294 math = __math__ (lat) 157 295 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)))) 159 297 diff = dy/100. 160 298 … … 172 310 lat : latitudes of the grid (2D) 173 311 ''' 312 math = __math__ (lat) 174 313 jpj, jpi = lat.shape[-2:] 175 try :176 if type (lat) == xr.core.dataarray.DataArray : math = xr177 else : math = np178 except : math = np179 314 180 315 dy = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2)))) … … 185 320 lat_ave = np.mean (lat, axis=-1) 186 321 187 if ( np.abs (lat_eq) < dy/100.) : # T, U or W grid322 if (np.abs (lat_eq) < dy/100.) : # T, U or W grid 188 323 dys = (90.-lat_reg) / (jpj-jreg-1)*0.5 189 324 yrange = 90.-dys-lat_reg 190 else 325 else : # V or F grid 191 326 yrange = 90. -lat_reg 192 327 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)) 194 329 195 330 if math == xr : lat1D.attrs = lat.attrs … … 206 341 207 342 def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) : 343 math = __math__ (ptab) 208 344 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() 211 347 except : pass 212 348 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) 218 353 219 354 return tab … … 231 366 232 367 ''' 233 if tab.shape[-1] == 1 : 234 extend = tab 368 math = __math__ (tab) 369 370 if tab.shape[-1] == 1 : extend = tab 235 371 236 372 else : 237 373 if jpi == None : jpi = tab.shape[-1] 238 239 try :240 if type (tab) == xr.core.dataarray.DataArray : math = xr241 else : math = np242 except : math = np243 374 244 375 if Lon : xplus = -360.0 … … 248 379 extend = tab 249 380 else : 381 if nperio == 0 or nperio == 4.2 : 382 istart = 0 ; le=jpi+1 ; la=0 250 383 if nperio == 1 : 251 384 istart = 0 ; le=jpi+1 ; la=0 252 if nperio > 1 : # OPA case with to halo points for periodicity385 if nperio == 4 or nperio == 6 : # OPA case with two halo points for periodicity 253 386 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) 260 400 return extend 261 401 … … 286 426 'time_counter', 'time', 'tbnds', 287 427 '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) 289 429 290 430 ff = ff.transpose (*coord_order) 291 431 return ff 292 432 293 def lbc_init (ptab, nperio=None , cd_type='T') :294 """433 def lbc_init (ptab, nperio=None) : 434 ''' 295 435 Prepare for all lbc calls 296 436 … … 305 445 306 446 See NEMO documentation for further details 307 """447 ''' 308 448 jpj, jpi = ptab.shape[-2:] 309 449 if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio) 310 450 311 451 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) 320 453 sys.exit () 321 454 322 return jpj, jpi, nperio , cd_type455 return jpj, jpi, nperio 323 456 324 457 325 458 def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) : 326 """459 ''' 327 460 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]) 330 462 nperio : Type of periodicity 331 463 cd_type : Grid specification : T, U, V or F … … 333 465 334 466 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) 337 469 psgn = ptab.dtype.type (psgn) 338 470 math = __math__ (ptab) 471 472 if math == xr : ztab = ptab.values.copy () 473 else : ztab = ptab.copy () 339 474 # 340 475 #> East-West boundary conditions … … 342 477 if nperio in [1, 4, 6] : 343 478 # ... cyclic 344 ptab [..., :, 0] = ptab [..., :, -2]345 ptab [..., :, -1] = ptab [..., :, 1]479 ztab [..., :, 0] = ztab [..., :, -2] 480 ztab [..., :, -1] = ztab [..., :, 1] 346 481 # 347 482 #> North-South boundary conditions … … 349 484 if nperio in [3, 4] : # North fold T-point pivot 350 485 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 ] 354 489 355 490 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 ] 359 494 360 495 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 ] 363 498 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] 365 500 366 501 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 ] 370 505 371 506 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 ] 376 511 377 512 if nperio in [4.2] : # North fold T-point pivot 378 513 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 ] 380 515 381 516 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] 383 518 384 519 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 ] 386 521 387 522 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 391 525 if nperio in [5, 6] : # North fold F-point pivot 392 526 if cd_type in ['T', 'W'] : 393 ptab [..., -1, 0: ] = psgn * ptab [..., -2, -1::-1 ]527 ztab [..., -1, 0: ] = psgn * ztab [..., -2, -1::-1 ] 394 528 395 529 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 ? 398 532 399 533 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 ] 402 536 403 537 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 424 542 # 425 543 #> East-West boundary conditions … … 427 545 if nperio in [1, 4, 6] : 428 546 # ... 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 555 def 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 431 575 432 576 # … … 434 578 # -------------------------------- 435 579 if nperio in [1, 3, 4, 5, 6] : 436 ptab [...,0,:] = sval580 ztab [..., 0, :] = sval 437 581 438 582 # … … 441 585 if nperio in [3, 4] : # North fold T-point pivot 442 586 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 446 589 447 590 if cd_type == 'U' : 448 ptab [..., -1, :] = sval449 ptab [..., -2, jpi//2: ] = sval591 ztab [..., -1, : ] = sval 592 ztab [..., -2, jpi//2+1: ] = sval 450 593 451 594 if cd_type == 'V' : 452 ptab [..., -2, : ] = sval453 ptab [..., -1, : ] = sval595 ztab [..., -2, : ] = sval 596 ztab [..., -1, : ] = sval 454 597 455 598 if cd_type == 'F' : 456 ptab [..., -2, :] = sval457 ptab [..., -1, :] = sval599 ztab [..., -2, : ] = sval 600 ztab [..., -1, : ] = sval 458 601 459 602 if nperio in [4.2] : # North fold T-point pivot 460 603 if cd_type in [ 'T', 'W' ] : # T-, W-point 461 ptab [..., -1, jpi//2: ] = sval604 ztab [..., -1, jpi//2 : ] = sval 462 605 463 606 if cd_type == 'U' : 464 ptab [..., -1, jpi//2-1:-1] = sval607 ztab [..., -1, jpi//2-1:-1] = sval 465 608 466 609 if cd_type == 'V' : 467 ptab [..., -1, 1: ] = sval610 ztab [..., -1, 1: ] = sval 468 611 469 612 if cd_type == 'F' : 470 ptab [..., -1, 0:-1 ] = sval613 ztab [..., -1, 0:-1 ] = sval 471 614 472 615 if nperio in [5, 6] : # North fold F-point pivot 473 616 if cd_type in ['T', 'W'] : 474 ptab [..., -1, 0: ] = sval617 ztab [..., -1, 0: ] = sval 475 618 476 619 if cd_type == 'U' : 477 ptab [..., -1, 0:-1 ] = sval478 ptab [..., -1, -1 ] = sval620 ztab [..., -1, 0:-1 ] = sval 621 ztab [..., -1, -1 ] = sval 479 622 480 623 if cd_type == 'V' : 481 ptab [..., -1, 0: ] = sval482 ptab [..., -2, jpi//2: ] = sval624 ztab [..., -1, 0: ] = sval 625 ztab [..., -2, jpi//2: ] = sval 483 626 484 627 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 491 633 492 634 def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) : 493 """635 ''' 494 636 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] 497 638 nperio : Type of periodicity 498 639 cd_type : Grid specification : T, U, V or F … … 500 641 501 642 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) 505 646 psgn = ptab.dtype.type (psgn) 506 647 ztab = ptab.copy () 507 648 # 508 649 #> East-West boundary conditions … … 510 651 if nperio in [1, 4, 6] : 511 652 # ... 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 514 659 515 660 # … … 518 663 if nperio in [3, 4] : # North fold T-point pivot 519 664 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 524 668 if cd_type == 'U' : 525 ptab [..., -1, : ] = sval669 ztab [..., -1, : ] = sval 526 670 527 671 if cd_type == 'V' : 528 ptab [..., -2, : ] = sval529 ptab [..., -1, : ] = sval672 ztab [..., -2, : ] = sval 673 ztab [..., -1, : ] = sval 530 674 531 675 if cd_type == 'F' : 532 ptab [..., -2, : ] = sval533 ptab [..., -1, : ] = sval676 ztab [..., -2, : ] = sval 677 ztab [..., -1, : ] = sval 534 678 535 679 if nperio in [4.2] : # North fold T-point pivot 536 680 if cd_type in [ 'T', 'W' ] : # T-, W-point 537 ptab [..., -1, jpi//2: ] = sval681 ztab [..., -1, jpi//2: ] = sval 538 682 539 683 if cd_type == 'U' : 540 ptab [..., -1, jpi//2-1:-1] = sval684 ztab [..., -1, jpi//2-1:-1] = sval 541 685 542 686 if cd_type == 'V' : 543 ptab [..., -1, 1: ] = sval687 ztab [..., -1, 1: ] = sval 544 688 545 689 if cd_type == 'F' : 546 ptab [..., -1, 0:-1 ] = sval690 ztab [..., -1, 0:-1 ] = sval 547 691 548 692 if nperio in [5, 6] : # North fold F-point pivot 549 693 if cd_type in ['T', 'W'] : 550 ptab [..., -1, : ] = sval694 ztab [..., -1, : ] = sval 551 695 552 696 if cd_type == 'U' : 553 ptab [..., -1, : ] = sval697 ztab [..., -1, : ] = sval 554 698 555 699 if cd_type == 'V' : 556 ptab [..., -1, :] = sval557 ptab [..., -2, jpi//2: ] = sval700 ztab [..., -1, : ] = sval 701 ztab [..., -2, jpi//2: ] = sval 558 702 559 703 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 709 def 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 563 747 return ptab 564 748 565 def geo2oce (pxx, pyy, pzz, glam, gphi) : 749 def 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 769 def 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 853 def geo2en (pxx, pyy, pzz, glam, gphi) : 566 854 ''' 567 855 Change vector from geocentric to east/north 568 856 569 Input 570 pxx, pyy, pzz : components on the geoce ,tric system571 glam, gphi : longitu e and latitude of the points572 573 ''' 857 Inputs : 858 pxx, pyy, pzz : components on the geocentric system 859 glam, gphi : longitude and latitude of the points 860 ''' 861 574 862 gsinlon = np.sin (rad * glam) 575 863 gcoslon = np.cos (rad * glam) … … 582 870 return pte, ptn 583 871 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 ##---------------------------------------------------------------------- 872 def 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 595 881 gsinlon = np.sin (rad * glam) 596 882 gcoslon = np.cos (rad * glam) … … 603 889 604 890 return pxx, pyy, pzz 891 892 def 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 943 def 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 955 def 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 1058 def 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 1090 def 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 1105 def 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 1119 def 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 1139 def 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 1159 def 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 1175 def 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 1191 def 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 1199 def 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 1215 def 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 1231 def 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 1247 def 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 1263 def 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 1271 def 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 1279 def 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 1294 def 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 1309 def 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 1332 def 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 1356 def 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 1405 def 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 605 1428 606 1429 ## =========================================================================== … … 609 1432 ## 610 1433 ## =========================================================================== 1434 1435 def 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 6 6 ### =========================================================================== 7 7 ## 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 ## 16 28 ## 17 29 ## SVN information
Note: See TracChangeset
for help on using the changeset viewer.