[6064] | 1 | # -*- coding: utf-8 -*- |
---|
| 2 | ## =========================================================================== |
---|
| 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. |
---|
| 12 | ## |
---|
| 13 | ''' |
---|
| 14 | Utilities to plot NEMO ORCA fields |
---|
| 15 | Periodicity and other stuff |
---|
| 16 | |
---|
| 17 | olivier.marti@lsce.ipsl.fr |
---|
| 18 | ''' |
---|
| 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,5,6] |
---|
| 27 | |
---|
| 28 | ## SVN information |
---|
| 29 | __Author__ = "$Author$" |
---|
| 30 | __Date__ = "$Date$" |
---|
| 31 | __Revision__ = "$Revision$" |
---|
| 32 | __Id__ = "$Id$" |
---|
| 33 | __HeadURL = "$HeadURL$" |
---|
| 34 | |
---|
| 35 | def __guessNperio__ (jpj, jpi, nperio=None) : |
---|
| 36 | ''' |
---|
| 37 | Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details) |
---|
| 38 | |
---|
| 39 | Inputs |
---|
| 40 | jpi : number of longitudes |
---|
| 41 | nperio : periodicity parameter |
---|
| 42 | ''' |
---|
| 43 | if nperio == None : |
---|
| 44 | if jpi == 182 : nperio = 4 # ORCA2. We choose legacy orca2. |
---|
| 45 | if jpi == 362 : nperio = 6 # eORCA1. |
---|
| 46 | if jpi == 1442 : nperio = 6 # ORCA025. |
---|
| 47 | |
---|
| 48 | if jpj == 149 : nperio = 4 # ORCA2. We choose legacy orca2. |
---|
| 49 | if jpj == 332 : nperio = 6 # eORCA1. |
---|
| 50 | # |
---|
| 51 | if nperio == None : |
---|
| 52 | sys.exit ('in nemo module : nperio not found, and cannot by guessed' ) |
---|
| 53 | else : |
---|
| 54 | print ('nperio set as {:d} (deduced from jpi={:d})'.format (nperio, jpi)) |
---|
| 55 | |
---|
| 56 | return nperio |
---|
| 57 | |
---|
| 58 | def __guessPoint__ (ptab) : |
---|
| 59 | ''' |
---|
| 60 | Tries to guess the grid point (periodicity parameter. See NEMO documentation for details) |
---|
| 61 | |
---|
| 62 | For array conforments with xgcm requirements |
---|
| 63 | |
---|
| 64 | Inputs |
---|
| 65 | ptab : xarray array |
---|
| 66 | |
---|
| 67 | Credits : who is the original author ? |
---|
| 68 | ''' |
---|
| 69 | gP = None |
---|
| 70 | if isinstance (ptab, xr.core.dataarray.DataArray) : |
---|
| 71 | if 'x_c' in ptab.dims and 'y_c' in ptab.dims : gP = 'T' |
---|
| 72 | if 'x_f' in ptab.dims and 'y_c' in ptab.dims : gP = 'U' |
---|
| 73 | if 'x_c' in ptab.dims and 'y_f' in ptab.dims : gP = 'V' |
---|
| 74 | if 'x_f' in ptab.dims and 'y_f' in ptab.dims : gP = 'F' |
---|
| 75 | if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_c' in ptab.dims : gP = 'T' |
---|
| 76 | if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'W' |
---|
| 77 | if 'x_f' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'U' |
---|
| 78 | if 'x_c' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'V' |
---|
| 79 | if 'x_f' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'F' |
---|
| 80 | |
---|
| 81 | if gP == None : |
---|
| 82 | sys.exit ('in nemo module : cd_type not found, and cannot by guessed' ) |
---|
| 83 | else : |
---|
| 84 | print ('Grid set as', gP, 'deduced from dims ', ptab.dims) |
---|
| 85 | return gP |
---|
| 86 | else : |
---|
| 87 | sys.exit ('in nemo module : cd_type not found, input is not an xarray data' ) |
---|
| 88 | #return gP |
---|
| 89 | |
---|
| 90 | def fixed_lon (lon) : |
---|
| 91 | ''' |
---|
| 92 | Returns corrected longitudes for nicer plots |
---|
| 93 | |
---|
| 94 | fixed_lon (lon) |
---|
| 95 | lon : longitudes of the grid. At least 2D. |
---|
| 96 | |
---|
| 97 | Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7 |
---|
| 98 | ''' |
---|
| 99 | fixed_lon = lon.copy () |
---|
| 100 | for i, start in enumerate (np.argmax (np.abs (np.diff (lon, axis=-1)) > 180, axis=-1)) : |
---|
| 101 | fixed_lon[..., i, start+1:] += 360 |
---|
| 102 | |
---|
| 103 | # Special case for eORCA025 |
---|
| 104 | if lon.shape[-1] == 1442 : lon [..., -2, :] = lon [..., -3, :] |
---|
| 105 | |
---|
| 106 | return fixed_lon |
---|
| 107 | |
---|
| 108 | def jeq (lat) : |
---|
| 109 | ''' |
---|
| 110 | Returns j index of equator in the grid |
---|
| 111 | |
---|
| 112 | lat : latitudes of the grid. At least 2D. |
---|
| 113 | ''' |
---|
| 114 | jeq = np.argmin (np.abs (np.float64(lat[:,0]))) |
---|
| 115 | return jeq |
---|
| 116 | |
---|
| 117 | def lon1D (lon, lat=None) : |
---|
| 118 | ''' |
---|
| 119 | Returns 1D longitude for simple plots. |
---|
| 120 | |
---|
| 121 | lon : longitudes of the grid |
---|
| 122 | lat (optionnal) : latitudes of the grid |
---|
| 123 | ''' |
---|
| 124 | if np.max (lat) != None : |
---|
| 125 | je = jeq (lat) |
---|
| 126 | lon1D = lon[..., je, :] |
---|
| 127 | else : |
---|
| 128 | jpj, jpi = lon.shape[-2:] |
---|
| 129 | lon1D = lon[..., jpj//3, :] |
---|
| 130 | |
---|
| 131 | start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1) |
---|
| 132 | lon1D[..., start+1:] += 360 |
---|
| 133 | |
---|
| 134 | return lon1D |
---|
| 135 | |
---|
| 136 | def latreg (lat, diff=0.1) : |
---|
| 137 | ''' |
---|
| 138 | Returns maximum j index where gridlines are along latitude in the northern hemisphere |
---|
| 139 | |
---|
| 140 | lat : latitudes of the grid (2D) |
---|
| 141 | diff [optional] : tolerance |
---|
| 142 | ''' |
---|
| 143 | if diff == None : |
---|
| 144 | dy = np.float64 (np.mean (np.abs (lat - np.roll(lat,shift=1,axis=-2)))) |
---|
| 145 | diff = dy/100. |
---|
| 146 | |
---|
| 147 | je = jeq (lat) |
---|
| 148 | jreg = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< diff)[-1][-1] + je |
---|
| 149 | latreg = np.float64 (lat[...,jreg,:].mean(axis=-1)) |
---|
| 150 | JREG = jreg |
---|
| 151 | |
---|
| 152 | return jreg, latreg |
---|
| 153 | |
---|
| 154 | def lat1D (lat) : |
---|
| 155 | ''' |
---|
| 156 | Returns 1D latitudes for zonal means and simple plots. |
---|
| 157 | |
---|
| 158 | lat : latitudes of the grid (2D) |
---|
| 159 | ''' |
---|
| 160 | jpj, jpi = lat.shape[-2:] |
---|
| 161 | try : |
---|
| 162 | if type (lat) == xr.core.dataarray.DataArray : math = xr |
---|
| 163 | else : math = np |
---|
| 164 | except : math = np |
---|
| 165 | |
---|
| 166 | dy = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2)))) |
---|
| 167 | je = jeq (lat) |
---|
| 168 | lat_eq = np.float64 (lat[...,je,:].mean(axis=-1)) |
---|
| 169 | |
---|
| 170 | jreg, lat_reg = latreg (lat) |
---|
| 171 | lat_ave = np.mean (lat, axis=-1) |
---|
| 172 | |
---|
| 173 | if ( np.abs (lat_eq) < dy/100. ) : # T, U or W grid |
---|
| 174 | dys = (90.-lat_reg) / (jpj-jreg-1)*0.5 |
---|
| 175 | yrange = 90.-dys-lat_reg |
---|
| 176 | else : # V or F grid |
---|
| 177 | yrange = 90. -lat_reg |
---|
| 178 | |
---|
| 179 | lat1D = math.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) ) |
---|
| 180 | |
---|
| 181 | if math == xr : lat1D.attrs = lat.attrs |
---|
| 182 | |
---|
| 183 | return lat1D |
---|
| 184 | |
---|
| 185 | def latlon1D (lat, lon) : |
---|
| 186 | ''' |
---|
| 187 | Returns simple latitude and longitude (1D) for simple plots. |
---|
| 188 | |
---|
| 189 | lat, lon : latitudes and longitudes of the grid (2D) |
---|
| 190 | ''' |
---|
| 191 | return lat1D (lat), lon1D (lon, lat) |
---|
| 192 | |
---|
| 193 | def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) : |
---|
| 194 | try : |
---|
| 195 | lon = lon.to_masked_array() |
---|
| 196 | lat = lat.to_masked_array() |
---|
| 197 | except : pass |
---|
| 198 | |
---|
| 199 | mask = np.logical_and ( np.logical_and(lat>y0, lat<y1), |
---|
| 200 | np.logical_or (np.logical_or (np.logical_and(lon>x0, lon<x1), np.logical_and(lon+360>x0, lon+360<x1)), |
---|
| 201 | np.logical_and(lon-360>x0, lon-360<x1)) ) |
---|
| 202 | try : tab = xr.where (mask, ptab, np.nan) |
---|
| 203 | except : tab = np.where (mask, ptab, np.nan) |
---|
| 204 | |
---|
| 205 | return tab |
---|
| 206 | |
---|
| 207 | def extend (tab, Lon=False, jplus=25, jpi=None, nperio=4) : |
---|
| 208 | ''' |
---|
| 209 | Returns extended field eastward to have better plots, and box average crossing the boundary |
---|
| 210 | Works only for xarray and numpy data (?) |
---|
| 211 | |
---|
| 212 | tab : field to extend. |
---|
| 213 | Lon : (optional, default=False) : if True, add 360 in the extended parts of the field |
---|
| 214 | jpi : normal longitude dimension of the field. exrtend does nothing it the actual |
---|
| 215 | size of the field != jpi (avoid to extend several times) |
---|
| 216 | jplus (optional, default=25) : number of points added on the east side of the field |
---|
| 217 | |
---|
| 218 | ''' |
---|
| 219 | if tab.shape[-1] == 1 : |
---|
| 220 | extend = tab |
---|
| 221 | |
---|
| 222 | else : |
---|
| 223 | if jpi == None : jpi = tab.shape[-1] |
---|
| 224 | |
---|
| 225 | try : |
---|
| 226 | if type (tab) == xr.core.dataarray.DataArray : math = xr |
---|
| 227 | else : math = np |
---|
| 228 | except : math = np |
---|
| 229 | |
---|
| 230 | if Lon : xplus = -360.0 |
---|
| 231 | else : xplus = 0.0 |
---|
| 232 | |
---|
| 233 | if tab.shape[-1] > jpi : |
---|
| 234 | extend = tab |
---|
| 235 | else : |
---|
| 236 | if nperio == 1 : |
---|
| 237 | istart = 0 ; le=jpi+1 ; la=0 |
---|
| 238 | if nperio > 1 : # OPA case with to halo points for periodicity |
---|
| 239 | istart = 1 ; le=jpi-2 ; la=1 # Perfect, except at the pole that should be masked by lbc_plot |
---|
| 240 | |
---|
| 241 | if math == xr : |
---|
| 242 | extend = xr.concat ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), dim=tab.dims[-1] ) |
---|
| 243 | if math == np : |
---|
| 244 | extend = np.concatenate ( (tab[..., istart:istart+le+1] + xplus, tab[..., istart+la:jplus]), axis=-1 ) |
---|
| 245 | |
---|
| 246 | return extend |
---|
| 247 | |
---|
| 248 | def orca2reg (ff, lat_name='nav_lat', lon_name='nav_lon', y_name='y', x_name='x') : |
---|
| 249 | ''' |
---|
| 250 | Assign ORCA dataset on a regular grid. |
---|
| 251 | For use in the tropical region. |
---|
| 252 | |
---|
| 253 | Inputs : |
---|
| 254 | ff : xarray dataset |
---|
| 255 | lat_name, lon_name : name of latitude and longitude 2D field in ff |
---|
| 256 | y_name, x_name : namex of dimensions in ff |
---|
| 257 | |
---|
| 258 | Returns : xarray dataset with rectangular grid. Incorrect above 20°N |
---|
| 259 | ''' |
---|
| 260 | # Compute 1D longitude and latitude |
---|
| 261 | (lat, lon) = latlon1D (ff[lat_name], ff[lon_name]) |
---|
| 262 | # Assign lon and lat as dimensions of the dataset |
---|
| 263 | if y_name in ff.dims : |
---|
| 264 | lat = xr.DataArray (lat, coords=[lat,], dims=['lat',]) |
---|
| 265 | ff = ff.rename_dims ({y_name: "lat",}).assign_coords (lat=lat) |
---|
| 266 | if x_name in ff.dims : |
---|
| 267 | lon = xr.DataArray (lon, coords=[lon,], dims=['lon',]) |
---|
| 268 | ff = ff.rename_dims ({x_name: "lon",}).assign_coords (lon=lon) |
---|
| 269 | # Force dimensions to be in the right order |
---|
| 270 | coord_order = ['lat', 'lon'] |
---|
| 271 | for dim in [ 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z', |
---|
| 272 | 'time_counter', 'time', 'tbnds', |
---|
| 273 | 'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] : |
---|
| 274 | if dim in ff.dims : coord_order.insert (0, dim ) |
---|
| 275 | |
---|
| 276 | ff = ff.transpose (*coord_order) |
---|
| 277 | return ff |
---|
| 278 | |
---|
| 279 | def lbc_init (ptab, nperio=None, cd_type='T') : |
---|
| 280 | """ |
---|
| 281 | Prepare for all lbc calls |
---|
| 282 | |
---|
| 283 | Set periodicity on input field |
---|
| 284 | nperio : Type of periodicity |
---|
| 285 | 0 : No periodicity |
---|
| 286 | 1, 4, 6 : Cyclic on i dimension (generaly longitudes) with 2 points halo |
---|
| 287 | 2 : Obsolete (was symmetric condition at southern boundary ?) |
---|
| 288 | 3, 4 : North fold T-point pivot (legacy ORCA2) |
---|
| 289 | 5, 6 : North fold F-point pivot (ORCA1, ORCA025, ORCA2 with new grid for paleo) |
---|
| 290 | cd_type : Grid specification : T, U, V or F |
---|
| 291 | |
---|
| 292 | See NEMO documentation for further details |
---|
| 293 | """ |
---|
| 294 | jpj, jpi = ptab.shape[-2:] |
---|
| 295 | if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio) |
---|
| 296 | |
---|
| 297 | if nperio not in nperio_valid_range : |
---|
| 298 | print ( 'nperio=', nperio, ' is not in the valid range', nperio_valid_range ) |
---|
| 299 | sys.exit() |
---|
| 300 | |
---|
| 301 | if cd_type == None and isinstance (ptab, xr.core.dataarray.DataArray) : |
---|
| 302 | cd_type = __guessPoint__ (ptab) |
---|
| 303 | |
---|
| 304 | if cd_type not in ['T', 'U', 'V', 'F', 'W', 'F' ]: |
---|
| 305 | print ( 'cd_type=', cd_type, ' is not in the list T, U, V, F, W, F' ) |
---|
| 306 | sys.exit () |
---|
| 307 | |
---|
| 308 | return jpj, jpi, nperio, cd_type |
---|
| 309 | |
---|
| 310 | |
---|
| 311 | def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) : |
---|
| 312 | """ |
---|
| 313 | Set periodicity on input field |
---|
| 314 | ptab : Input array (works |
---|
| 315 | rank 2 at least : ptab[...., lat, lon] |
---|
| 316 | nperio : Type of periodicity |
---|
| 317 | cd_type : Grid specification : T, U, V or F |
---|
| 318 | psgn : For change of sign for vector components (1 for scalars, -1 for vector components) |
---|
| 319 | |
---|
| 320 | See NEMO documentation for further details |
---|
| 321 | """ |
---|
| 322 | jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) |
---|
| 323 | psgn = ptab.dtype.type (psgn) |
---|
| 324 | |
---|
| 325 | # |
---|
| 326 | #> East-West boundary conditions |
---|
| 327 | # ------------------------------ |
---|
| 328 | if nperio in [1, 4, 6] : |
---|
| 329 | # ... cyclic |
---|
| 330 | ptab [..., :, 0] = ptab [..., :, -2] |
---|
| 331 | ptab [..., :, -1] = ptab [..., :, 1] |
---|
| 332 | # |
---|
| 333 | #> North-South boundary conditions |
---|
| 334 | # -------------------------------- |
---|
| 335 | if nperio in [3, 4] : # North fold T-point pivot |
---|
| 336 | if cd_type in [ 'T', 'W' ] : # T-, W-point |
---|
| 337 | ptab [..., -1, 1: ] = psgn * ptab [..., -3, -1:0:-1] |
---|
| 338 | ptab [..., -1, 0 ] = psgn * ptab [..., -3, 2 ] |
---|
| 339 | ptab [..., -2, jpi//2: ] = psgn * ptab [..., -2, jpi//2:0:-1 ] |
---|
| 340 | |
---|
| 341 | if cd_type == 'U' : |
---|
| 342 | ptab [..., -1, 0:-1 ] = psgn * ptab [..., -3, -1:0:-1 ] |
---|
| 343 | ptab [..., -1, 0 ] = psgn * ptab [..., -3, 1 ] |
---|
| 344 | ptab [..., -1, -1 ] = psgn * ptab [..., -3, -2 ] |
---|
| 345 | |
---|
| 346 | if nemo_4U_bug : |
---|
| 347 | ptab [..., -2, jpi//2+1:-1] = psgn * ptab [..., -2, jpi//2-2:0:-1] |
---|
| 348 | ptab [..., -2, jpi//2-1 ] = psgn * ptab [..., -2, jpi//2 ] |
---|
| 349 | else : |
---|
| 350 | ptab [..., -2, jpi//2-1:-1] = psgn * ptab [..., -2, jpi//2:0:-1] |
---|
| 351 | |
---|
| 352 | if cd_type == 'V' : |
---|
| 353 | ptab [..., -2, 1: ] = psgn * ptab [..., -3, jpi-1:0:-1 ] |
---|
| 354 | ptab [..., -1, 1: ] = psgn * ptab [..., -4, -1:0:-1 ] |
---|
| 355 | ptab [..., -1, 0 ] = psgn * ptab [..., -4, 2 ] |
---|
| 356 | |
---|
| 357 | if cd_type == 'F' : |
---|
| 358 | ptab [..., -2, 0:-1 ] = psgn * ptab [..., -3, -1:0:-1 ] |
---|
| 359 | ptab [..., -1, 0:-1 ] = psgn * ptab [..., -4, -1:0:-1 ] |
---|
| 360 | ptab [..., -1, 0 ] = psgn * ptab [..., -4, 1 ] |
---|
| 361 | ptab [..., -1, -1 ] = psgn * ptab [..., -4, -2 ] |
---|
| 362 | |
---|
| 363 | if nperio in [5, 6] : # North fold F-point pivot |
---|
| 364 | if cd_type in ['T', 'W'] : |
---|
| 365 | ptab [..., -1, 0: ] = psgn * ptab [..., -2, -1::-1 ] |
---|
| 366 | |
---|
| 367 | if cd_type == 'U' : |
---|
| 368 | ptab [..., -1, 0:-1 ] = psgn * ptab [..., -2, -2::-1 ] |
---|
| 369 | ptab [..., -1, -1 ] = psgn * ptab [..., -2, 0 ] # Bug ? |
---|
| 370 | |
---|
| 371 | if cd_type == 'V' : |
---|
| 372 | ptab [..., -1, 0: ] = psgn * ptab [..., -3, -1::-1 ] |
---|
| 373 | ptab [..., -2, jpi//2: ] = psgn * ptab [..., -2, jpi//2-1::-1 ] |
---|
| 374 | |
---|
| 375 | if cd_type == 'F' : |
---|
| 376 | ptab [..., -1, 0:-1 ] = psgn * ptab [..., -3, -2::-1 ] |
---|
| 377 | ptab [..., -1, -1 ] = psgn * ptab [..., -3, 0 ] |
---|
| 378 | ptab [..., -2, jpi//2:-1] = psgn * ptab [..., -2, jpi//2-2::-1 ] |
---|
| 379 | |
---|
| 380 | return ptab |
---|
| 381 | |
---|
| 382 | def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) : |
---|
| 383 | # |
---|
| 384 | """ |
---|
| 385 | Mask fields on duplicated points |
---|
| 386 | ptab : Input array |
---|
| 387 | rank 2 at least : ptab[...., lat, lon] |
---|
| 388 | nperio : Type of periodicity |
---|
| 389 | cd_type : Grid specification : T, U, V or F |
---|
| 390 | |
---|
| 391 | See NEMO documentation for further details |
---|
| 392 | """ |
---|
| 393 | |
---|
| 394 | jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) |
---|
| 395 | |
---|
| 396 | # |
---|
| 397 | #> East-West boundary conditions |
---|
| 398 | # ------------------------------ |
---|
| 399 | if nperio in [1, 4, 6] : |
---|
| 400 | # ... cyclic |
---|
| 401 | ptab [...,:, 0] = sval |
---|
| 402 | ptab [...,:, -1] = sval |
---|
| 403 | |
---|
| 404 | # |
---|
| 405 | #> South (in which nperio cases ?) |
---|
| 406 | # -------------------------------- |
---|
| 407 | if nperio in [1, 3, 4, 5, 6] : |
---|
| 408 | ptab [...,0,:] = sval |
---|
| 409 | |
---|
| 410 | # |
---|
| 411 | #> North-South boundary conditions |
---|
| 412 | # -------------------------------- |
---|
| 413 | if nperio in [3, 4] : # North fold T-point pivot |
---|
| 414 | if cd_type in [ 'T', 'W' ] : # T-, W-point |
---|
| 415 | ptab [..., -1, 1: ] = sval |
---|
| 416 | ptab [..., -1, 0 ] = sval |
---|
| 417 | ptab [..., -2, jpi//2: ] = sval |
---|
| 418 | |
---|
| 419 | if cd_type == 'U' : |
---|
| 420 | ptab [..., -1, : ] = sval |
---|
| 421 | ptab [..., -2, jpi//2: ] = sval |
---|
| 422 | |
---|
| 423 | if cd_type == 'V' : |
---|
| 424 | ptab [..., -2, : ] = sval |
---|
| 425 | ptab [..., -1, : ] = sval |
---|
| 426 | |
---|
| 427 | if cd_type == 'F' : |
---|
| 428 | ptab [..., -2, : ] = sval |
---|
| 429 | ptab [..., -1, : ] = sval |
---|
| 430 | |
---|
| 431 | if nperio in [5, 6] : # North fold F-point pivot |
---|
| 432 | if cd_type in ['T', 'W'] : |
---|
| 433 | ptab [..., -1, 0: ] = sval |
---|
| 434 | |
---|
| 435 | if cd_type == 'U' : |
---|
| 436 | ptab [..., -1, 0:-1 ] = sval |
---|
| 437 | ptab [..., -1, -1 ] = sval |
---|
| 438 | |
---|
| 439 | if cd_type == 'V' : |
---|
| 440 | ptab [..., -1, 0: ] = sval |
---|
| 441 | ptab [..., -2, jpi//2: ] = sval |
---|
| 442 | |
---|
| 443 | if cd_type == 'F' : |
---|
| 444 | ptab [..., -1, 0:-1 ] = sval |
---|
| 445 | ptab [..., -1, -1 ] = sval |
---|
| 446 | ptab [..., -2, jpi//2:-1] = sval |
---|
| 447 | |
---|
| 448 | return ptab |
---|
| 449 | |
---|
| 450 | |
---|
| 451 | def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) : |
---|
| 452 | """ |
---|
| 453 | Set periodicity on input field, adapted for plotting for any cartopy projection |
---|
| 454 | ptab : Input array (works |
---|
| 455 | rank 2 at least : ptab[...., lat, lon] |
---|
| 456 | nperio : Type of periodicity |
---|
| 457 | cd_type : Grid specification : T, U, V or F |
---|
| 458 | psgn : For change of sign for vector components (1 for scalars, -1 for vector components) |
---|
| 459 | |
---|
| 460 | See NEMO documentation for further details |
---|
| 461 | """ |
---|
| 462 | |
---|
| 463 | jpj, jpi, nperio, cd_type = lbc_init (ptab, nperio, cd_type) |
---|
| 464 | psgn = ptab.dtype.type (psgn) |
---|
| 465 | |
---|
| 466 | # |
---|
| 467 | #> East-West boundary conditions |
---|
| 468 | # ------------------------------ |
---|
| 469 | if nperio in [1, 4, 6] : |
---|
| 470 | # ... cyclic |
---|
| 471 | ptab [..., :, 0] = ptab [..., :, -2] |
---|
| 472 | ptab [..., :, -1] = ptab [..., :, 1] |
---|
| 473 | |
---|
| 474 | # |
---|
| 475 | #> North-South boundary conditions |
---|
| 476 | # -------------------------------- |
---|
| 477 | if nperio in [3, 4] : # North fold T-point pivot |
---|
| 478 | if cd_type in [ 'T', 'W' ] : # T-, W-point |
---|
| 479 | ptab [..., -1, :] = sval |
---|
| 480 | ptab [..., -2, :jpi//2] = sval |
---|
| 481 | #ptab [..., -2, :] = sval |
---|
| 482 | |
---|
| 483 | if cd_type == 'U' : |
---|
| 484 | ptab [..., -1, : ] = sval |
---|
| 485 | |
---|
| 486 | if cd_type == 'V' : |
---|
| 487 | ptab [..., -2, : ] = sval |
---|
| 488 | ptab [..., -1, : ] = sval |
---|
| 489 | |
---|
| 490 | if cd_type == 'F' : |
---|
| 491 | ptab [..., -2, : ] = sval |
---|
| 492 | ptab [..., -1, : ] = sval |
---|
| 493 | |
---|
| 494 | if nperio in [5, 6] : # North fold F-point pivot |
---|
| 495 | if cd_type in ['T', 'W'] : |
---|
| 496 | ptab [..., -1, : ] = sval |
---|
| 497 | |
---|
| 498 | if cd_type == 'U' : |
---|
| 499 | ptab [..., -1, : ] = sval |
---|
| 500 | |
---|
| 501 | if cd_type == 'V' : |
---|
| 502 | ptab [..., -1, : ] = sval |
---|
| 503 | ptab [..., -2, jpi//2: ] = sval |
---|
| 504 | |
---|
| 505 | if cd_type == 'F' : |
---|
| 506 | ptab [..., -1, : ] = sval |
---|
| 507 | ptab [..., -2, jpi//2:-1] = sval |
---|
| 508 | |
---|
| 509 | return ptab |
---|
| 510 | |
---|
| 511 | def geo2oce (pxx, pyy, pzz, glam, gphi) : |
---|
| 512 | ''' |
---|
| 513 | Change vector from geocentric to east/north |
---|
| 514 | |
---|
| 515 | Input |
---|
| 516 | pxx, pyy, pzz : components on the geoce,tric system |
---|
| 517 | glam, gphi : longitue and latitude of the points |
---|
| 518 | |
---|
| 519 | ''' |
---|
| 520 | gsinlon = np.sin (rad * glam) |
---|
| 521 | gcoslon = np.cos (rad * glam) |
---|
| 522 | gsinlat = np.sin (rad * gphi) |
---|
| 523 | gcoslat = np.cos (rad * gphi) |
---|
| 524 | |
---|
| 525 | pte = - pxx * gsinlon + pyy * gcoslon |
---|
| 526 | ptn = - pxx * gcoslon * gsinlat - pyy * gsinlon * gsinlat + pzz * gcoslat |
---|
| 527 | |
---|
| 528 | return pte, ptn |
---|
| 529 | |
---|
| 530 | def oce2geo (pte, ptn, glam, gphi) : |
---|
| 531 | ##---------------------------------------------------------------------- |
---|
| 532 | ## *** ROUTINE oce2geo *** |
---|
| 533 | ## |
---|
| 534 | ## ** Purpose : |
---|
| 535 | ## |
---|
| 536 | ## ** Method : Change vector from east/north to geocentric |
---|
| 537 | ## |
---|
| 538 | ## History : |
---|
| 539 | ## # (A. Caubel) oce2geo - Original code |
---|
| 540 | ##---------------------------------------------------------------------- |
---|
| 541 | gsinlon = np.sin (rad * glam) |
---|
| 542 | gcoslon = np.cos (rad * glam) |
---|
| 543 | gsinlat = np.sin (rad * gphi) |
---|
| 544 | gcoslat = np.cos (rad * gphi) |
---|
| 545 | |
---|
| 546 | pxx = - pte * gsinlon - ptn * gcoslon * gsinlat |
---|
| 547 | pyy = pte * gcoslon - ptn * gsinlon * gsinlat |
---|
| 548 | pzz = ptn * gcoslat |
---|
| 549 | |
---|
| 550 | return pxx, pyy, pzz |
---|
| 551 | |
---|
| 552 | ## =========================================================================== |
---|
| 553 | ## |
---|
| 554 | ## That's all folk's !!! |
---|
| 555 | ## |
---|
| 556 | ## =========================================================================== |
---|