- Timestamp:
- 10/25/23 17:11:15 (6 months ago)
- Location:
- TOOLS/MOSAIX
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
TOOLS/MOSAIX/CalvingWeights.py
r6190 r6666 25 25 ## software by incorrectly or partially configured 26 26 ## 27 ## 28 import numpy as np, xarray as xr 29 import sys, os, platform, argparse, textwrap, time 27 '''Python code to generates weights for calving. Launched by `CreateWeights.bash`. 28 29 More information with python `CalvingWeights.py -h`. 30 31 - The atmosphere model integrate the calving over on specific 32 regions. Presently the regions are three latitudinal bands with 33 limits 90°S/40°S/50°N/90°N. 34 35 - The weights distribute the calving in the ocean in the 36 corresponding latitudinal bands. By default, the distribution is 37 uniform. 38 39 - For the southernmost band, it is possible to use a geographical 40 distribution read in a file. These files are currently available 41 for eORCA1 and eORCA025. 42 43 ## SVN information 44 Author = "$Author$" 45 Date = "$Date$" 46 Revision = "$Revision$" 47 Id = "$Id$" 48 HeadURL = "$HeadURL$" 49 ''' 50 import sys 51 import os 52 import getpass 53 import platform 54 import argparse 55 import time 56 import numpy as np 57 import xarray as xr 30 58 from scipy import ndimage 31 59 import nemo 32 60 33 61 ## SVN information 34 __Author__ = "$Author$" 35 __Date__ = "$Date$" 36 __Revision__ = "$Revision$" 37 __Id__ = "$Id$" 38 __HeadURL__ = "$HeadURL$" 39 __SVN_Date__ = "$SVN_Date: $" 40 62 __SVN__ = ( { 63 'Author' : "$Author$", 64 'Date' : '$Date$', 65 'Revision' : '$Revision$', 66 'Id' : '$Id$', 67 'HeadURL' : '$HeadURL$', 68 'SVN_Date' : '$SVN_Date: $', 69 } ) 41 70 ### 42 71 43 ### ===== Handling command line parameters ========================== ========================72 ### ===== Handling command line parameters ========================== 44 73 # Creating a parser 45 74 parser = argparse.ArgumentParser ( 46 description = """Compute calving weights""",75 description = '''Compute calving weights''', 47 76 epilog='-------- End of the help message --------') 48 77 49 78 # Adding arguments 50 79 parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', 51 choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 52 parser.add_argument ('--atm' , help='atm model name (ICO* or LMD*)', type=str, default='ICO40' ) 53 parser.add_argument ('--type' , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full' ) 80 choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 81 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 82 parser.add_argument ('--atm' , type=str, default='ICO40', 83 help='atm model name (ICO* or LMD*)' ) 84 parser.add_argument ('--type' , help='Type of distribution', type=str, 85 choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full' ) 54 86 ## parser.add_argument ('--dir' , help='Directory of input file', type=str, default='./' ) 55 parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' ) 56 parser.add_argument ('--repartition_var' , help='Variable name for iceshelf' , type=str, default=None) 87 parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, 88 default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' ) 89 parser.add_argument ('--repartition_var' , type=str, default=None, 90 help='Variable name for iceshelf' ) 57 91 parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) 58 92 parser.add_argument ('--areas' , help='masks file name', default='areas.nc' ) 59 93 parser.add_argument ('--masks' , help='areas file name', default='masks.nc' ) 60 94 parser.add_argument ('--o2a' , help='o2a file name' , default='o2a.nc' ) 61 parser.add_argument ('--output' , help='output rmp file name', default='rmp_tlmd_to_torc_calving.nc' ) 62 parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', 63 choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 64 parser.add_argument ('--ocePerio', help='periodicity of ocean grid', type=float, default=0, choices=nemo.nperio_valid_range ) 95 parser.add_argument ('--output', help='output rmp file name', 96 default='rmp_tlmd_to_torc_calving.nc' ) 97 parser.add_argument ('--fmt' , default='netcdf4' , 98 help='NetCDF file format, using nco syntax', 99 choices=['classic', 'netcdf3', '64bit', '64bit_data', 100 '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 101 parser.add_argument ('--ocePerio', type=float, default=0, choices=nemo.NPERIO_VALID_RANGE, 102 help='periodicity of ocean grid' ) 103 parser.add_argument ('--modelName' , help='Name of model', type=str, default=None) 65 104 66 105 # Parse command line 67 myargs = parser.parse_args ()106 myargs = parser.parse_args () 68 107 69 108 # Model Names 70 src_Name = myargs.atm ; dst_Name = myargs.oce 71 109 src_Name = myargs.atm 110 dst_Name = myargs.oce 111 model_name = myargs.modelName 112 72 113 # Default vars 73 if myargs.repartition_var ==None :114 if myargs.repartition_var is None : 74 115 # Runoff from Antarctica iceshelves (Depoorter, 2013) 75 if myargs.type == 'iceshelf' : repartitionVar = 'sornfisf'76 116 if myargs.type == 'iceshelf' : 117 repartitionVar = 'sornfisf' 77 118 # Runoff from Antarctica iceberg (Depoorter, 2013) 78 if myargs.type == 'iceberg' : repartitionVar = 'Icb_Flux' 79 80 else : repartitionVar = myargs.repartition_var 119 if myargs.type == 'iceberg' : 120 repartitionVar = 'Icb_Flux' 121 else : 122 repartitionVar = myargs.repartition_var 81 123 82 124 # Latitude limits of each calving zone … … 92 134 if myargs.fmt == 'netcdf4' : FmtNetcdf = 'NETCDF4' 93 135 if myargs.fmt == 'netcdf4_classic' : FmtNetcdf = 'NETCDF4_CLASSIC' 94 ### ===================================================================== =====================136 ### ===================================================================== 95 137 96 138 # Model short names 97 src_name = None ; dst_name = None 139 src_name = None 140 dst_name = None 98 141 if src_Name.count('ICO') != 0 : src_name = 'ico' 99 142 if src_Name.count('LMD') != 0 : src_name = 'lmd' … … 138 181 if dst_Name in ['eORCA1.2', 'eORCA1.3'] : nperio_dst = 6 139 182 if dst_Name in ['ORCA025', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] : nperio_dst = 6 140 183 141 184 print ('nperio_dst: ' + str(nperio_dst) ) 142 185 dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T', sval=0 ) … … 161 204 # Set Hudson Strait to 0 to fill Hudson Bay 162 205 dst_mskutil[997:1012,907] = 0 163 206 164 207 if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4', 'eORCA1.4.2' ] : 165 208 print ('Filling some seas for ' + dst_Name) … … 174 217 # Set Hudson Strait to 0 to fill Hudson Bay 175 218 dst_mskutil[284,222:225] = 0 176 219 177 220 if dst_Name in [ 'ORCA2.3', 'ORCA2.4' ] : 178 221 print ('Filling some seas for ' + dst_Name) … … 193 236 # Fill closed seas with image processing library 194 237 # ---------------------------------------------- 195 dst_mskutil = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(dst_mskutil, nperio=nperio_dst, cd_type='T')), nperio=nperio_dst, cd_type='T', sval=0) 238 dst_mskutil = nemo.lbc_mask ( 1-ndimage.binary_fill_holes ( 239 1-nemo.lbc(dst_mskutil, nperio=nperio_dst, cd_type='T')), 240 nperio=nperio_dst, cd_type='T', sval=0) 196 241 197 242 # Surfaces … … 208 253 209 254 # Indices 210 ( src_jpj, src_jpi) = src_lat.shape ; src_grid_size = src_jpj*src_jpi 211 ( dst_jpj, dst_jpi) = dst_lat.shape ; dst_grid_size = dst_jpj*dst_jpi 255 ( src_jpj, src_jpi) = src_lat.shape 256 src_grid_size = src_jpj*src_jpi 257 ( dst_jpj, dst_jpi) = dst_lat.shape 258 dst_grid_size = dst_jpj*dst_jpi 212 259 orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 # Fortran indices (starting at one) 213 260 214 ### ===== Reading needed data ======================================== ==========215 if myargs.type in ['iceberg', 'iceshelf' ]: 261 ### ===== Reading needed data ======================================== 262 if myargs.type in ['iceberg', 'iceshelf' ]: 216 263 # Reading data file for calving or iceberg geometry around Antarctica 217 264 print ( 'Opening ' + myargs.repartition_file) … … 226 273 print (' ') 227 274 228 ### ===== Starting loop on basins===================================== =========275 ### ===== Starting loop on basins===================================== 229 276 230 277 # Initialise some fields … … 238 285 ## Loop 239 286 for n_bas in np.arange ( nb_zone ) : 240 south = False ; ok_run = False 241 lat_south = np.min(limit_lat[n_bas]) ; lat_north = np.max(limit_lat[n_bas]) 287 south = False 288 ok_run = False 289 lat_south = np.min(limit_lat[n_bas]) 290 lat_north = np.max(limit_lat[n_bas]) 242 291 if lat_south <= -60.0 : south = True 243 244 print ( 'basin: {:2d} -- Latitudes: {:+.1f} {:+.1f} --'.format(n_bas, lat_south, lat_north))292 293 print ( f'basin: {n_bas:2d} -- Latitudes: {lat_south:+.1f} {lat_north:+.1f} --' ) 245 294 ## 246 if myargs.type == 'iceberg' and south : ok_run = True ; print ('Applying iceberg to south' ) 247 if myargs.type == 'iceshelf' and south : ok_run = True ; print ('Applying iceshelf to south' ) 248 if myargs.type == 'iceberg' and not south : ok_run = False ; print ('Skipping iceberg: not south ') 249 if myargs.type == 'iceshelf' and not south : ok_run = False ; print ('Skipping iceshelf: not south ') 250 if myargs.type == 'nosouth' and south : ok_run = False ; print ('Skipping south: nosouth case' ) 251 if myargs.type == 'nosouth' and not south : ok_run = True ; print ('Running: not in south, uniform repartition') 252 if myargs.type == 'full' : ok_run = True ; print ('Running general trivial case, uniform repartition' ) 295 if myargs.type == 'iceberg' and south : 296 ok_run = True 297 print ('Applying iceberg to south' ) 298 if myargs.type == 'iceshelf' and south : 299 ok_run = True 300 print ('Applying iceshelf to south' ) 301 if myargs.type == 'iceberg' and not south : 302 ok_run = False 303 print ('Skipping iceberg: not south ') 304 if myargs.type == 'iceshelf' and not south : 305 ok_run = False 306 print ('Skipping iceshelf: not south ') 307 if myargs.type == 'nosouth' and south : 308 ok_run = False 309 print ('Skipping south: nosouth case' ) 310 if myargs.type == 'nosouth' and not south : 311 ok_run = True 312 print ('Running: not in south, uniform repartition') 313 if myargs.type == 'full' : 314 ok_run = True 315 print ('Running general trivial case, uniform repartition' ) 253 316 254 317 if ok_run : 255 318 # Calving integral send to one point per latitude bands 256 319 index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1 257 320 258 321 # Basin mask 259 basin_msk[n_bas] = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 ) 322 basin_msk[n_bas] = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), 323 dst_mskutil, 0 ) 260 324 261 325 # Repartition pattern 262 if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition[n_bas] = repartition * basin_msk[n_bas] 263 else : key_repartition[n_bas] = basin_msk[n_bas] 326 if myargs.type in ['iceberg', 'iceshelf' ] : 327 key_repartition[n_bas] = repartition * basin_msk[n_bas] 328 else : 329 key_repartition[n_bas] = basin_msk[n_bas] 264 330 265 331 # Integral and normalisation 266 332 sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil ).values 267 333 key_repartition = key_repartition / sum_repartition 268 334 269 335 print ( 'Sum of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] )) ) 270 336 print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] * dst_srfutil ).values) ) 271 337 272 338 # Basin surface (modulated by repartition key) 273 339 basin_srfutil = np.sum ( key_repartition[n_bas] * dst_srfutil ) 274 340 275 341 # Weights and links 276 342 poids = 1.0 / basin_srfutil 277 matrix_local = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids.values, 0. ) 343 matrix_local = np.where ( basin_msk[n_bas].ravel() > 0.5, 344 key_repartition[n_bas].ravel()*poids.values, 0. ) 278 345 matrix_local = matrix_local[matrix_local > 0.0] # Keep only non zero values 279 346 num_links = np.count_nonzero (matrix_local) … … 289 356 # 290 357 #print ( 'Sum of remap_matrix : {:12.3e}'.format(np.sum(matrix_local)) ) 291 print ( 'Point in atmospheric grid : { :4d} -- num_links: {:6d}'.format(index_src, num_links) )358 print ( 'Point in atmospheric grid : {index_src:4d} -- num_links: {num_links:6d}' ) 292 359 print (' ') 293 294 295 360 296 361 ## End of loop 297 362 print (' ') 298 363 299 # 364 # 300 365 num_links = np.count_nonzero (remap_matrix) 301 print ( 'Total num_links : {:10d}'.format(num_links))366 print ( f'Total num_links : {num_links:10d}' ) 302 367 303 368 # Diag : interpolate uniform field 304 src_field = np.zeros(shape=(src_grid_size ))369 src_field = np.zeros(shape=(src_grid_size,)) 305 370 for n_bas in np.arange ( nb_zone ) : 306 371 index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1 … … 311 376 dst_field [dst_address [link]-1] += remap_matrix[link] * src_field[src_address[link]-1] 312 377 313 ### ===== Writing the weights file, for OASIS MCT ======================= ===================378 ### ===== Writing the weights file, for OASIS MCT ======================= 314 379 315 380 # Output file 316 381 calving = myargs.output 317 382 318 383 print ('Output file: ' + calving ) 319 384 320 remap_matrix = xr.DataArray (np.reshape(remap_matrix, (num_links, 1)), dims = ['num_links', 'num_wgts'] ) 321 src_address = xr.DataArray (src_address.astype(np.int32) , dims = ['num_links',] ) 322 src_address.attrs['convention'] = "Fortran style addressing, starting at 1" 323 dst_address = xr.DataArray (dst_address.astype(np.int32) , dims = ['num_links',] ) 324 dst_address.attrs['convention'] = "Fortran style addressing, starting at 1" 325 326 src_grid_dims = xr.DataArray ( np.array (( src_jpi, src_jpi ), dtype=np.int32), dims = ['src_grid_rank',] ) 385 remap_matrix = xr.DataArray (np.reshape(remap_matrix, (num_links, 1)), 386 dims = ['num_links', 'num_wgts'] ) 387 src_address = xr.DataArray (src_address.astype(np.int32) , 388 dims = ['num_links',] ) 389 src_address.attrs['convention'] = 'Fortran style addressing, starting at 1' 390 dst_address = xr.DataArray (dst_address.astype(np.int32) , 391 dims = ['num_links',] ) 392 dst_address.attrs['convention'] = 'Fortran style addressing, starting at 1' 393 394 src_grid_dims = xr.DataArray ( np.array (( src_jpi, src_jpi ), dtype=np.int32), 395 dims = ['src_grid_rank',] ) 327 396 src_grid_center_lon = xr.DataArray ( src_lon.values.ravel(), dims = ['src_grid_size',] ) 328 397 src_grid_center_lat = xr.DataArray ( src_lat.values.ravel(), dims = ['src_grid_size',] ) 329 src_grid_center_lon.attrs['units']='degrees_east' ; src_grid_center_lon.attrs['long_name']='Longitude' 330 src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" 331 src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' 332 src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" 333 src_grid_corner_lon = xr.DataArray ( src_clo.values.reshape( (src_jpi*src_jpj, src_clo.shape[0]) ), dims = ['src_grid_size', 'src_grid_corners'] ) 334 src_grid_corner_lat = xr.DataArray ( src_cla.values.reshape( (src_jpi*src_jpj, src_clo.shape[0]) ), dims = ['src_grid_size', 'src_grid_corners' ] ) 335 src_grid_corner_lon.attrs['units']="degrees_east" 336 src_grid_corner_lat.attrs['units']="degrees_north" 398 src_grid_center_lon.attrs.update ( {'units':'degrees_east ', 'long_name':'Longitude', 399 'standard_name':'longitude', 400 'bounds':'src_grid_corner_lon'}) 401 src_grid_center_lat.attrs.update ( {'units':'degrees_north', 'long_name':'Latitude' , 402 'standard_name':'latitude' , 403 'bounds':'src_grid_corner_lat'}) 404 src_grid_corner_lon = xr.DataArray ( src_clo.values.reshape( (src_jpi*src_jpj, src_clo.shape[0]) ), 405 dims = ['src_grid_size', 'src_grid_corners'] ) 406 src_grid_corner_lat = xr.DataArray ( src_cla.values.reshape( (src_jpi*src_jpj, src_clo.shape[0]) ), 407 dims = ['src_grid_size', 'src_grid_corners' ] ) 408 src_grid_corner_lon.attrs['units']='degrees_east' 409 src_grid_corner_lat.attrs['units']='degrees_north' 337 410 src_grid_area = xr.DataArray ( src_srf.values.ravel(), dims = ['src_grid_size',] ) 338 src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" 339 src_grid_imask = xr.DataArray ( src_msk.values.ravel().astype(np.int32) , dims = ['src_grid_size',] ) 340 src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" 411 src_grid_area.attrs.update ({'long_name':'Grid area', 'standard_name':'cell_area', 'units':'m2'}) 412 src_grid_imask = xr.DataArray ( src_msk.values.ravel().astype(np.int32) , 413 dims = ['src_grid_size',] ) 414 src_grid_imask.attrs.update ( {'long_name':'-sea mask', 'units':'Land:1, Ocean:0' }) 341 415 342 416 # -- 343 dst_grid_dims = xr.DataArray ( np.array(( dst_jpi, dst_jpi), dtype=np.int32) , dims = ['dst_grid_rank', ]) 417 dst_grid_dims = xr.DataArray ( np.array(( dst_jpi, dst_jpi), dtype=np.int32) , 418 dims = ['dst_grid_rank', ]) 344 419 dst_grid_center_lon = xr.DataArray ( dst_lon.values.ravel(), dims = ['dst_grid_size', ]) 345 420 dst_grid_center_lat = xr.DataArray ( dst_lat.values.ravel(), dims = ['dst_grid_size', ]) 346 dst_grid_center_lon.attrs['units']='degrees_east' ; dst_grid_center_lon.attrs['long_name']='Longitude' 347 dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" 348 dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' 349 dst_grid_center_lat.attrs['long_name']='latitude' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" 350 dst_grid_corner_lon = xr.DataArray ( dst_clo.values.reshape(dst_jpi*dst_jpj, dst_clo.shape[0] ), dims = ['dst_grid_size', 'dst_grid_corners'] ) 351 dst_grid_corner_lat = xr.DataArray ( dst_cla.values.reshape(dst_jpi*dst_jpj, dst_clo.shape[0] ), dims = ['dst_grid_size', 'dst_grid_corners'] ) 352 dst_grid_corner_lon.attrs['units']="degrees_east" 353 dst_grid_corner_lat.attrs['units']="degrees_north" 421 dst_grid_center_lon.attrs.update ( {'units':'degrees_east ', 'long_name':'Longitude', 422 'standard_name':'longitude', 423 'bounds':'dst_grid_corner_lon'}) 424 dst_grid_center_lat.attrs.update ( {'units':'degrees_north', 'long_name':'Latitude' , 425 'standard_name':'latitude' , 426 'bounds':'dst_grid_corner_lat'}) 427 dst_grid_corner_lon = xr.DataArray ( dst_clo.values.reshape(dst_jpi*dst_jpj, dst_clo.shape[0] ), 428 dims = ['dst_grid_size', 'dst_grid_corners'] ) 429 dst_grid_corner_lat = xr.DataArray ( dst_cla.values.reshape(dst_jpi*dst_jpj, dst_clo.shape[0] ), 430 dims = ['dst_grid_size', 'dst_grid_corners'] ) 431 dst_grid_corner_lon.attrs['units']='degrees_east' 432 dst_grid_corner_lat.attrs['units']='degrees_north' 354 433 dst_grid_area = xr.DataArray ( dst_srf.values.ravel(), dims = ['dst_grid_size', ] ) 355 dst_grid_area.attrs ['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2"434 dst_grid_area.attrs.update ({'long_name':'Grid area', 'standard_name':'cell_area', 'units':'m2'}) 356 435 dst_grid_imask = xr.DataArray ( dst_msk.values.ravel(), dims = ['dst_grid_size',] ) 357 dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" 358 359 dst_bassin = xr.DataArray ( basin_msk.reshape ( (nb_zone,dst_grid_size) ) , dims = ['nb_zone', 'dst_grid_size', ] ) 360 dst_repartition = xr.DataArray ( key_repartition.reshape( (nb_zone,dst_grid_size) ) , dims = ['nb_zone', 'dst_grid_size', ] ) 361 362 dst_southLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , dims = ['nb_zone', ] ) 363 dst_northLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , dims = ['nb_zone', ] ) 436 dst_grid_imask.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0' } ) 437 438 dst_bassin = xr.DataArray ( basin_msk.reshape ( (nb_zone,dst_grid_size) ) , 439 dims = ['nb_zone', 'dst_grid_size', ] ) 440 dst_repartition = xr.DataArray ( key_repartition.reshape( (nb_zone,dst_grid_size) ) , 441 dims = ['nb_zone', 'dst_grid_size', ] ) 442 443 dst_southLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , 444 dims = ['nb_zone', ] ) 445 dst_northLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , 446 dims = ['nb_zone', ] ) 364 447 365 448 # Additionnal fields for debugging 366 449 # ================================ 367 dst_grid_imaskutil = xr.DataArray ( dst_mskutil.ravel().astype(np.float32), dims = ['dst_grid_size',] ) 368 dst_grid_imaskutil.attrs['long_name']="Land-sea mask" ; dst_grid_imaskutil.attrs['units']="Land:1, Ocean:0" 450 dst_grid_imaskutil = xr.DataArray ( dst_mskutil.ravel().astype(np.float32), 451 dims = ['dst_grid_size',] ) 452 dst_grid_imaskutil.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0' }) 369 453 dst_grid_imaskclose = xr.DataArray ( dst_closed.values.ravel(), dims = ['dst_grid_size',] ) 370 dst_grid_imaskclose.attrs['long_name']="Land-sea mask" ; dst_grid_imaskclose.attrs['units']="Land:1, Ocean:0" 371 372 dst_lon_addressed = xr.DataArray ( dst_lon.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 373 dst_lat_addressed = xr.DataArray ( dst_lat.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 374 dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" 375 dst_lat_addressed.attrs['long_name']="Latitude" ; dst_lat_addressed.attrs['standard_name']="latitude" ; dst_lat_addressed.attrs['units']="degrees_north" 376 dst_area_addressed = xr.DataArray ( dst_srf.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 377 dst_area_addressed.attrs['long_name']="Grid area" ; dst_area_addressed.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 378 dst_imask_addressed = xr.DataArray ( dst_msk.values.ravel()[dst_address-1].ravel(), dims = ['num_links', ] ) 379 dst_imask_addressed.attrs['long_name']="Land-sea mask" ; dst_imask_addressed.attrs['units']="Land:1, Ocean:0" 380 dst_imaskutil_addressed = xr.DataArray ( dst_mskutil.ravel()[dst_address-1].astype(np.float32), dims = ['num_links'] ) 381 dst_imaskutil_addressed.attrs['long_name']="Land-sea mask" ; dst_imaskutil_addressed.attrs['units']="Land:1, Ocean:0" 454 dst_grid_imaskclose.attrs.update ({'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0' }) 455 456 dst_lon_addressed = xr.DataArray ( dst_lon.values.ravel()[dst_address-1].ravel(), 457 dims = ['num_links',] ) 458 dst_lat_addressed = xr.DataArray ( dst_lat.values.ravel()[dst_address-1].ravel(), 459 dims = ['num_links',] ) 460 dst_lon_addressed.attrs.update ( {'units':'degrees_east ', 'long_name':'Longitude', 461 'standard_name':'longitude' }) 462 dst_lat_addressed.attrs.update ( {'units':'degrees_north', 'long_name':'Latitude' , 463 'standard_name':'latitude' }) 464 dst_area_addressed = xr.DataArray ( dst_srf.values.ravel()[dst_address-1].ravel(), 465 dims = ['num_links',] ) 466 dst_area_addressed.attrs.update ( {'long_name':'Grid area0', 'standard_name':'cell_area', 467 'units':'m2' }) 468 dst_imask_addressed = xr.DataArray ( dst_msk.values.ravel()[dst_address-1].ravel(), 469 dims = ['num_links', ] ) 470 dst_imask_addressed.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0' }) 471 dst_imaskutil_addressed = xr.DataArray ( dst_mskutil.ravel()[dst_address-1].astype(np.float32), 472 dims = ['num_links'] ) 473 dst_imaskutil_addressed.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0' }) 382 474 383 475 src_field = xr.DataArray ( src_field, dims = ['src_grid_size',] ) 384 476 dst_field = xr.DataArray ( dst_field, dims = ['dst_grid_size',] ) 385 477 386 f_calving = xr.Dataset ( { 478 f_calving = xr.Dataset ( { 387 479 'remap_matrix' : remap_matrix, 388 480 'src_address' : src_address, … … 417 509 } ) 418 510 419 f_calving.attrs['Conventions'] = "CF-1.6" 420 f_calving.attrs['source'] = "IPSL Earth system model" 421 f_calving.attrs['group'] = "ICMC IPSL Climate Modelling Center" 422 f_calving.attrs['Institution'] = "IPSL https.//www.ipsl.fr" 423 f_calving.attrs['Ocean'] = dst_Name + " https://www.nemo-ocean.eu" 424 f_calving.attrs['Atmosphere'] = src_Name + " http://lmdz.lmd.jussieu.fr" 425 if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.attrs['originalFiles'] = myargs.repartition_file 426 f_calving.attrs['associatedFiles'] = grids + " " + areas + " " + masks 427 428 f_calving.attrs['description'] = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX" 429 f_calving.attrs['title'] = calving 430 f_calving.attrs['Program'] = "Generated by " + sys.argv[0] + " with flags " + ' '.join (sys.argv[1:]) 431 432 f_calving.attrs['repartitionType'] = myargs.type 511 f_calving.attrs.update ( { 512 'Conventions' : 'CF-1.6', 513 'source' : 'IPSL Earth system model', 514 'group' : 'ICMC IPSL Climate Modelling Center', 515 'Institution' : 'IPSL https.//www.ipsl.fr', 516 'Ocean' : f'{dst_Name} https://www.nemo-ocean.eu', 517 'Atmosphere' : f'{src_Name} http://lmdz.lmd.jussieu.fr', 518 'associatedFiles' : f'{grids} {areas} {masks}', 519 'description' : 'Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX', 520 'title' : calving, 521 'Program' : f'Generated by {sys.argv[0]} with flags ' + ' '.join (sys.argv[1:]), 522 'repartitionType' : myargs.type, 523 'gridsFile' : grids, 524 'areasFile' : areas, 525 'masksFile' : masks, 526 'timeStamp' : time.asctime(), 527 'directory' : os.getcwd (), 528 'HOSTNAME' : platform.node (), 529 'LOGNAME' : getpass.getuser(), 530 'Python' : 'Python version ' + platform.python_version (), 531 'OS' : platform.system (), 532 'release' : platform.release (), 533 'hardware' : platform.machine (), 534 'conventions' : 'SCRIP', 535 'dest_grid' : 'curvilinear', 536 'Model' : model_name, 537 'SVN_Author' : '$Author$', 538 'SVN_Date' : '$Date$', 539 'SVN_Revision' : '$Revision$', 540 'SVN_Id' : '$Id$', 541 'SVN_HeadURL' : '$HeadURL$', 542 }) 543 433 544 if myargs.type in [ 'iceberg', 'iceshelf' ] : 434 f_calving.attrs['repartitionFile'] = myargs.repartition_file 435 f_calving.attrs['repartitionVar'] = repartitionVar 436 f_calving.attrs['gridsFile'] = grids 437 f_calving.attrs['areasFile'] = areas 438 f_calving.attrs['masksFile'] = masks 439 f_calving.attrs['timeStamp'] = time.asctime() 440 try : f_calving.attrs['directory'] = os.getcwd () 441 except : pass 442 try : f_runoff.attrs['HOSTNAME'] = platform.node () 443 except : pass 444 try : f_runoff.attrs['LOGNAME'] = os.getlogin () 445 except : pass 446 try : f_runoff.attrs['Python'] = "Python version " + platform.python_version () 447 except : pass 448 try : f_runoff.attrs['OS'] = platform.system () 449 except : pass 450 try : f_runoff.attrs['release'] = platform.release () 451 except : pass 452 try : f_runoff.attrs['hardware'] = platform.machine () 453 except : pass 454 f_calving.attrs['conventions'] = "SCRIP" 455 if src_name == 'lmd' : f_calving.attrs['source_grid'] = "curvilinear" 456 if src_name == 'ico' : f_calving.attrs['source_grid'] = "unstructured" 457 f_calving.attrs['dest_grid'] = "curvilinear" 458 f_calving.attrs['Model'] = "IPSL CM6" 459 f_calving.attrs['SVN_Author'] = "$Author$" 460 f_calving.attrs['SVN_Date'] = "$Date$" 461 f_calving.attrs['SVN_Revision'] = "$Revision$" 462 f_calving.attrs['SVN_Id'] = "$Id$" 463 f_calving.attrs['SVN_HeadURL'] = "$HeadURL$" 545 f_calving.attrs.update ( { 546 'originalFiles' : myargs.repartition_file, 547 'repartitionFile' : myargs.repartition_file, 548 'repartitionVar' : repartitionVar, 549 } ) 464 550 465 551 ## … … 467 553 f_calving.close() 468 554 469 print ( 'That''s all folks !')555 print ("That''s all folks !") 470 556 ## ====================================================================================== -
TOOLS/MOSAIX/CreateWeightsMask.bash
r6360 r6666 3 3 #MSUB -o Out_WeightsMask # Standard output 4 4 #MSUB -e Out_WeightsMask # Error output 5 #MSUB -n 1# Number of processors6 #MSUB -T 7200 # Time limit (seconds)7 #MSUB -Q normal5 #MSUB -n 4 # Number of processors 6 #MSUB -T 1800 # Time limit (seconds) 7 #MSUB -Q test 8 8 #MSUB -q rome 9 9 #MSUB -p devcmip6 … … 21 21 ### =========================================================================== 22 22 ## 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 ## 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. 43 31 ## 44 32 ## SVN information … … 50 38 # 51 39 ## Tested with : 52 # CplModel=eORCA1. 2xLMD144142 ; qsub -r ${CplModel} -o Out_${CplModel} -e Out_${CplModel} CreateWeightsMask.bash40 # CplModel=eORCA1.4.2xICO60 ; ccc_msub -r ${CplModel} -o Out_${CplModel} -e Out_${CplModel} CreateWeightsMask.bash 53 41 # CplModel=ORCA2.3xLMD9695 54 42 # CplModel=ORCA2.3xICO30 … … 59 47 # CplModel=eORCA025.1xLMD144142 60 48 # CplModel=eORCA025.1xLMD256256 61 # 49 # CplModel=eORCA1.4.2xICO40 50 # CplModel=eORCA1.4.2xICO60 51 # ccc_msub -r ${CplModel} -o Out_${CplModel} -e Out_${CplModel} CreateWeightsMask.bash 52 53 echo 'Test' 54 ls -alt /ccc/work/cont003/gencmip6/p86mart/TMP/MOSAIX 62 55 63 56 set +vx … … 82 75 echo ${Titre}"Defines model"${Norm} 83 76 # ================================= 84 #CplModel=ORCA2.3xLMD969577 CplModel=ORCA2.3xLMD9695 85 78 #CplModel=ORCA2.3xICO30 86 79 #CplModel=ORCA2.3xICO40 … … 90 83 #CplModel=eORCA1.2xLMD256256 91 84 #CplModel=eORCA1.2xICO40 85 #CplModel=eORCA1.2xICO60 92 86 #CplModel=eORCA1.4.2xICO40 93 87 #CplModel=eORCA1.2xICO450 94 88 #CplModel=eORCA025.1xLMD256256 95 CplModel=eORCA1.4.2xICO8096 89 97 90 #Version="v0" ; Comment="Fully tested in IPSLCM6 eORCA1.2 x LMD 144x142" … … 144 137 # \!/ No spaces in any analysis \!/ 145 138 # 146 # Specific commands : 'Runoff', 'Calving' 139 # Specific commands : 'Runoff', 'Calving', 'Grids' 147 140 # 148 141 # Keywords : … … 194 187 #CommandList=( Grids ) 195 188 189 if [[ ${Version} = test_runoff_* ]] ; then 190 CommandList=( Runoff ) 191 fi 192 196 193 ## =========================================================================== 197 194 ## … … 286 283 set +e 287 284 R_IN=$(ccc_home -u igcmg --cccwork)/IGCM 288 TMPDIR=${CCCWORKDIR}/TMP 285 TMPDIR=${CCCWORKDIR}/TMP/MOSAIX 289 286 SUBMIT_DIR=${BRIDGE_MSUB_PWD:-${SUBMIT_DIR}} 290 287 MpiRun="time ccc_mprun" … … 292 289 module purge 293 290 source ${SUBMIT_DIR}/arch.env 294 module load nco #/4.9.1 295 module load cdo #/1.9.5 291 module load nco 296 292 module load python3 297 293 module load datadir/igcmg 298 294 module list 299 cp ${SUBMIT_DIR}/arch.env . ${TMPDIR}300 295 set -e 301 296 ;; … … 311 306 312 307 set -x ; set -e 313 308 ls -alt ${TMPDIR} 314 309 mkdir -p ${TMPDIR}/${CplModel} || exit 1 315 310 cd ${TMPDIR}/${CplModel} || exit 2 … … 346 341 cp ${SUBMIT_DIR}/iodef_atm_to_oce.xml . 347 342 cp ${SUBMIT_DIR}/iodef_oce_to_atm.xml . 343 cp ${SUBMIT_DIR}/arch.env . 348 344 349 345 [[ -f ${R_IN}/OCE/NEMO/${OCE}/${OCE}_coordinates_mask.nc ]] && cp ${R_IN}/OCE/NEMO/${OCE}/${OCE}_coordinates_mask.nc . … … 604 600 --attribute Ocean,global,o,c,"${OCE} https://www.nemo-ocean.eu" \ 605 601 --attribute Atmosphere,global,o,c,"${ATM} http://lmdz.lmd.jussieu.fr" \ 606 --attribute production,global,o,c,"$(finger ${LOGNAME} | head -1 | awk '{print $4}')" \607 602 --attribute originalFiles,global,o,c,"${OCE}_coordinates_mask.nc ${ATM}_grid_mask.nc" \ 608 603 --attribute associatedFiles,global,o,c,"grids_${CplModel}.nc areas_${CplModel}.nc masks_${CplModel}.nc" \ … … 822 817 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --output=rmp_t${atm}_to_t${oce}_runoff_${runOff_atmQuantity}_to_${runOff_oceQuantity}.nc \ 823 818 --fmt=${FMT_XIOS} \ 824 --atmQuantity=${runOff_atmQuantity} --oceQuantity=${runOff_oceQuantity} --ocePerio=${OcePerio} 819 --atmQuantity=${runOff_atmQuantity} --oceQuantity=${runOff_oceQuantity} --ocePerio=${OcePerio} --modelName="${ModelName}" 825 820 fi 826 821 … … 835 830 --oce=${OCE} --atm=${ATM} --type=nosouth \ 836 831 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 837 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 832 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --modelName="${ModelName}" 838 833 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceberg.nc --fmt=${FMT_XIOS} \ 839 834 --oce=${OCE} --atm=${ATM} --type=iceberg --repartition_file=eORCA_R025_runoff_v1.2.nc --repartition_var=Icb_flux \ 840 835 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 841 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 836 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --modelName="${ModelName}" 842 837 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceshelf.nc --fmt=${FMT_XIOS} \ 843 838 --oce=${OCE} --atm=${ATM} --type=iceshelf --repartition_file=eORCA_R025_runoff_v1.2.nc --repartition_var=sornfisf \ 844 839 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 845 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 840 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --modelName="${ModelName}" 846 841 ;; 847 842 … … 858 853 --oce=${OCE} --atm=${ATM} --type=nosouth \ 859 854 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 860 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 855 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --modelName="${ModelName}" 861 856 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceberg.nc --fmt=${FMT_XIOS} \ 862 857 --oce=${OCE} --atm=${ATM} --type=iceberg --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=Icb_flux \ 863 858 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 864 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 859 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --modelName="${ModelName}" 865 860 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceshelf.nc --fmt=${FMT_XIOS} \ 866 861 --oce=${OCE} --atm=${ATM} --type=iceshelf --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=sornfisf \ 867 862 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 868 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 863 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --modelName="${ModelName}" 869 864 ;; 870 865 … … 873 868 --oce=${OCE} --atm=${ATM} --type=full --ocePerio=${OcePerio} \ 874 869 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 875 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 870 --o2a=${ATM}_grid_maskFrom_${OCE}.nc --modelName="${ModelName}" 876 871 ;; 877 872 esac … … 957 952 Python version : ${PYTHON_VER} 958 953 959 EOF 960 961 echo 'SVN Information : ' >> README.txt 962 echo '$Author$ ' >> README.txt 963 echo '$Date$ ' >> README.txt 964 echo '$Revision$ ' >> README.txt 965 echo '$Id$ ' >> README.txt 966 echo '$HeadURL$ ' >> README.txt 954 SVN Information : 955 $Author$ 956 $Date$ 957 $Revision$ 958 $Id$ 959 $HeadURL$ 960 961 EOF 962 963 967 964 968 965 echo ${Titre}"Compute checksums and add them to README"${Norm} … … 996 993 ## 997 994 ## =========================================================================== 998 -
TOOLS/MOSAIX/RunoffWeights.py
r6314 r6666 29 29 ## 30 30 ## 31 '''Python code to generates weights for run-off. Launched by `CreateWeights.bash`. 32 33 More information with python `RunOffWeights.py -h`. 34 35 - Defines a coastal band on land in the atmosphere model. For LMDZ 36 rectilinear grids, the width of this band is parametrized, in grid 37 points. For ico grids, the band contains only point with a 38 fraction of ocean in ]0,1[. 39 40 - Defines a coastal band on ocean in the ocean model. The width of 41 this band is parametrized. 42 43 - An atmosphere point of the coastal band is linked to an ocean 44 point in the coastal band if the distance between the two is less 45 than searchRadius. 46 47 ## SVN information 48 Author = "$Author$" 49 Date = "$Date$" 50 Revision = "$Revision$" 51 Id = "$Id$" 52 HeadURL = "$HeadURL$" 53 ''' 54 55 ## Modules 56 import sys 57 import os 58 import platform 59 import getpass 60 import argparse 61 import time 62 import numpy as np 63 import xarray as xr 64 from scipy import ndimage 65 import nemo 66 31 67 # SVN information 32 __Author__ = "$Author$" 33 __Date__ = "$Date$" 34 __Revision__ = "$Revision$" 35 __Id__ = "$Id$" 36 __HeadURL__ = "$HeadURL$" 37 __SVN_Date__ = "$SVN_Date: $" 38 ## 39 40 ## Modules 41 import numpy as np, xarray as xr 42 import nemo 43 from scipy import ndimage 44 import sys, os, platform, argparse, textwrap, time 68 __SVN__ = ({ 69 'Author' : '$Author$', 70 'Date' : '$Date$', 71 'Revision' : '$Revision$', 72 'Id' : '$Id$', 73 'HeadURL' : '$HeadURL$', 74 'SVN_Date' : '$SVN_Date: $', 75 }) 76 ## 45 77 46 78 ## Useful constants … … 54 86 ## Functions 55 87 def geodist (plon1, plat1, plon2, plat2) : 56 """Distance between two points (on sphere)""" 57 zs = np.sin (rad*plat1) * np.sin (rad*plat2) + np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) 58 zs = np.maximum (-zone, np.minimum (zone, zs)) 59 geodist = np.arccos (zs) 60 return geodist 88 '''Distance between two points (on sphere)''' 89 zs = ( np.sin (rad*plat1) * np.sin (rad*plat2) + 90 np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) ) 91 zs = np.maximum (-zone, np.minimum (zone, zs)) 92 zgeodist = np.arccos (zs) 93 return zgeodist 61 94 62 95 ### ===== Reading command line parameters ================================================== 63 96 # Creating a parser 64 97 parser = argparse.ArgumentParser ( 65 description = """Compute calving weights""",98 description = '''Compute calving weights''', 66 99 epilog='-------- End of the help message --------') 67 100 68 101 # Adding arguments 69 102 parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', 70 choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 103 choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 104 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 71 105 parser.add_argument ('--atm' , help='atm model name', type=str, default='LMD9695' ) 72 parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', type=int, default=1 ) 73 parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)' , type=int, default=2 ) 74 parser.add_argument ('--atmQuantity' , help='Quantity if atm provides quantities (m/s), Surfacic if atm provided flux (m/s/m2)' , type=str, 75 default='Quantity', choices=['Quantity', 'Surfacic'] ) 76 parser.add_argument ('--oceQuantity' , help='Quantity if oce requires quantities (ks/s), Surfacic if oce requires flux (m/s/m2)', type=str, 77 default='Surfacic', choices=['Quantity', 'Surfacic'] ) 78 parser.add_argument ('--searchRadius' , help='max distance to connect a land point to an ocean point (in km)', type=float, default=600.0 ) 106 parser.add_argument ('--atmCoastWidth', type=int, default=1, 107 help='width of the coastal band in atmosphere (in grid points)' ) 108 parser.add_argument ('--oceCoastWidth', 109 help='width of the coastal band in ocean (in grid points)' , 110 type=int, default=2 ) 111 parser.add_argument ('--atmQuantity' , type=str, default='Quantity', 112 choices=['Quantity', 'Surfacic'], 113 help='Quantity if atm provides quantities (m/s), Surfacic if atm provided flux (m/s/m2)' ) 114 parser.add_argument ('--oceQuantity' , type=str, default='Surfacic', 115 choices=['Quantity', 'Surfacic'], 116 help='Quantity if oce requires quantities (ks/s), Surfacic if oce requires flux (m/s/m2)' ) 117 parser.add_argument ('--searchRadius' , type=float, default=600.0, 118 help='max distance to connect a land point to an ocean point (in km)' ) 79 119 parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) 80 120 parser.add_argument ('--areas' , help='masks file name', default='areas.nc' ) 81 121 parser.add_argument ('--masks' , help='areas file name', default='masks.nc' ) 82 122 parser.add_argument ('--o2a' , help='o2a file name' , default='o2a.nc' ) 83 parser.add_argument ('--output', help='output rmp file name', default='rmp_tlmd_to_torc_runoff.nc' ) 84 parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', 85 choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 86 parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=float, default=0, choices=nemo.nperio_valid_range) 123 parser.add_argument ('--output', help='output rmp file name', 124 default='rmp_tlmd_to_torc_runoff.nc' ) 125 parser.add_argument ('--fmt' , default='netcdf4', 126 help='NetCDF file format, using nco syntax', 127 choices=['classic', 'netcdf3', '64bit', '64bit_data', 128 '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 129 parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=float, 130 default=0, choices=nemo.NPERIO_VALID_RANGE) 131 parser.add_argument ('--modelName' , help='Name of model', type=str, default=None) 87 132 88 133 # Parse command line … … 96 141 97 142 # Model Names 98 atm_Name = myargs.atm 99 oce_Name = myargs.oce 143 atm_Name = myargs.atm 144 oce_Name = myargs.oce 145 model_name = myargs.modelName 100 146 # Width of the coastal band (land points) in the atmopshere 101 147 atmCoastWidth = myargs.atmCoastWidth 102 # Width of the coastal band (ocean points) in the ocean 148 # Width of the coastal band (ocean points) in the ocean 103 149 oceCoastWidth = myargs.oceCoastWidth 104 150 searchRadius = myargs.searchRadius * 1000.0 # From km to meters … … 152 198 atm_grid_area = atm_grid_area.rename ({'ycell':'cell'}) 153 199 atm_grid_imask = atm_grid_imask.rename ({'ycell':'cell'}) 154 200 155 201 if atmDomainType == 'rectilinear' : 156 202 atm_grid_center_lat = atm_grid_center_lat.stack (cell=['y', 'x']) … … 164 210 atm_grid_pmask = atm_grid_imask 165 211 166 167 212 atm_address = np.arange(atm_jpj*atm_jpi) 168 213 … … 186 231 # if oce_jpi == 362 : oce_perio = 6 # ORCA 1 187 232 # if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 188 oce_preio = nemo.__guessNperio__ (oce_jpj, oce_jpi, nperio=oce_perio) 189 190 print ("Oce NPERIO parameter : {:}".format(oce_perio)) 191 oce_grid_imask = nemo.lbc (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T' ).ravel() 192 oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask , (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel() 233 oce_preio = nemo.__guess_nperio__ (oce_jpj, oce_jpi, nperio=oce_perio) 234 235 print ( f'Oce NPERIO parameter : {oce_perio}' ) 236 oce_grid_imask = nemo.lbc (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)), 237 nperio=oce_perio, cd_type='T' ).ravel() 238 oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask , (oce_jpj,oce_jpi)), 239 nperio=oce_perio, cd_type='T', sval=0).ravel() 193 240 oce_address = np.arange(oce_jpj*oce_jpi) 194 241 195 print ( "Fill closed seas with image processing library")242 print ('Fill closed seas with image processing library') 196 243 oce_grid_imask2D = np.reshape (oce_grid_imask, (oce_jpj,oce_jpi)) 197 244 oce_grid_imask2D = 1-nemo.fill_closed_seas (1-oce_grid_imask2D, nperio=oce_perio, cd_type='T') … … 200 247 oce_grid_imask = oce_grid_imask2D.ravel () 201 248 202 ## 203 print ( "Computes an ocean coastal band with image processing library")249 ## 250 print ('Computes an ocean coastal band with image processing library') 204 251 oceLand2D = np.reshape ( np.where (oce_grid_imask < 0.5, True, False), (oce_jpj, oce_jpi) ) 205 252 oceOcean2D = np.reshape ( np.where (oce_grid_imask > 0.5, True, False), (oce_jpj, oce_jpi) ) … … 208 255 oceOceanFiltered2D = ndimage.uniform_filter (oceOcean2D.astype(float), size=NNocean) 209 256 oceCoast2D = np.where (oceOceanFiltered2D<(1.0-0.5/(NNocean**2)),True,False) & oceOcean2D 210 oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0.).ravel() 257 oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D, (oce_jpj,oce_jpi)), 258 nperio=oce_perio, cd_type='T', sval=0.).ravel() 211 259 212 260 oceOceanFiltered = oceOceanFiltered2D.ravel() … … 215 263 oceCoast = oceCoast2D.ravel() 216 264 217 print ( 'Number of points in oceLand : {:8d}'.format (oceLand.sum()))218 print ( 'Number of points in oceOcean : {:8d}'.format (oceOcean.sum()))219 print ( 'Number of points in oceCoast : {:8d}'.format (oceCoast.sum()))265 print ( f'Number of points in oceLand : {oceLand.sum():8d}' ) 266 print ( f'Number of points in oceOcean : {oceOcean.sum():8d}' ) 267 print ( f'Number of points in oceCoast : {oceCoast.sum():8d}' ) 220 268 221 269 # Arrays with coastal points only … … 227 275 oceCoast_address = oce_address [oceCoast] 228 276 229 print ( "Computes an atmosphere coastal band")277 print ('Computes an atmosphere coastal band' ) 230 278 atmLand = np.where (o2aFrac[:] < epsfrac , True, False) 231 279 atmLandFrac = np.where (o2aFrac[:] < zone-epsfrac , True, False) 232 atmFrac = np.where (o2aFrac[:] > epsfrac , True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) 280 atmFrac = np.where (o2aFrac[:] > epsfrac , True, False) & \ 281 np.where (o2aFrac[:] < (zone-epsfrac), True, False) 233 282 atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False) 234 283 atmOceanFrac = np.where (o2aFrac[:] > epsfrac , True, False) … … 236 285 ## For LMDZ only !! 237 286 if atmDomainType == 'rectilinear' : 238 print ( "Extend coastal band ")287 print ('Extend coastal band ' ) 239 288 NNatm = 1+2*atmCoastWidth 240 289 atmLand2D = np.reshape ( atmLand, ( atm_jpj, atm_jpi) ) … … 242 291 atmLandFiltered2D = ndimage.uniform_filter(atmLand2D.astype(float), size=NNatm) 243 292 atmCoast2D = np.where (atmLandFiltered2D<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac 244 293 245 294 atmLandFiltered = atmLandFiltered2D.ravel() 246 295 atmCoast = atmCoast2D.ravel() 247 296 248 print ( 'Number of points in atmLand : {:8d}'.format (atmLand.sum()))249 print ( 'Number of points in atmOcean : {:8d}'.format (atmOcean.sum()))250 print ( 'Number of points in atmCoast : {:8d}'.format (atmCoast.sum()))297 print ( f'Number of points in atmLand : {atmLand.sum():8d}' ) 298 print ( f'Number of points in atmOcean : {atmOcean.sum():8d}' ) 299 print ( f'Number of points in atmCoast : {atmCoast.sum():8d}' ) 251 300 252 301 else : 253 302 atmCoast = atmFrac 254 255 303 304 256 305 # Arrays with coastal points only 257 306 atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast] … … 269 318 ## Loop on atmosphere coastal points 270 319 if searchRadius > 0. : 271 print ( "Loop on atmosphere coastal points")320 print ('Loop on atmosphere coastal points') 272 321 for ja in np.arange (len(atmCoast_grid_pmask)) : 273 z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) 322 z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], 323 oceCoast_grid_center_lon, oceCoast_grid_center_lat) 274 324 z_mask = np.where (z_dist*ra < searchRadius, True, False) 275 325 num_links = int (z_mask.sum()) … … 280 330 if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] 281 331 if ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print 282 print ( 'ja:{:8d}, num_links:{:8d}, z_area:{:8.4e}, atm area:{:8.4e}, weights sum:{:8.4e} '283 .format(ja, num_links, z_area, atm_grid_area[ja].values, poids.sum() ))332 print ( f'ja:{ja:8d} num_links:{num_links:8d} z_area:{z_area:8.4e} '+ 333 f'atm area:{atm_grid_area[ja].values:8.4e} weights sum:{poids.sum():8.4e}' ) 284 334 # 285 335 matrix_local = poids … … 296 346 num_links = remap_matrix.shape[0] 297 347 298 print ( "Write output file")348 print ('Write output file') 299 349 runoff = myargs.output 300 350 print ('Output file: ' + runoff ) 301 351 302 303 remap_matrix = xr.DataArray ( np.reshape(remap_matrix, (num_links, 1)),dims = ['num_links', 'num_wgts'] )352 remap_matrix = xr.DataArray ( np.reshape(remap_matrix, (num_links, 1)), 353 dims = ['num_links', 'num_wgts'] ) 304 354 305 355 # OASIS uses Fortran style indexing, starting at one 306 356 src_address = xr.DataArray ( atm_address.astype(np.int32)+1, dims = ['num_links'], 307 attrs={ "convention": "Fortran style addressing, starting at 1"})357 attrs={'convention': 'Fortran style addressing, starting at 1'}) 308 358 dst_address = xr.DataArray ( oce_address.astype(np.int32)+1, dims = ['num_links'], 309 attrs={"convention": "Fortran style addressing, starting at 1"}) 310 311 src_grid_dims = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32), dims = ['src_grid_rank',] ) 359 attrs={'convention': 'Fortran style addressing, starting at 1'}) 360 361 src_grid_dims = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32), 362 dims = ['src_grid_rank',] ) 312 363 src_grid_center_lon = xr.DataArray (atm_grid_center_lon.values , dims = ['src_grid_size',] ) 313 364 src_grid_center_lat = xr.DataArray (atm_grid_center_lat.values , dims = ['src_grid_size',] ) 314 365 315 src_grid_center_lon.attrs['units']='degrees_east' ; src_grid_center_lon.attrs['long_name']='Longitude' 316 src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" 317 src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' 318 src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" 319 320 src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) 321 src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) 322 src_grid_corner_lon.attrs['units']="degrees_east" 323 src_grid_corner_lat.attrs['units']="degrees_north" 324 325 src_grid_area = xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] ) 326 src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" 327 328 src_grid_imask = xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] ) 329 src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" 330 331 src_grid_pmask = xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] ) 332 src_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; src_grid_pmask.attrs['units']="Land:1, Ocean:0" 366 src_grid_center_lon.attrs.update ( {'units':'degrees_east' , 'long_name':'Longitude', 367 'standard_name':'longitude', 368 'bounds':'src_grid_corner_lon'} ) 369 src_grid_center_lat.attrs.update ( {'units':'degrees_north', 'long_name':'Latitude' , 370 'standard_name':'latitude' , 371 'bounds':'src_grid_corner_lat'} ) 372 373 src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(), 374 dims = [ 'src_grid_size', 'src_grid_corners'] ) 375 src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(), 376 dims = [ 'src_grid_size', 'src_grid_corners'] ) 377 src_grid_corner_lon.attrs['units']='degrees_east' 378 src_grid_corner_lat.attrs['units']='degrees_north' 379 380 src_grid_area = xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] ) 381 src_grid_area.attrs.update ( {'long_name':'Grid area', 'standard_name':'cell_area', 382 'units':'m2' } ) 383 src_grid_imask = xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] ) 384 src_grid_imask.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0'} ) 385 386 src_grid_pmask = xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] ) 387 src_grid_pmask.attrs.update ( {'long_name':'Land-sea mask (periodicity removed)', 388 'units':'Land:1, Ocean:0' } ) 333 389 334 390 # -- 335 dst_grid_dims = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32), dims = ['dst_grid_rank',] ) 391 dst_grid_dims = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32), 392 dims = ['dst_grid_rank',] ) 336 393 dst_grid_center_lon = xr.DataArray (oce_grid_center_lon.values, dims = ['dst_grid_size',] ) 337 394 dst_grid_center_lat = xr.DataArray (oce_grid_center_lat.values, dims = ['dst_grid_size',] ) 338 395 339 dst_grid_center_lon.attrs['units']='degrees_east' ; dst_grid_center_lon.attrs['long_name']='Longitude' 340 dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" 341 dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' 342 dst_grid_center_lat.attrs['long_name']='latitude ' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" 343 344 dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 345 dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 346 dst_grid_corner_lon.attrs['units']="degrees_east" 347 dst_grid_corner_lat.attrs['units']="degrees_north" 348 349 dst_grid_area = xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] ) 350 dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 351 352 dst_grid_imask = xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] ) 353 dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" 354 355 dst_grid_pmask = xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] ) 356 dst_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; dst_grid_pmask.attrs['units']="Land:1, Ocean:0" 357 358 src_lon_addressed = xr.DataArray (atm_grid_center_lon.values[atm_address] , dims = ['num_links'] ) 359 src_lat_addressed = xr.DataArray (atm_grid_center_lat.values[atm_address] , dims = ['num_links'] ) 360 src_area_addressed = xr.DataArray (atm_grid_area .values[atm_address] , dims = ['num_links'] ) 361 src_imask_addressed = xr.DataArray (1-atm_grid_imask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) 362 src_pmask_addressed = xr.DataArray (1-atm_grid_pmask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) 363 364 dst_lon_addressed = xr.DataArray (oce_grid_center_lon.values[atm_address], dims = ['num_links'] ) 365 dst_lat_addressed = xr.DataArray (oce_grid_center_lat.values[oce_address], dims = ['num_links'] ) 366 dst_area_addressed = xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32) , dims = ['num_links'] ) 367 dst_imask_addressed = xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32) , dims = ['num_links'] ) 368 dst_pmask_addressed = xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32) , dims = ['num_links'] ) 369 370 src_lon_addressed.attrs['long_name']="Longitude" ; src_lon_addressed.attrs['standard_name']="longitude" ; src_lon_addressed.attrs['units']="degrees_east" 371 src_lat_addressed.attrs['long_name']="Latitude" ; src_lat_addressed.attrs['standard_name']="latitude" ; src_lat_addressed.attrs['units']="degrees_north" 372 373 dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" 374 dst_lat_addressed.attrs['long_name']="Latitude" ; dst_lat_addressed.attrs['standard_name']="latitude" ; dst_lat_addressed.attrs['units']="degrees_north" 396 dst_grid_center_lon.attrs.update ( {'units':'degrees_east' , 'long_name':'Longitude', 397 'standard_name':'longitude', 398 'bounds':'dst_grid_corner_lon' } ) 399 dst_grid_center_lat.attrs.update ( {'units':'degrees_north', 'long_name':'Latitude' , 400 'standard_name':'latitude' , 401 'bounds':'dst_grid_corner_lat' } ) 402 403 dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values), 404 dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 405 dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values), 406 dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 407 dst_grid_corner_lon.attrs['units']='degrees_east' 408 dst_grid_corner_lat.attrs['units']='degrees_north' 409 410 dst_grid_area = xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] ) 411 dst_grid_area.attrs.update ( {'long_name':'Grid area', 'standard_name':'cell_area', 'units':'m2' }) 412 413 dst_grid_imask = xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] ) 414 dst_grid_imask.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0'} ) 415 416 dst_grid_pmask = xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] ) 417 dst_grid_pmask.attrs.update ( {'long_name':'Land-sea mask (periodicity removed)', 418 'units':'Land:1, Ocean:0'} ) 419 420 src_lon_addressed = xr.DataArray (atm_grid_center_lon.values[atm_address] , 421 dims = ['num_links'] ) 422 src_lat_addressed = xr.DataArray (atm_grid_center_lat.values[atm_address] , 423 dims = ['num_links'] ) 424 src_area_addressed = xr.DataArray (atm_grid_area .values[atm_address] , 425 dims = ['num_links'] ) 426 src_imask_addressed = xr.DataArray (1-atm_grid_imask .values[atm_address].astype(np.int32) , 427 dims = ['num_links'] ) 428 src_pmask_addressed = xr.DataArray (1-atm_grid_pmask .values[atm_address].astype(np.int32) , 429 dims = ['num_links'] ) 430 431 dst_lon_addressed = xr.DataArray (oce_grid_center_lon.values[atm_address], 432 dims = ['num_links'] ) 433 dst_lat_addressed = xr.DataArray (oce_grid_center_lat.values[oce_address], 434 dims = ['num_links'] ) 435 dst_area_addressed = xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32) , 436 dims = ['num_links'] ) 437 dst_imask_addressed = xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32) , 438 dims = ['num_links'] ) 439 dst_pmask_addressed = xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32) , 440 dims = ['num_links'] ) 441 442 src_lon_addressed.attrs.update ({'units':'degrees_east' , 'long_name':'Longitude', 443 'standard_name':'longitude' }) 444 src_lat_addressed.attrs.update ({'units':'degrees_north', 'long_name':'Latitude' , 445 'standard_name':'latitude' }) 446 447 dst_lon_addressed.attrs.update ({'units':'degrees_east' , 'long_name':'Longitude', 448 'standard_name':'longitude' }) 449 dst_lat_addressed.attrs.update ({'units':'degrees_north', 'long_name':'Latitude' , 450 'standard_name':'latitude' }) 375 451 376 452 if atmDomainType == 'rectilinear' : … … 390 466 391 467 f_runoff = xr.Dataset ( { 392 'remap_matrix' 468 'remap_matrix' : remap_matrix, 393 469 'src_address' : src_address, 394 470 'dst_address' : dst_address, … … 398 474 'src_grid_corner_lon' : src_grid_corner_lon, 399 475 'src_grid_corner_lat' : src_grid_corner_lat, 400 'src_grid_area' : src_grid_area,401 476 'src_grid_area' : src_grid_area, 402 477 'src_grid_imask' : src_grid_imask, … … 425 500 } ) 426 501 427 f_runoff.attrs['Conventions'] = "CF-1.6" 428 f_runoff.attrs['source'] = "IPSL Earth system model" 429 f_runoff.attrs['group'] = "ICMC IPSL Climate Modelling Center" 430 f_runoff.attrs['Institution'] = "IPSL https.//www.ipsl.fr" 431 f_runoff.attrs['Ocean'] = oce_Name + " https://www.nemo-ocean.eu" 432 f_runoff.attrs['Atmosphere'] = atm_Name + " http://lmdz.lmd.jussieu.fr" 433 f_runoff.attrs['associatedFiles'] = grids + " " + areas + " " + masks 434 f_runoff.attrs['description'] = "Generated with RunoffWeights.py" 435 f_runoff.attrs['title'] = runoff 436 f_runoff.attrs['Program'] = "Generated by " + sys.argv[0] + " with flags " + ' '.join (sys.argv[1:]) 437 f_runoff.attrs['atmCoastWidth'] = "{:d} grid points".format(atmCoastWidth) 438 f_runoff.attrs['oceCoastWidth'] = "{:d} grid points".format(oceCoastWidth) 439 f_runoff.attrs['searchRadius'] = "{:.0f} km".format(searchRadius/1000.) 440 f_runoff.attrs['atmQuantity'] = myargs.atmQuantity 441 f_runoff.attrs['oceQuantity'] = myargs.oceQuantity 442 f_runoff.attrs['gridsFile'] = grids 443 f_runoff.attrs['areasFile'] = areas 444 f_runoff.attrs['masksFile'] = masks 445 f_runoff.attrs['o2aFile'] = o2a 446 f_runoff.attrs['timeStamp'] = time.asctime () 447 try : f_calving.attrs['directory'] = os.getcwd () 448 except : pass 449 try : f_runoff.attrs['HOSTNAME'] = platform.node () 450 except : pass 451 try : f_runoff.attrs['LOGNAME'] = os.getlogin () 452 except : pass 453 try : f_runoff.attrs['Python'] = "Python version " + platform.python_version () 454 except : pass 455 try : f_runoff.attrs['OS'] = platform.system () 456 except : pass 457 try : f_runoff.attrs['release'] = platform.release () 458 except : pass 459 try : f_runoff.attrs['hardware'] = platform.machine () 460 except : pass 461 f_runoff.attrs['conventions'] = "SCRIP" 462 f_runoff.attrs['source_grid'] = "curvilinear" 463 f_runoff.attrs['dest_grid'] = "curvilinear" 464 f_runoff.attrs['Model'] = "IPSL CM6" 465 f_runoff.attrs['SVN_Author'] = "$Author$" 466 f_runoff.attrs['SVN_Date'] = "$Date$" 467 f_runoff.attrs['SVN_Revision'] = "$Revision$" 468 f_runoff.attrs['SVN_Id'] = "$Id$" 469 f_runoff.attrs['SVN_HeadURL'] = "$HeadURL$" 502 f_runoff.attrs.update ( { 503 'Conventions' : 'CF-1.6', 504 'source' : 'IPSL Earth system model', 505 'group' : 'ICMC IPSL Climate Modelling Center', 506 'Institution' : 'IPSL https.//www.ipsl.fr', 507 'Ocean' : f'{oce_Name} https://www.nemo-ocean.eu', 508 'Atmosphere' : f'{atm_Name} http://lmdz.lmd.jussieu.fr', 509 'associatedFiles' : f'{grids} {areas} {masks}', 510 'description' : 'Generated with RunoffWeights.py', 511 'title' : runoff, 512 'Program' : f'Generated by {sys.argv[0]} with flags ' + ' '.join (sys.argv[1:]), 513 'atmCoastWidth' : f'{atmCoastWidth:d} grid points', 514 'oceCoastWidth' : f'{oceCoastWidth:d} grid points', 515 'searchRadius' : f'{searchRadius/1000.:.0f} km', 516 'atmQuantity' : myargs.atmQuantity, 517 'oceQuantity' : myargs.oceQuantity, 518 'gridsFile' : grids, 519 'areasFile' : areas, 520 'masksFile' : masks, 521 'o2aFile' : o2a, 522 'timeStamp' : time.asctime (), 523 'directory' : os.getcwd (), 524 'HOSTNAME' : platform.node (), 525 'LOGNAME' : getpass.getuser(), 526 'Python' : f'Python version {platform.python_version ()}', 527 'OS' : platform.system (), 528 'release' : platform.release (), 529 'hardware' : platform.machine (), 530 'conventions' : 'SCRIP', 531 'source_grid' : 'curvilinear', 532 'dest_grid' : 'curvilinear', 533 'Model' : model_name, 534 'SVN_Author' : '$Author$', 535 'SVN_Date' : '$Date$', 536 'SVN_Revision' : '$Revision$', 537 'SVN_Id' : '$Id$', 538 'SVN_HeadURL' : '$HeadURL$', 539 } ) 470 540 471 541 f_runoff.to_netcdf ( runoff, mode='w', format=FmtNetcdf ) … … 474 544 ## 475 545 476 print ( 'That''s all folks !')546 print ("That''s all folks !") 477 547 ## ====================================================================================== -
TOOLS/MOSAIX/nemo.py
r6245 r6666 18 18 ## 19 19 ## =========================================================================== 20 '''Utilities to plot NEMO ORCA fields, 21 22 Handles periodicity and other stuff 23 24 - Lots of tests for xarray object 25 - Not much tested for numpy objects 26 27 Author: olivier.marti@lsce.ipsl.fr 28 29 ## SVN information 30 Author = "$Author$" 31 Date = "$Date$" 32 Revision = "$Revision$" 33 Id = "$Id$" 34 HeadURL = "$HeadURL$" 20 35 ''' 21 Utilities to plot NEMO ORCA fields22 Periodicity and other stuff23 24 olivier.marti@lsce.ipsl.fr25 '''26 27 ## SVN information28 __Author__ = "$Author$"29 __Date__ = "$Date$"30 __Revision__ = "$Revision$"31 __Id__ = "$Id$"32 __HeadURL = "$HeadURL$"33 36 34 37 import numpy as np 35 try : import xarray as xr 36 except ImportError : pass 37 38 try : import f90nml 39 except : pass 40 41 try : from sklearn.impute import SimpleImputer 42 except : pass 43 44 try : import numba 45 except : pass 46 47 rpi = np.pi ; rad = np.deg2rad (1.0) ; dar = np.rad2deg (1.0) 48 49 nperio_valid_range = [0, 1, 4, 4.2, 5, 6, 6.2] 50 51 rday = 24.*60.*60. # Day length [s] 52 rsiyea = 365.25 * rday * 2. * rpi / 6.283076 # Sideral year length [s] 53 rsiday = rday / (1. + rday / rsiyea) 54 raamo = 12. # Number of months in one year 55 rjjhh = 24. # Number of hours in one day 56 rhhmm = 60. # Number of minutes in one hour 57 rmmss = 60. # Number of seconds in one minute 58 omega = 2. * rpi / rsiday # Earth rotation parameter [s-1] 59 ra = 6371229. # Earth radius [m] 60 grav = 9.80665 # Gravity [m/s2] 61 repsi = np.finfo (1.0).eps 62 63 xList = [ 'x', 'X', 'lon' , 'longitude' ] 64 yList = [ 'y', 'Y', 'lat' , 'latitude' ] 65 zList = [ 'z', 'Z', 'depth' , ] 66 tList = [ 't', 'T', 'time' , ] 38 import xarray as xr 39 40 # Tries to import some moldules that are available in all Python installations 41 # and used in only a few specific functions 42 try : 43 from sklearn.impute import SimpleImputer 44 except ImportError as err : 45 print ( f'===> Warning : Module nemo : Import error of sklearn.impute.SimpleImputer : {err}' ) 46 SimpleImputer = None 47 48 try : 49 import f90nml 50 except ImportError as err : 51 print ( f'===> Warning : Module nemo : Import error of f90nml : {err}' ) 52 f90nml = None 53 54 # SVN information 55 __SVN__ = ({ 56 'Author' : '$Author$', 57 'Date' : '$Date$', 58 'Revision' : '$Revision$', 59 'Id' : '$Id$', 60 'HeadURL' : '$HeadURL$', 61 'SVN_Date' : '$SVN_Date: $', 62 }) 63 ## 64 65 RPI = np.pi 66 RAD = np.deg2rad (1.0) 67 DAR = np.rad2deg (1.0) 68 REPSI = np.finfo (1.0).eps 69 70 NPERIO_VALID_RANGE = [0, 1, 4, 4.2, 5, 6, 6.2] 71 72 RAAMO = 12 # Number of months in one year 73 RJJHH = 24 # Number of hours in one day 74 RHHMM = 60 # Number of minutes in one hour 75 RMMSS = 60 # Number of seconds in one minute 76 RA = 6371229.0 # Earth radius [m] 77 GRAV = 9.80665 # Gravity [m/s2] 78 RT0 = 273.15 # Freezing point of fresh water [Kelvin] 79 RAU0 = 1026.0 # Volumic mass of sea water [kg/m3] 80 SICE = 6.0 # Salinity of ice (for pisces) [psu] 81 SOCE = 34.7 # Salinity of sea (for pisces and isf) [psu] 82 RLEVAP = 2.5e+6 # Latent heat of evaporation (water) [J/K] 83 VKARMN = 0.4 # Von Karman constant 84 STEFAN = 5.67e-8 # Stefan-Boltzmann constant [W/m2/K4] 85 RHOS = 330. # Volumic mass of snow [kg/m3] 86 RHOI = 917. # Volumic mass of sea ice [kg/m3] 87 RHOW = 1000. # Volumic mass of freshwater in melt ponds [kg/m3] 88 RCND_I = 2.034396 # Thermal conductivity of fresh ice [W/m/K] 89 RCPI = 2067.0 # Specific heat of fresh ice [J/kg/K] 90 RLSUB = 2.834e+6 # Pure ice latent heat of sublimation [J/kg] 91 RLFUS = 0.334e+6 # Latent heat of fusion of fresh ice [J/kg] 92 RTMLT = 0.054 # Decrease of seawater meltpoint with salinity 93 94 RDAY = RJJHH * RHHMM * RMMSS # Day length [s] 95 RSIYEA = 365.25 * RDAY * 2. * RPI / 6.283076 # Sideral year length [s] 96 RSIDAY = RDAY / (1. + RDAY / RSIYEA) # Sideral day length [s] 97 OMEGA = 2. * RPI / RSIDAY # Earth rotation parameter [s-1] 98 99 ## Default names of dimensions 100 UDIMS = {'x':'x', 'y':'y', 'z':'olevel', 't':'time_counter'} 101 102 ## All possibles name of dimensions in Nemo files 103 XNAME = [ 'x', 'X', 'X1', 'xx', 'XX', 104 'x_grid_T', 'x_grid_U', 'x_grid_V', 'x_grid_F', 'x_grid_W', 105 'lon', 'nav_lon', 'longitude', 'X1', 'x_c', 'x_f', ] 106 YNAME = [ 'y', 'Y', 'Y1', 'yy', 'YY', 107 'y_grid_T', 'y_grid_U', 'y_grid_V', 'y_grid_F', 'y_grid_W', 108 'lat', 'nav_lat', 'latitude' , 'Y1', 'y_c', 'y_f', ] 109 ZNAME = [ 'z', 'Z', 'Z1', 'zz', 'ZZ', 'depth', 'tdepth', 'udepth', 110 'vdepth', 'wdepth', 'fdepth', 'deptht', 'depthu', 111 'depthv', 'depthw', 'depthf', 'olevel', 'z_c', 'z_f', ] 112 TNAME = [ 't', 'T', 'tt', 'TT', 'time', 'time_counter', 'time_centered', ] 113 114 ## All possibles name of units of dimensions in Nemo files 115 XUNIT = [ 'degrees_east', ] 116 YUNIT = [ 'degrees_north', ] 117 ZUNIT = [ 'm', 'meter', ] 118 TUNIT = [ 'second', 'minute', 'hour', 'day', 'month', 'year', ] 119 120 ## All possibles size of dimensions in Orca files 121 XLENGTH = [ 180, 182, 360, 362, 1440 ] 122 YLENGTH = [ 148, 149, 331, 332 ] 123 ZLENGTH = [ 31, 75] 67 124 68 125 ## =========================================================================== 69 def __mmath__ (tab, default=None) : 126 def __mmath__ (ptab, default=None) : 127 '''Determines the type of tab : xarray, numpy or numpy.ma object ? 128 129 Returns type 130 ''' 70 131 mmath = default 71 try : 72 if type (tab) == xr.core.dataarray.DataArray : mmath = xr 73 except : 74 pass 75 76 try : 77 if type (tab) == np.ndarray : mmath = np 78 except : 79 pass 80 132 if isinstance (ptab, xr.core.dataarray.DataArray) : 133 mmath = xr 134 if isinstance (ptab, np.ndarray) : 135 mmath = np 136 if isinstance (ptab, np.ma.MaskType) : 137 mmath = np.ma 138 81 139 return mmath 82 140 83 84 def __guessNperio__ (jpj, jpi, nperio=None, out='nperio') : 85 ''' 86 Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details) 87 141 def __guess_nperio__ (jpj, jpi, nperio=None, out='nperio') : 142 '''Tries to guess the value of nperio (periodicity parameter. 143 144 See NEMO documentation for details) 88 145 Inputs 89 146 jpj : number of latitudes … … 91 148 nperio : periodicity parameter 92 149 ''' 93 if nperio == None : 94 nperio = __guessConfig__ (jpj, jpi, nperio=None, out='nperio') 95 150 if nperio is None : 151 nperio = __guess_config__ (jpj, jpi, nperio=None, out=out) 96 152 return nperio 97 153 98 def __guess Config__ (jpj, jpi, nperio=None, config=None, out='nperio') :99 ''' 100 Tries to guess the value of nperio (periodicity parameter. See NEMO documentation for details) 101 154 def __guess_config__ (jpj, jpi, nperio=None, config=None, out='nperio') : 155 '''Tries to guess the value of nperio (periodicity parameter). 156 157 See NEMO documentation for details) 102 158 Inputs 103 159 jpj : number of latitudes … … 105 161 nperio : periodicity parameter 106 162 ''' 107 if nperio == None : 163 print ( jpi, jpj) 164 if nperio is None : 108 165 ## Values for NEMO version < 4.2 109 if jpj == 149 and jpi == 182 : 110 config = 'ORCA2.3' 111 nperio = 4 # ORCA2. We choose legacy orca2. 112 Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'T' 113 if jpj == 332 and jpi == 362 : # eORCA1. 114 config = 'eORCA1.2' 115 nperio = 6 116 Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 166 if ( (jpj == 149 and jpi == 182) or (jpj is None and jpi == 182) or 167 (jpj == 149 or jpi is None) ) : 168 # ORCA2. We choose legacy orca2. 169 config, nperio, iperio, jperio, nfold, nftype = 'ORCA2.3' , 4, 1, 0, 1, 'T' 170 if ((jpj == 332 and jpi == 362) or (jpj is None and jpi == 362) or 171 (jpj == 332 and jpi is None) ) : # eORCA1. 172 config, nperio, iperio, jperio, nfold, nftype = 'eORCA1.2', 6, 1, 0, 1, 'F' 117 173 if jpi == 1442 : # ORCA025. 118 config = 'ORCA025' 119 nperio = 6 120 Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 174 config, nperio, iperio, jperio, nfold, nftype = 'ORCA025' , 6, 1, 0, 1, 'F' 121 175 if jpj == 294 : # ORCA1 122 config = 'ORCA1' 123 nperio = 6 124 Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 125 176 config, nperio, iperio, jperio, nfold, nftype = 'ORCA1' , 6, 1, 0, 1, 'F' 177 126 178 ## Values for NEMO version >= 4.2. No more halo points 127 if jpj == 148 and jpi == 180 : 128 config = 'ORCA2.4' 129 nperio = 4.2 # ORCA2. We choose legacy orca2. 130 Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 131 if jpj == 331 and jpi == 360 : # eORCA1. 132 config = 'eORCA1.4' 133 nperio = 6.2 134 Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 179 if (jpj == 148 and jpi == 180) or (jpj is None and jpi == 180) or \ 180 (jpj == 148 and jpi is None) : # ORCA2. We choose legacy orca2. 181 config, nperio, iperio, jperio, nfold, nftype = 'ORCA2.4' , 4.2, 1, 0, 1, 'F' 182 if (jpj == 331 and jpi == 360) or (jpj is None and jpi == 360) or \ 183 (jpj == 331 and jpi is None) : # eORCA1. 184 config, nperio, iperio, jperio, nfold, nftype = 'eORCA1.4', 6.2, 1, 0, 1, 'F' 135 185 if jpi == 1440 : # ORCA025. 136 config = 'ORCA025' 137 nperio = 6.2 138 Iperio = 1 ; Jperio = 0 ; NFold = 1 ; NFtype = 'F' 139 140 if nperio == None : 141 raise Exception ('in nemo module : nperio not found, and cannot by guessed') 186 config, nperio, iperio, jperio, nfold, nftype = 'ORCA025' , 6.2, 1, 0, 1, 'F' 187 188 if nperio is None : 189 raise ValueError ('in nemo module : nperio not found, and cannot by guessed') 190 191 if nperio in NPERIO_VALID_RANGE : 192 print ( f'nperio set as {nperio} (deduced from {jpj=} and {jpi=})' ) 142 193 else : 143 if nperio in nperio_valid_range : 144 print ('nperio set as {:} (deduced from jpj={:d} jpi={:d})'.format (nperio, jpj, jpi)) 145 else : 146 raise ValueError ('nperio set as {:} (deduced from jpi={:d}) : nemo.py is not ready for this value'.format (nperio, jpi)) 147 148 if out == 'nperio' : return nperio 149 if out == 'config' : return config 150 if out == 'perio' : return Iperio, Jperio, NFold, NFtype 151 if out in ['full', 'all'] : return {'nperio':nperio, 'Iperio':Iperio, 'Jperio':Jperio, 'NFold':NFold, 'NFtype':NFtype} 152 153 def __guessPoint__ (ptab) : 154 ''' 155 Tries to guess the grid point (periodicity parameter. See NEMO documentation for details) 156 194 raise ValueError ( f'nperio set as {nperio} (deduced from {jpi=} and {jpj=}) : \n'+ 195 'nemo.py is not ready for this value' ) 196 197 if out == 'nperio' : 198 return nperio 199 if out == 'config' : 200 return config 201 if out == 'perio' : 202 return iperio, jperio, nfold, nftype 203 if out in ['full', 'all'] : 204 return {'nperio':nperio, 'iperio':iperio, 'jperio':jperio, 'nfold':nfold, 'nftype':nftype} 205 206 def __guess_point__ (ptab) : 207 '''Tries to guess the grid point (periodicity parameter. 208 209 See NEMO documentation for details) 157 210 For array conforments with xgcm requirements 158 211 … … 162 215 Credits : who is the original author ? 163 216 ''' 164 gP = None 217 218 gp = None 165 219 mmath = __mmath__ (ptab) 166 220 if mmath == xr : 167 if 'x_c' in ptab.dims and 'y_c' in ptab.dims : gP = 'T' 168 if 'x_f' in ptab.dims and 'y_c' in ptab.dims : gP = 'U' 169 if 'x_c' in ptab.dims and 'y_f' in ptab.dims : gP = 'V' 170 if 'x_f' in ptab.dims and 'y_f' in ptab.dims : gP = 'F' 171 if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_c' in ptab.dims : gP = 'T' 172 if 'x_c' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'W' 173 if 'x_f' in ptab.dims and 'y_c' in ptab.dims and 'z_f' in ptab.dims : gP = 'U' 174 if 'x_c' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'V' 175 if 'x_f' in ptab.dims and 'y_f' in ptab.dims and 'z_f' in ptab.dims : gP = 'F' 176 177 if gP == None : 178 raise Exception ('in nemo module : cd_type not found, and cannot by guessed') 179 else : 180 print ('Grid set as', gP, 'deduced from dims ', ptab.dims) 181 return gP 182 else : 183 raise Exception ('in nemo module : cd_type not found, input is not an xarray data') 221 if ('x_c' in ptab.dims and 'y_c' in ptab.dims ) : 222 gp = 'T' 223 if ('x_f' in ptab.dims and 'y_c' in ptab.dims ) : 224 gp = 'U' 225 if ('x_c' in ptab.dims and 'y_f' in ptab.dims ) : 226 gp = 'V' 227 if ('x_f' in ptab.dims and 'y_f' in ptab.dims ) : 228 gp = 'F' 229 if ('x_c' in ptab.dims and 'y_c' in ptab.dims 230 and 'z_c' in ptab.dims ) : 231 gp = 'T' 232 if ('x_c' in ptab.dims and 'y_c' in ptab.dims 233 and 'z_f' in ptab.dims ) : 234 gp = 'W' 235 if ('x_f' in ptab.dims and 'y_c' in ptab.dims 236 and 'z_f' in ptab.dims ) : 237 gp = 'U' 238 if ('x_c' in ptab.dims and 'y_f' in ptab.dims 239 and 'z_f' in ptab.dims ) : 240 gp = 'V' 241 if ('x_f' in ptab.dims and 'y_f' in ptab.dims 242 and 'z_f' in ptab.dims ) : 243 gp = 'F' 244 245 if gp is None : 246 raise AttributeError ('in nemo module : cd_type not found, and cannot by guessed') 247 print ( f'Grid set as {gp} deduced from dims {ptab.dims}' ) 248 return gp 249 else : 250 raise AttributeError ('in nemo module : cd_type not found, input is not an xarray data') 251 252 def get_shape ( ptab ) : 253 '''Get shape of ptab return a string with axes names 254 255 shape may contain X, Y, Z or T 256 Y is missing for a latitudinal slice 257 X is missing for on longitudinal slice 258 etc ... 259 ''' 260 261 g_shape = '' 262 if __find_axis__ (ptab, 'x')[0] : 263 g_shape = 'X' 264 if __find_axis__ (ptab, 'y')[0] : 265 g_shape = 'Y' + g_shape 266 if __find_axis__ (ptab, 'z')[0] : 267 g_shape = 'Z' + g_shape 268 if __find_axis__ (ptab, 't')[0] : 269 g_shape = 'T' + g_shape 270 return g_shape 184 271 185 272 def lbc_diag (nperio) : 186 lperio = nperio ; aperio = False 273 '''Useful to switch between field with and without halo''' 274 lperio, aperio = nperio, False 187 275 if nperio == 4.2 : 188 lperio = 4 ; aperio =True276 lperio, aperio = 4, True 189 277 if nperio == 6.2 : 190 lperio = 6 ; aperio = True 191 192 return lperio, aperio 193 194 def __findAxis__ (tab, axis='z') : 195 ''' 196 Find number and name of the requested axis 197 ''' 198 mmath = __mmath__ (tab) 199 ix = None ; ax = None 200 201 if axis in xList : 202 axList = [ 'x', 'X', 203 'lon', 'nav_lon', 'nav_lon_T', 'nav_lon_U', 'nav_lon_V', 'nav_lon_F', 'nav_lon_W', 204 'x_grid_T', 'x_grid_U', 'x_grid_V', 'x_grid_F', 'x_grid_W', 205 'glam', 'glamt', 'glamu', 'glamv', 'glamf', 'glamw' ] 206 unList = [ 'degrees_east' ] 207 if axis in yList : 208 axList = [ 'y', 'Y', 'lat', 209 'nav_lat', 'nav_lat_T', 'nav_lat_U', 'nav_lat_V', 'nav_lat_F', 'nav_lat_W', 210 'y_grid_T', 'y_grid_U', 'y_grid_V', 'y_grid_F', 'y_grid_W', 211 'gphi', 'gphi', 'gphiu', 'gphiv', 'gphif', 'gphiw'] 212 unList = [ 'degrees_north' ] 213 if axis in zList : 214 axList = [ 'z', 'Z', 215 'depth', 'deptht', 'depthu', 'depthv', 'depthf', 'depthw', 216 'olevel' ] 217 unList = [ 'm', 'meter' ] 218 if axis in tList : 219 axList = [ 't', 'T', 'time', 'time_counter' ] 220 unList = [ 'second', 'minute', 'hour', 'day', 'month' ] 221 278 lperio, aperio = 6, True 279 return lperio, aperio 280 281 def __find_axis__ (ptab, axis='z', back=True) : 282 '''Returns name and name of the requested axis''' 283 mmath = __mmath__ (ptab) 284 ax, ix = None, None 285 286 if axis in XNAME : 287 ax_name, unit_list, length = XNAME, XUNIT, XLENGTH 288 if axis in YNAME : 289 ax_name, unit_list, length = YNAME, YUNIT, YLENGTH 290 if axis in ZNAME : 291 ax_name, unit_list, length = ZNAME, ZUNIT, ZLENGTH 292 if axis in TNAME : 293 ax_name, unit_list, length = TNAME, TUNIT, None 294 222 295 if mmath == xr : 223 for Name in axList : 224 try : 225 ix = tab.dims.index (Name) 226 ax = Name 227 except : pass 228 229 for i, dim in enumerate (tab.dims) : 230 if 'units' in tab.coords[dim].attrs.keys() : 231 for name in unList : 232 if name in tab.coords[dim].attrs['units'] : 233 ix = i 234 ax = dim 235 else : 236 if axis in xList : ix=-1 237 if axis in yList : 238 if len(tab.shape) >= 2 : ix=-2 239 if axis in zList : 240 if len(tab.shape) >= 3 : ix=-3 241 if axis in tList : 242 if len(tab.shape) >=3 : ix=-3 243 if len(tab.shape) >=4 : ix=-4 244 245 return ix, ax 246 247 #@numba.jit(forceobj=True) 248 def fixed_lon (lon, center_lon=0.0) : 249 ''' 250 Returns corrected longitudes for nicer plots 296 # Try by name 297 for dim in ax_name : 298 if dim in ptab.dims : 299 ix, ax = ptab.dims.index (dim), dim 300 301 # If not found, try by axis attributes 302 if not ix : 303 for i, dim in enumerate (ptab.dims) : 304 if 'axis' in ptab.coords[dim].attrs.keys() : 305 l_axis = ptab.coords[dim].attrs['axis'] 306 if axis in ax_name and l_axis == 'X' : 307 ix, ax = (i, dim) 308 if axis in ax_name and l_axis == 'Y' : 309 ix, ax = (i, dim) 310 if axis in ax_name and l_axis == 'Z' : 311 ix, ax = (i, dim) 312 if axis in ax_name and l_axis == 'T' : 313 ix, ax = (i, dim) 314 315 # If not found, try by units 316 if not ix : 317 for i, dim in enumerate (ptab.dims) : 318 if 'units' in ptab.coords[dim].attrs.keys() : 319 for name in unit_list : 320 if name in ptab.coords[dim].attrs['units'] : 321 ix, ax = i, dim 322 323 # If numpy array or dimension not found, try by length 324 if mmath != xr or not ix : 325 if length : 326 l_shape = ptab.shape 327 for nn in np.arange ( len(l_shape) ) : 328 if l_shape[nn] in length : 329 ix = nn 330 331 if ix and back : 332 ix -= len(ptab.shape) 333 334 return ax, ix 335 336 def find_axis ( ptab, axis='z', back=True ) : 337 '''Version of find_axis with no __''' 338 ix, xx = __find_axis__ (ptab, axis, back) 339 return xx, ix 340 341 def fixed_lon (plon, center_lon=0.0) : 342 '''Returns corrected longitudes for nicer plots 251 343 252 344 lon : longitudes of the grid. At least 2D. 253 345 center_lon : center longitude. Default=0. 254 346 255 Designed by Phil Pelson. See https://gist.github.com/pelson/79cf31ef324774c97ae7 256 ''' 257 mmath = __mmath__ (lon) 347 Designed by Phil Pelson. 348 See https://gist.github.com/pelson/79cf31ef324774c97ae7 349 ''' 350 mmath = __mmath__ (plon) 351 352 f_lon = plon.copy () 353 354 f_lon = mmath.where (f_lon > center_lon+180., f_lon-360.0, f_lon) 355 f_lon = mmath.where (f_lon < center_lon-180., f_lon+360.0, f_lon) 356 357 for i, start in enumerate (np.argmax (np.abs (np.diff (f_lon, axis=-1)) > 180., axis=-1)) : 358 f_lon [..., i, start+1:] += 360. 359 360 # Special case for eORCA025 361 if f_lon.shape [-1] == 1442 : 362 f_lon [..., -2, :] = f_lon [..., -3, :] 363 if f_lon.shape [-1] == 1440 : 364 f_lon [..., -1, :] = f_lon [..., -2, :] 365 366 if f_lon.min () > center_lon : 367 f_lon += -360.0 368 if f_lon.max () < center_lon : 369 f_lon += 360.0 370 371 if f_lon.min () < center_lon-360.0 : 372 f_lon += 360.0 373 if f_lon.max () > center_lon+360.0 : 374 f_lon += -360.0 375 376 return f_lon 377 378 def bounds_clolon ( pbounds_lon, plon, rad=False, deg=True) : 379 '''Choose closest to lon0 longitude, adding/substacting 360° if needed 380 ''' 381 382 if rad : 383 lon_range = 2.0*np.pi 384 if deg : 385 lon_range = 360.0 386 b_clolon = pbounds_lon.copy () 387 388 b_clolon = xr.where ( b_clolon < plon-lon_range/2., 389 b_clolon+lon_range, 390 b_clolon ) 391 b_clolon = xr.where ( b_clolon > plon+lon_range/2., 392 b_clolon-lon_range, 393 b_clolon ) 394 return b_clolon 395 396 def unify_dims ( dd, x='x', y='y', z='olevel', t='time_counter', verbose=False ) : 397 '''Rename dimensions to unify them between NEMO versions 398 ''' 399 for xx in XNAME : 400 if xx in dd.dims and xx != x : 401 if verbose : 402 print ( f"{xx} renamed to {x}" ) 403 dd = dd.rename ( {xx:x}) 404 405 for yy in YNAME : 406 if yy in dd.dims and yy != y : 407 if verbose : 408 print ( f"{yy} renamed to {y}" ) 409 dd = dd.rename ( {yy:y} ) 410 411 for zz in ZNAME : 412 if zz in dd.dims and zz != z : 413 if verbose : 414 print ( f"{zz} renamed to {z}" ) 415 dd = dd.rename ( {zz:z} ) 416 417 for tt in TNAME : 418 if tt in dd.dims and tt != t : 419 if verbose : 420 print ( f"{tt} renamed to {t}" ) 421 dd = dd.rename ( {tt:t} ) 422 423 return dd 424 425 426 if SimpleImputer : 427 def fill_empty (ptab, sval=np.nan, transpose=False) : 428 '''Fill empty values 429 430 Useful when NEMO has run with no wet points options : 431 some parts of the domain, with no ocean points, have no 432 values 433 ''' 434 mmath = __mmath__ (ptab) 435 436 imp = SimpleImputer (missing_values=sval, strategy='mean') 437 if transpose : 438 imp.fit (ptab.T) 439 ztab = imp.transform (ptab.T).T 440 else : 441 imp.fit (ptab) 442 ztab = imp.transform (ptab) 443 444 if mmath == xr : 445 ztab = xr.DataArray (ztab, dims=ztab.dims, coords=ztab.coords) 446 ztab.attrs.update (ptab.attrs) 447 448 return ztab 449 258 450 259 fixed_lon = lon.copy () 260 261 fixed_lon = mmath.where (fixed_lon > center_lon+180., fixed_lon-360.0, fixed_lon) 262 fixed_lon = mmath.where (fixed_lon < center_lon-180., fixed_lon+360.0, fixed_lon) 263 264 for i, start in enumerate (np.argmax (np.abs (np.diff (fixed_lon, axis=-1)) > 180., axis=-1)) : 265 fixed_lon [..., i, start+1:] += 360. 266 267 # Special case for eORCA025 268 if fixed_lon.shape [-1] == 1442 : fixed_lon [..., -2, :] = fixed_lon [..., -3, :] 269 if fixed_lon.shape [-1] == 1440 : fixed_lon [..., -1, :] = fixed_lon [..., -2, :] 270 271 if fixed_lon.min () > center_lon : fixed_lon += -360.0 272 if fixed_lon.max () < center_lon : fixed_lon += 360.0 273 274 if fixed_lon.min () < center_lon-360.0 : fixed_lon += 360.0 275 if fixed_lon.max () > center_lon+360.0 : fixed_lon += -360.0 276 277 return fixed_lon 278 279 #@numba.jit(forceobj=True) 280 def fill_empty (ztab, sval=np.nan, transpose=False) : 281 ''' 282 Fill values 283 284 Useful when NEMO has run with no wet points options : 285 some parts of the domain, with no ocean points, has no 451 else : 452 print ("Import error of sklearn.impute.SimpleImputer") 453 def fill_empty (ptab, sval=np.nan, transpose=False) : 454 '''Void version of fill_empy, because module sklearn.impute.SimpleImputer is not available 455 456 fill_empty : 457 Fill values 458 459 Useful when NEMO has run with no wet points options : 460 some parts of the domain, with no ocean points, have no 461 values 462 ''' 463 print ( 'Error : module sklearn.impute.SimpleImputer not found' ) 464 print ( 'Can not call fill_empty' ) 465 print ( 'Call arguments where : ' ) 466 print ( f'{ptab.shape=} {sval=} {transpose=}' ) 467 468 def fill_lonlat (plon, plat, sval=-1) : 469 '''Fill longitude/latitude values 470 471 Useful when NEMO has run with no wet points options : 472 some parts of the domain, with no ocean points, have no 286 473 lon/lat values 287 474 ''' 288 mmath = __mmath__ (ztab) 475 from sklearn.impute import SimpleImputer 476 mmath = __mmath__ (plon) 289 477 290 478 imp = SimpleImputer (missing_values=sval, strategy='mean') 291 if transpose : 292 imp.fit (ztab.T) 293 ptab = imp.transform (ztab.T).T 294 else : 295 imp.fit (ztab) 296 ptab = imp.transform (ztab) 297 479 imp.fit (plon) 480 zlon = imp.transform (plon) 481 imp.fit (plat.T) 482 zlat = imp.transform (plat.T).T 483 298 484 if mmath == xr : 299 ptab = xr.DataArray (ptab, dims=ztab.dims, coords=ztab.coords) 300 ptab.attrs = ztab.attrs 301 302 return ptab 303 304 #@numba.jit(forceobj=True) 305 def fill_lonlat (lon, lat, sval=-1) : 306 ''' 307 Fill longitude/latitude values 308 309 Useful when NEMO has run with no wet points options : 485 zlon = xr.DataArray (zlon, dims=plon.dims, coords=plon.coords) 486 zlat = xr.DataArray (zlat, dims=plat.dims, coords=plat.coords) 487 zlon.attrs.update (plon.attrs) 488 zlat.attrs.update (plat.attrs) 489 490 zlon = fixed_lon (zlon) 491 492 return zlon, zlat 493 494 def fill_bounds_lonlat (pbounds_lon, pbounds_lat, sval=-1) : 495 '''Fill longitude/latitude bounds values 496 497 Useful when NEMO has run with no wet points options : 310 498 some parts of the domain, with no ocean points, as no 311 499 lon/lat values 312 500 ''' 313 mmath = __mmath__ (lon) 501 mmath = __mmath__ (pbounds_lon) 502 503 z_bounds_lon = np.empty ( pbounds_lon.shape ) 504 z_bounds_lat = np.empty ( pbounds_lat.shape ) 314 505 315 506 imp = SimpleImputer (missing_values=sval, strategy='mean') 316 imp.fit (lon) 317 plon = imp.transform (lon) 318 imp.fit (lat.T) 319 plat = imp.transform (lat.T).T 507 508 for n in np.arange (4) : 509 imp.fit (pbounds_lon[:,:,n]) 510 z_bounds_lon[:,:,n] = imp.transform (pbounds_lon[:,:,n]) 511 imp.fit (pbounds_lat[:,:,n].T) 512 z_bounds_lat[:,:,n] = imp.transform (pbounds_lat[:,:,n].T).T 320 513 321 514 if mmath == xr : 322 plon = xr.DataArray (plon, dims=lon.dims, coords=lon.coords) 323 plat = xr.DataArray (plat, dims=lat.dims, coords=lat.coords) 324 plon.attrs = lon.attrs ; plat.attrs = lat.attrs 325 326 plon = fixed_lon (plon) 327 328 return plon, plat 329 330 #@numba.jit(forceobj=True) 331 def jeq (lat) : 332 ''' 333 Returns j index of equator in the grid 334 515 z_bounds_lon = xr.DataArray (pbounds_lon, dims=pbounds_lon.dims, 516 coords=pbounds_lon.coords) 517 z_bounds_lat = xr.DataArray (pbounds_lat, dims=pbounds_lat.dims, 518 coords=pbounds_lat.coords) 519 z_bounds_lon.attrs.update (pbounds_lat.attrs) 520 z_bounds_lat.attrs.update (pbounds_lat.attrs) 521 522 return z_bounds_lon, z_bounds_lat 523 524 def jeq (plat) : 525 '''Returns j index of equator in the grid 526 335 527 lat : latitudes of the grid. At least 2D. 336 528 ''' 337 mmath = __mmath__ (lat) 338 ix, ax = __findAxis__ (lat, 'x') 339 iy, ay = __findAxis__ (lat, 'y') 529 mmath = __mmath__ (plat) 530 jy = __find_axis__ (plat, 'y')[-1] 340 531 341 532 if mmath == xr : 342 jeq = int ( np.mean ( np.argmin (np.abs (np.float64 (lat)), axis=iy) ) ) 343 else : 344 jeq = np.argmin (np.abs (np.float64 (lat[...,:, 0]))) 345 return jeq 346 347 #@numba.jit(forceobj=True) 348 def lon1D (lon, lat=None) : 349 ''' 350 Returns 1D longitude for simple plots. 351 352 lon : longitudes of the grid 353 lat (optionnal) : latitudes of the grid 354 ''' 355 mmath = __mmath__ (lon) 356 if np.max (lat) != None : 357 je = jeq (lat) 358 lon1D = lon.copy() [..., je, :] 359 else : 360 jpj, jpi = lon.shape [-2:] 361 lon1D = lon.copy() [..., jpj//3, :] 362 363 start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1) 364 lon1D [..., start+1:] += 360 533 jj = int ( np.mean ( np.argmin (np.abs (np.float64 (plat)), 534 axis=jy) ) ) 535 else : 536 jj = np.argmin (np.abs (np.float64 (plat[...,:, 0]))) 537 538 return jj 539 540 def lon1d (plon, plat=None) : 541 '''Returns 1D longitude for simple plots. 542 543 plon : longitudes of the grid 544 plat (optionnal) : latitudes of the grid 545 ''' 546 mmath = __mmath__ (plon) 547 jpj, jpi = plon.shape [-2:] 548 if np.max (plat) : 549 je = jeq (plat) 550 lon0 = plon [..., je, 0].copy() 551 dlon = plon [..., je, 1].copy() - plon [..., je, 0].copy() 552 lon_1d = np.linspace ( start=lon0, stop=lon0+360.+2*dlon, num=jpi ) 553 else : 554 lon0 = plon [..., jpj//3, 0].copy() 555 dlon = plon [..., jpj//3, 1].copy() - plon [..., jpj//3, 0].copy() 556 lon_1d = np.linspace ( start=lon0, stop=lon0+360.+2*dlon, num=jpi ) 557 558 #start = np.argmax (np.abs (np.diff (lon1D, axis=-1)) > 180.0, axis=-1) 559 #lon1D [..., start+1:] += 360 365 560 366 561 if mmath == xr : 367 lon1D.attrs = lon.attrs 368 lon1D = lon1D.assign_coords ( {'x':lon1D} ) 369 370 return lon1D 371 372 #@numba.jit(forceobj=True) 373 def latreg (lat, diff=0.1) : 374 ''' 375 Returns maximum j index where gridlines are along latitudes in the northern hemisphere 376 562 lon_1d = xr.DataArray( lon_1d, dims=('lon',), coords=(lon_1d,)) 563 lon_1d.attrs.update (plon.attrs) 564 lon_1d.attrs['units'] = 'degrees_east' 565 lon_1d.attrs['standard_name'] = 'longitude' 566 lon_1d.attrs['long_name :'] = 'Longitude' 567 568 return lon_1d 569 570 def latreg (plat, diff=0.1) : 571 '''Returns maximum j index where gridlines are along latitudes 572 in the northern hemisphere 573 377 574 lat : latitudes of the grid (2D) 378 575 diff [optional] : tolerance 379 576 ''' 380 mmath = __mmath__ (lat) 381 if diff == None : 382 dy = np.float64 (np.mean (np.abs (lat - np.roll (lat,shift=1,axis=-2, roll_coords=False)))) 577 #mmath = __mmath__ (plat) 578 if diff is None : 579 dy = np.float64 (np.mean (np.abs (plat - 580 np.roll (plat,shift=1,axis=-2, roll_coords=False)))) 581 print ( f'{dy=}' ) 383 582 diff = dy/100. 384 385 je = jeq (lat) 386 jreg = np.where (lat[...,je:,:].max(axis=-1) - lat[...,je:,:].min(axis=-1)< diff)[-1][-1] + je 387 latreg = np.float64 (lat[...,jreg,:].mean(axis=-1)) 388 JREG = jreg 389 390 return jreg, latreg 391 392 #@numba.jit(forceobj=True) 393 def lat1D (lat) : 394 ''' 395 Returns 1D latitudes for zonal means and simple plots. 396 397 lat : latitudes of the grid (2D) 398 ''' 399 mmath = __mmath__ (lat) 400 jpj, jpi = lat.shape[-2:] 401 402 dy = np.float64 (np.mean (np.abs (lat - np.roll (lat, shift=1,axis=-2)))) 403 je = jeq (lat) 404 lat_eq = np.float64 (lat[...,je,:].mean(axis=-1)) 405 406 jreg, lat_reg = latreg (lat) 407 lat_ave = np.mean (lat, axis=-1) 408 409 if (np.abs (lat_eq) < dy/100.) : # T, U or W grid 410 dys = (90.-lat_reg) / (jpj-jreg-1)*0.5 583 584 je = jeq (plat) 585 jreg = np.where (plat[...,je:,:].max(axis=-1) - 586 plat[...,je:,:].min(axis=-1)< diff)[-1][-1] + je 587 lareg = np.float64 (plat[...,jreg,:].mean(axis=-1)) 588 589 return jreg, lareg 590 591 def lat1d (plat) : 592 '''Returns 1D latitudes for zonal means and simple plots. 593 594 plat : latitudes of the grid (2D) 595 ''' 596 mmath = __mmath__ (plat) 597 iy = __find_axis__ (plat, 'y')[-1] 598 jpj = plat.shape[iy] 599 600 dy = np.float64 (np.mean (np.abs (plat - np.roll (plat, shift=1,axis=-2)))) 601 je = jeq (plat) 602 lat_eq = np.float64 (plat[...,je,:].mean(axis=-1)) 603 604 jreg, lat_reg = latreg (plat) 605 lat_ave = np.mean (plat, axis=-1) 606 607 if np.abs (lat_eq) < dy/100. : # T, U or W grid 608 if jpj-1 > jreg : 609 dys = (90.-lat_reg) / (jpj-jreg-1)*0.5 610 else : 611 dys = (90.-lat_reg) / 2.0 411 612 yrange = 90.-dys-lat_reg 412 613 else : # V or F grid 413 yrange = 90. -lat_reg 414 415 lat1D = mmath.where (lat_ave<lat_reg, lat_ave, lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1)) 416 614 yrange = 90.-lat_reg 615 616 if jpj-1 > jreg : 617 lat_1d = mmath.where (lat_ave<lat_reg, 618 lat_ave, 619 lat_reg + yrange * (np.arange(jpj)-jreg)/(jpj-jreg-1) ) 620 else : 621 lat_1d = lat_ave 622 lat_1d[-1] = 90.0 623 417 624 if mmath == xr : 418 lat1D.attrs = lat.attrs 419 lat1D = lat1D.assign_coords ( {'y':lat1D} ) 420 421 return lat1D 422 423 #@numba.jit(forceobj=True) 424 def latlon1D (lat, lon) : 425 ''' 426 Returns simple latitude and longitude (1D) for simple plots. 427 428 lat, lon : latitudes and longitudes of the grid (2D) 429 ''' 430 return lat1D (lat), lon1D (lon, lat) 431 432 ##@numba.jit(forceobj=True) 625 lat_1d = xr.DataArray( lat_1d.values, dims=('lat',), coords=(lat_1d,)) 626 lat_1d.attrs.update (plat.attrs) 627 lat_1d.attrs ['units'] = 'degrees_north' 628 lat_1d.attrs ['standard_name'] = 'latitude' 629 lat_1d.attrs ['long_name :'] = 'Latitude' 630 631 return lat_1d 632 633 def latlon1d (plat, plon) : 634 '''Returns simple latitude and longitude (1D) for simple plots. 635 636 plat, plon : latitudes and longitudes of the grid (2D) 637 ''' 638 return lat1d (plat), lon1d (plon, plat) 639 640 def ff (plat) : 641 '''Returns Coriolis factor 642 ''' 643 zff = np.sin (RAD * plat) * OMEGA 644 return zff 645 646 def beta (plat) : 647 '''Return Beta factor (derivative of Coriolis factor) 648 ''' 649 zbeta = np.cos (RAD * plat) * OMEGA / RA 650 return zbeta 651 433 652 def mask_lonlat (ptab, x0, x1, y0, y1, lon, lat, sval=np.nan) : 653 '''Returns masked values outside a lat/lon box 654 ''' 434 655 mmath = __mmath__ (ptab) 435 try:656 if mmath == xr : 436 657 lon = lon.copy().to_masked_array() 437 658 lat = lat.copy().to_masked_array() 438 except : pass 439 440 mask = np.logical_and (np.logical_and(lat>y0, lat<y1), 441 np.logical_or (np.logical_or (np.logical_and(lon>x0, lon<x1), np.logical_and(lon+360>x0, lon+360<x1)), 442 np.logical_and(lon-360>x0, lon-360<x1))) 443 tab = mmath.where (mask, ptab, np.nan) 444 659 660 mask = np.logical_and (np.logical_and(lat>y0, lat<y1), 661 np.logical_or (np.logical_or ( 662 np.logical_and(lon>x0, lon<x1), 663 np.logical_and(lon+360>x0, lon+360<x1)), 664 np.logical_and(lon-360>x0, lon-360<x1))) 665 tab = mmath.where (mask, ptab, sval) 666 445 667 return tab 446 668 447 #@numba.jit(forceobj=True) 448 def extend (tab, Lon=False, jplus=25, jpi=None, nperio=4) : 449 '''450 Returns extended field eastward to have better plots, and box average crossing the boundary 669 def extend (ptab, blon=False, jplus=25, jpi=None, nperio=4) : 670 '''Returns extended field eastward to have better plots, 671 and box average crossing the boundary 672 451 673 Works only for xarray and numpy data (?) 452 453 tab : field to extend. 454 Lon : (optional, default=False) : if True, add 360 in the extended parts of the field 455 jpi : normal longitude dimension of the field. exrtend does nothing it the actual 456 size of the field != jpi (avoid to extend several times) 457 jplus (optional, default=25) : number of points added on the east side of the field 458 459 ''' 460 mmath = __mmath__ (tab) 461 462 if tab.shape[-1] == 1 : extend = tab 463 464 else : 465 if jpi == None : jpi = tab.shape[-1] 466 467 if Lon : xplus = -360.0 468 else : xplus = 0.0 469 470 if tab.shape[-1] > jpi : 471 extend = tab 674 Useful for plotting vertical sections in OCE and ATM. 675 676 ptab : field to extend. 677 blon : (optional, default=False) : if True, add 360 in the extended 678 parts of the field 679 jpi : normal longitude dimension of the field. extend does nothing 680 if the actual size of the field != jpi 681 (avoid to extend several times in notebooks) 682 jplus (optional, default=25) : number of points added on 683 the east side of the field 684 685 ''' 686 mmath = __mmath__ (ptab) 687 688 if ptab.shape[-1] == 1 : 689 tabex = ptab 690 691 else : 692 if jpi is None : 693 jpi = ptab.shape[-1] 694 695 if blon : 696 xplus = -360.0 697 else : 698 xplus = 0.0 699 700 if ptab.shape[-1] > jpi : 701 tabex = ptab 472 702 else : 473 if nperio == 0 or nperio == 4.2:474 istart = 0 ; le=jpi+1 ; la=0703 if nperio in [ 0, 4.2 ] : 704 istart, le, la = 0, jpi+1, 0 475 705 if nperio == 1 : 476 istart = 0 ; le=jpi+1 ; la=0477 if nperio == 4 or nperio == 6: # OPA case with two halo points for periodicity478 istart = 1 ; le=jpi-2 ; la=1# Perfect, except at the pole that should be masked by lbc_plot479 706 istart, le, la = 0, jpi+1, 0 707 if nperio in [4, 6] : # OPA case with two halo points for periodicity 708 # Perfect, except at the pole that should be masked by lbc_plot 709 istart, le, la = 1, jpi-2, 1 480 710 if mmath == xr : 481 extend = np.concatenate ((tab.values[..., istart :istart+le+1 ] + xplus, 482 tab.values[..., istart+la:istart+la+jplus] ), axis=-1) 483 lon = tab.dims[-1] 711 tabex = np.concatenate ( 712 (ptab.values[..., istart :istart+le+1 ] + xplus, 713 ptab.values[..., istart+la:istart+la+jplus] ), 714 axis=-1) 715 lon = ptab.dims[-1] 484 716 new_coords = [] 485 for coord in tab.dims : 486 if coord == lon : new_coords.append ( np.arange( extend.shape[-1])) 487 else : new_coords.append ( tab.coords[coord].values) 488 extend = xr.DataArray ( extend, dims=tab.dims, coords=new_coords ) 489 else : 490 extend = np.concatenate ((tab [..., istart :istart+le+1 ] + xplus, 491 tab [..., istart+la:istart+la+jplus] ), axis=-1) 492 return extend 493 494 def orca2reg (ff, lat_name='nav_lat', lon_name='nav_lon', y_name='y', x_name='x') : 495 ''' 496 Assign an ORCA dataset on a regular grid. 717 for coord in ptab.dims : 718 if coord == lon : 719 new_coords.append ( np.arange( tabex.shape[-1])) 720 else : 721 new_coords.append ( ptab.coords[coord].values) 722 tabex = xr.DataArray ( tabex, dims=ptab.dims, 723 coords=new_coords ) 724 else : 725 tabex = np.concatenate ( 726 (ptab [..., istart :istart+le+1 ] + xplus, 727 ptab [..., istart+la:istart+la+jplus] ), 728 axis=-1) 729 return tabex 730 731 def orca2reg (dd, lat_name='nav_lat', lon_name='nav_lon', 732 y_name='y', x_name='x') : 733 '''Assign an ORCA dataset on a regular grid. 734 497 735 For use in the tropical region. 498 499 Inputs : 736 Inputs : 500 737 ff : xarray dataset 501 738 lat_name, lon_name : name of latitude and longitude 2D field in ff 502 739 y_name, x_name : namex of dimensions in ff 503 740 504 741 Returns : xarray dataset with rectangular grid. Incorrect above 20°N 505 742 ''' 506 743 # Compute 1D longitude and latitude 507 (lat, lon) = latlon1D (ff[lat_name], ff[lon_name]) 508 744 (zlat, zlon) = latlon1d ( dd[lat_name], dd[lon_name]) 745 746 zdd = dd 509 747 # Assign lon and lat as dimensions of the dataset 510 if y_name in ff.dims :511 lat = xr.DataArray (lat, coords=[lat,], dims=['lat',])512 ff = ff.rename_dims ({y_name: "lat",}).assign_coords (lat=lat)513 if x_name in ff.dims :514 lon = xr.DataArray (lon, coords=[lon,], dims=['lon',])515 ff = ff.rename_dims ({x_name: "lon",}).assign_coords (lon=lon)748 if y_name in zdd.dims : 749 zlat = xr.DataArray (zlat, coords=[zlat,], dims=['lat',]) 750 zdd = zdd.rename_dims ({y_name: "lat",}).assign_coords (lat=zlat) 751 if x_name in zdd.dims : 752 zlon = xr.DataArray (zlon, coords=[zlon,], dims=['lon',]) 753 zdd = zdd.rename_dims ({x_name: "lon",}).assign_coords (lon=zlon) 516 754 # Force dimensions to be in the right order 517 755 coord_order = ['lat', 'lon'] 518 756 for dim in [ 'depthw', 'depthv', 'depthu', 'deptht', 'depth', 'z', 519 'time_counter', 'time', 'tbnds', 757 'time_counter', 'time', 'tbnds', 520 758 'bnds', 'axis_nbounds', 'two2', 'two1', 'two', 'four',] : 521 if dim in ff.dims : coord_order.insert (0, dim) 522 523 ff = ff.transpose (*coord_order) 524 return ff 759 if dim in zdd.dims : 760 coord_order.insert (0, dim) 761 762 zdd = zdd.transpose (*coord_order) 763 return zdd 525 764 526 765 def lbc_init (ptab, nperio=None) : 527 ''' 528 Prepare for all lbc calls 529 766 '''Prepare for all lbc calls 767 530 768 Set periodicity on input field 531 769 nperio : Type of periodicity … … 539 777 See NEMO documentation for further details 540 778 ''' 541 jpj, jpi = ptab.shape[-2:] 542 if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio) 543 544 if nperio not in nperio_valid_range : 545 raise Exception ('nperio=', nperio, ' is not in the valid range', nperio_valid_range) 779 jpi, jpj = None, None 780 ax, ix = __find_axis__ (ptab, 'x') 781 ay, jy = __find_axis__ (ptab, 'y') 782 if ax : 783 jpi = ptab.shape[ix] 784 if ay : 785 jpj = ptab.shape[jy] 786 787 if nperio is None : 788 nperio = __guess_nperio__ (jpj, jpi, nperio) 789 790 if nperio not in NPERIO_VALID_RANGE : 791 raise AttributeError ( f'{nperio=} is not in the valid range {NPERIO_VALID_RANGE}' ) 546 792 547 793 return jpj, jpi, nperio 548 549 #@numba.jit(forceobj=True) 550 def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4U_bug=False) : 551 ''' 552 Set periodicity on input field 794 795 def lbc (ptab, nperio=None, cd_type='T', psgn=1.0, nemo_4u_bug=False) : 796 '''Set periodicity on input field 797 553 798 ptab : Input array (works for rank 2 at least : ptab[...., lat, lon]) 554 799 nperio : Type of periodicity 555 800 cd_type : Grid specification : T, U, V or F 556 801 psgn : For change of sign for vector components (1 for scalars, -1 for vector components) 557 802 558 803 See NEMO documentation for further details 559 804 ''' 560 jpj, jpi, nperio = lbc_init (ptab, nperio) 805 jpi, nperio = lbc_init (ptab, nperio)[1:] 806 ax = __find_axis__ (ptab, 'x')[0] 807 ay = __find_axis__ (ptab, 'y')[0] 561 808 psgn = ptab.dtype.type (psgn) 562 mmath = __mmath__ (ptab) 563 564 if mmath == xr : ztab = ptab.values.copy () 565 else : ztab = ptab.copy () 566 567 # 568 #> East-West boundary conditions 569 # ------------------------------ 570 if nperio in [1, 4, 6] : 809 mmath = __mmath__ (ptab) 810 811 if mmath == xr : 812 ztab = ptab.values.copy () 813 else : 814 ztab = ptab.copy () 815 816 if ax : 817 # 818 #> East-West boundary conditions 819 # ------------------------------ 820 if nperio in [1, 4, 6] : 571 821 # ... cyclic 572 ztab [..., :, 0] = ztab [..., :, -2] 573 ztab [..., :, -1] = ztab [..., :, 1] 574 # 575 #> North-South boundary conditions 576 # -------------------------------- 577 if nperio in [3, 4] : # North fold T-point pivot 578 if cd_type in [ 'T', 'W' ] : # T-, W-point 579 ztab [..., -1, 1: ] = psgn * ztab [..., -3, -1:0:-1 ] 580 ztab [..., -1, 0 ] = psgn * ztab [..., -3, 2 ] 581 ztab [..., -2, jpi//2: ] = psgn * ztab [..., -2, jpi//2:0:-1 ] 582 583 if cd_type == 'U' : 584 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -3, -1:0:-1 ] 585 ztab [..., -1, 0 ] = psgn * ztab [..., -3, 1 ] 586 ztab [..., -1, -1 ] = psgn * ztab [..., -3, -2 ] 587 588 if nemo_4U_bug : 589 ztab [..., -2, jpi//2+1:-1] = psgn * ztab [..., -2, jpi//2-2:0:-1] 590 ztab [..., -2, jpi//2-1 ] = psgn * ztab [..., -2, jpi//2 ] 591 else : 592 ztab [..., -2, jpi//2-1:-1] = psgn * ztab [..., -2, jpi//2:0:-1] 593 594 if cd_type == 'V' : 595 ztab [..., -2, 1: ] = psgn * ztab [..., -3, jpi-1:0:-1 ] 596 ztab [..., -1, 1: ] = psgn * ztab [..., -4, -1:0:-1 ] 597 ztab [..., -1, 0 ] = psgn * ztab [..., -4, 2 ] 598 599 if cd_type == 'F' : 600 ztab [..., -2, 0:-1 ] = psgn * ztab [..., -3, -1:0:-1 ] 601 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -4, -1:0:-1 ] 602 ztab [..., -1, 0 ] = psgn * ztab [..., -4, 1 ] 603 ztab [..., -1, -1 ] = psgn * ztab [..., -4, -2 ] 604 605 if nperio in [4.2] : # North fold T-point pivot 606 if cd_type in [ 'T', 'W' ] : # T-, W-point 607 ztab [..., -1, jpi//2: ] = psgn * ztab [..., -1, jpi//2:0:-1 ] 608 609 if cd_type == 'U' : 610 ztab [..., -1, jpi//2-1:-1] = psgn * ztab [..., -1, jpi//2:0:-1] 611 612 if cd_type == 'V' : 613 ztab [..., -1, 1: ] = psgn * ztab [..., -2, jpi-1:0:-1 ] 614 615 if cd_type == 'F' : 616 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -2, -1:0:-1 ] 617 618 if nperio in [5, 6] : # North fold F-point pivot 619 if cd_type in ['T', 'W'] : 620 ztab [..., -1, 0: ] = psgn * ztab [..., -2, -1::-1 ] 621 622 if cd_type == 'U' : 623 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -2, -2::-1 ] 624 ztab [..., -1, -1 ] = psgn * ztab [..., -2, 0 ] # Bug ? 625 626 if cd_type == 'V' : 627 ztab [..., -1, 0: ] = psgn * ztab [..., -3, -1::-1 ] 628 ztab [..., -2, jpi//2: ] = psgn * ztab [..., -2, jpi//2-1::-1 ] 629 630 if cd_type == 'F' : 631 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -3, -2::-1 ] 632 ztab [..., -1, -1 ] = psgn * ztab [..., -3, 0 ] 633 ztab [..., -2, jpi//2:-1] = psgn * ztab [..., -2, jpi//2-2::-1 ] 634 635 # 636 #> East-West boundary conditions 637 # ------------------------------ 638 if nperio in [1, 4, 6] : 639 # ... cyclic 640 ztab [..., :, 0] = ztab [..., :, -2] 641 ztab [..., :, -1] = ztab [..., :, 1] 822 ztab [..., 0] = ztab [..., -2] 823 ztab [..., -1] = ztab [..., 1] 824 825 if ay : 826 # 827 #> North-South boundary conditions 828 # -------------------------------- 829 if nperio in [3, 4] : # North fold T-point pivot 830 if cd_type in [ 'T', 'W' ] : # T-, W-point 831 ztab [..., -1, 1: ] = psgn * ztab [..., -3, -1:0:-1 ] 832 ztab [..., -1, 0 ] = psgn * ztab [..., -3, 2 ] 833 ztab [..., -2, jpi//2: ] = psgn * ztab [..., -2, jpi//2:0:-1 ] 834 835 if cd_type == 'U' : 836 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -3, -1:0:-1 ] 837 ztab [..., -1, 0 ] = psgn * ztab [..., -3, 1 ] 838 ztab [..., -1, -1 ] = psgn * ztab [..., -3, -2 ] 839 840 if nemo_4u_bug : 841 ztab [..., -2, jpi//2+1:-1] = psgn * ztab [..., -2, jpi//2-2:0:-1] 842 ztab [..., -2, jpi//2-1 ] = psgn * ztab [..., -2, jpi//2 ] 843 else : 844 ztab [..., -2, jpi//2-1:-1] = psgn * ztab [..., -2, jpi//2:0:-1] 845 846 if cd_type == 'V' : 847 ztab [..., -2, 1: ] = psgn * ztab [..., -3, jpi-1:0:-1 ] 848 ztab [..., -1, 1: ] = psgn * ztab [..., -4, -1:0:-1 ] 849 ztab [..., -1, 0 ] = psgn * ztab [..., -4, 2 ] 850 851 if cd_type == 'F' : 852 ztab [..., -2, 0:-1 ] = psgn * ztab [..., -3, -1:0:-1 ] 853 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -4, -1:0:-1 ] 854 ztab [..., -1, 0 ] = psgn * ztab [..., -4, 1 ] 855 ztab [..., -1, -1 ] = psgn * ztab [..., -4, -2 ] 856 857 if nperio in [4.2] : # North fold T-point pivot 858 if cd_type in [ 'T', 'W' ] : # T-, W-point 859 ztab [..., -1, jpi//2: ] = psgn * ztab [..., -1, jpi//2:0:-1 ] 860 861 if cd_type == 'U' : 862 ztab [..., -1, jpi//2-1:-1] = psgn * ztab [..., -1, jpi//2:0:-1] 863 864 if cd_type == 'V' : 865 ztab [..., -1, 1: ] = psgn * ztab [..., -2, jpi-1:0:-1 ] 866 867 if cd_type == 'F' : 868 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -2, -1:0:-1 ] 869 870 if nperio in [5, 6] : # North fold F-point pivot 871 if cd_type in ['T', 'W'] : 872 ztab [..., -1, 0: ] = psgn * ztab [..., -2, -1::-1 ] 873 874 if cd_type == 'U' : 875 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -2, -2::-1 ] 876 ztab [..., -1, -1 ] = psgn * ztab [..., -2, 0 ] # Bug ? 877 878 if cd_type == 'V' : 879 ztab [..., -1, 0: ] = psgn * ztab [..., -3, -1::-1 ] 880 ztab [..., -2, jpi//2: ] = psgn * ztab [..., -2, jpi//2-1::-1 ] 881 882 if cd_type == 'F' : 883 ztab [..., -1, 0:-1 ] = psgn * ztab [..., -3, -2::-1 ] 884 ztab [..., -1, -1 ] = psgn * ztab [..., -3, 0 ] 885 ztab [..., -2, jpi//2:-1] = psgn * ztab [..., -2, jpi//2-2::-1 ] 886 887 # 888 #> East-West boundary conditions 889 # ------------------------------ 890 if nperio in [1, 4, 6] : 891 # ... cyclic 892 ztab [..., 0] = ztab [..., -2] 893 ztab [..., -1] = ztab [..., 1] 642 894 643 895 if mmath == xr : 644 ztab = xr.DataArray ( ztab, dims=ptab.dims, coords=ptab.coords )896 ztab = xr.DataArray ( ztab, dims=ptab.dims, coords=ptab.coords ) 645 897 ztab.attrs = ptab.attrs 646 898 647 899 return ztab 648 900 649 #@numba.jit(forceobj=True)650 901 def lbc_mask (ptab, nperio=None, cd_type='T', sval=np.nan) : 651 # 652 ''' 653 Mask fields on duplicated points 902 '''Mask fields on duplicated points 903 654 904 ptab : Input array. Rank 2 at least : ptab [...., lat, lon] 655 905 nperio : Type of periodicity 656 906 cd_type : Grid specification : T, U, V or F 657 907 658 908 See NEMO documentation for further details 659 909 ''' 660 jpj, jpi, nperio = lbc_init (ptab, nperio) 910 jpi, nperio = lbc_init (ptab, nperio)[1:] 911 ax = __find_axis__ (ptab, 'x')[0] 912 ay = __find_axis__ (ptab, 'y')[0] 661 913 ztab = ptab.copy () 662 914 663 # 664 #> East-West boundary conditions 665 # ------------------------------ 666 if nperio in [1, 4, 6] : 915 if ax : 916 # 917 #> East-West boundary conditions 918 # ------------------------------ 919 if nperio in [1, 4, 6] : 667 920 # ... cyclic 668 ztab [..., :, 0] = sval 669 ztab [..., :, -1] = sval 670 671 # 672 #> South (in which nperio cases ?) 673 # -------------------------------- 674 if nperio in [1, 3, 4, 5, 6] : 675 ztab [..., 0, :] = sval 676 677 # 678 #> North-South boundary conditions 679 # -------------------------------- 680 if nperio in [3, 4] : # North fold T-point pivot 681 if cd_type in [ 'T', 'W' ] : # T-, W-point 682 ztab [..., -1, : ] = sval 683 ztab [..., -2, :jpi//2 ] = sval 684 685 if cd_type == 'U' : 686 ztab [..., -1, : ] = sval 687 ztab [..., -2, jpi//2+1: ] = sval 688 689 if cd_type == 'V' : 690 ztab [..., -2, : ] = sval 691 ztab [..., -1, : ] = sval 692 693 if cd_type == 'F' : 694 ztab [..., -2, : ] = sval 695 ztab [..., -1, : ] = sval 696 697 if nperio in [4.2] : # North fold T-point pivot 698 if cd_type in [ 'T', 'W' ] : # T-, W-point 699 ztab [..., -1, jpi//2 : ] = sval 700 701 if cd_type == 'U' : 702 ztab [..., -1, jpi//2-1:-1] = sval 703 704 if cd_type == 'V' : 705 ztab [..., -1, 1: ] = sval 706 707 if cd_type == 'F' : 708 ztab [..., -1, 0:-1 ] = sval 709 710 if nperio in [5, 6] : # North fold F-point pivot 711 if cd_type in ['T', 'W'] : 712 ztab [..., -1, 0: ] = sval 713 714 if cd_type == 'U' : 715 ztab [..., -1, 0:-1 ] = sval 716 ztab [..., -1, -1 ] = sval 717 718 if cd_type == 'V' : 719 ztab [..., -1, 0: ] = sval 720 ztab [..., -2, jpi//2: ] = sval 721 722 if cd_type == 'F' : 723 ztab [..., -1, 0:-1 ] = sval 724 ztab [..., -1, -1 ] = sval 725 ztab [..., -2, jpi//2+1:-1] = sval 921 ztab [..., 0] = sval 922 ztab [..., -1] = sval 923 924 if ay : 925 # 926 #> South (in which nperio cases ?) 927 # -------------------------------- 928 if nperio in [1, 3, 4, 5, 6] : 929 ztab [..., 0, :] = sval 930 931 # 932 #> North-South boundary conditions 933 # -------------------------------- 934 if nperio in [3, 4] : # North fold T-point pivot 935 if cd_type in [ 'T', 'W' ] : # T-, W-point 936 ztab [..., -1, : ] = sval 937 ztab [..., -2, :jpi//2 ] = sval 938 939 if cd_type == 'U' : 940 ztab [..., -1, : ] = sval 941 ztab [..., -2, jpi//2+1: ] = sval 942 943 if cd_type == 'V' : 944 ztab [..., -2, : ] = sval 945 ztab [..., -1, : ] = sval 946 947 if cd_type == 'F' : 948 ztab [..., -2, : ] = sval 949 ztab [..., -1, : ] = sval 950 951 if nperio in [4.2] : # North fold T-point pivot 952 if cd_type in [ 'T', 'W' ] : # T-, W-point 953 ztab [..., -1, jpi//2 : ] = sval 954 955 if cd_type == 'U' : 956 ztab [..., -1, jpi//2-1:-1] = sval 957 958 if cd_type == 'V' : 959 ztab [..., -1, 1: ] = sval 960 961 if cd_type == 'F' : 962 ztab [..., -1, 0:-1 ] = sval 963 964 if nperio in [5, 6] : # North fold F-point pivot 965 if cd_type in ['T', 'W'] : 966 ztab [..., -1, 0: ] = sval 967 968 if cd_type == 'U' : 969 ztab [..., -1, 0:-1 ] = sval 970 ztab [..., -1, -1 ] = sval 971 972 if cd_type == 'V' : 973 ztab [..., -1, 0: ] = sval 974 ztab [..., -2, jpi//2: ] = sval 975 976 if cd_type == 'F' : 977 ztab [..., -1, 0:-1 ] = sval 978 ztab [..., -1, -1 ] = sval 979 ztab [..., -2, jpi//2+1:-1] = sval 726 980 727 981 return ztab 728 982 729 #@numba.jit(forceobj=True)730 983 def lbc_plot (ptab, nperio=None, cd_type='T', psgn=1.0, sval=np.nan) : 731 ''' 732 Set periodicity on input field, adapted for plotting for any cartopy projection 984 '''Set periodicity on input field, for plotting for any cartopy projection 985 986 Points at the north fold are masked 987 Points for zonal periodicity are kept 733 988 ptab : Input array. Rank 2 at least : ptab[...., lat, lon] 734 989 nperio : Type of periodicity 735 990 cd_type : Grid specification : T, U, V or F 736 psgn : For change of sign for vector components (1 for scalars, -1 for vector components) 737 991 psgn : For change of sign for vector components 992 (1 for scalars, -1 for vector components) 993 738 994 See NEMO documentation for further details 739 995 ''' 740 741 jpj, jpi, nperio = lbc_init (ptab, nperio) 996 jpi, nperio = lbc_init (ptab, nperio)[1:] 997 ax = __find_axis__ (ptab, 'x')[0] 998 ay = __find_axis__ (ptab, 'y')[0] 742 999 psgn = ptab.dtype.type (psgn) 743 1000 ztab = ptab.copy () 744 # 745 #> East-West boundary conditions 746 # ------------------------------ 747 if nperio in [1, 4, 6] : 748 # ... cyclic 749 ztab [..., :, 0] = ztab [..., :, -2] 750 ztab [..., :, -1] = ztab [..., :, 1] 751 752 #> Masks south 753 # ------------ 754 if nperio in [4, 6] : ztab [..., 0, : ] = sval 755 756 # 757 #> North-South boundary conditions 758 # -------------------------------- 759 if nperio in [3, 4] : # North fold T-point pivot 760 if cd_type in [ 'T', 'W' ] : # T-, W-point 761 ztab [..., -1, : ] = sval 762 #ztab [..., -2, jpi//2: ] = sval 763 ztab [..., -2, :jpi//2 ] = sval # Give better plots than above 764 if cd_type == 'U' : 765 ztab [..., -1, : ] = sval 766 767 if cd_type == 'V' : 768 ztab [..., -2, : ] = sval 769 ztab [..., -1, : ] = sval 770 771 if cd_type == 'F' : 772 ztab [..., -2, : ] = sval 773 ztab [..., -1, : ] = sval 774 775 if nperio in [4.2] : # North fold T-point pivot 776 if cd_type in [ 'T', 'W' ] : # T-, W-point 777 ztab [..., -1, jpi//2: ] = sval 778 779 if cd_type == 'U' : 780 ztab [..., -1, jpi//2-1:-1] = sval 781 782 if cd_type == 'V' : 783 ztab [..., -1, 1: ] = sval 784 785 if cd_type == 'F' : 786 ztab [..., -1, 0:-1 ] = sval 787 788 if nperio in [5, 6] : # North fold F-point pivot 789 if cd_type in ['T', 'W'] : 790 ztab [..., -1, : ] = sval 791 792 if cd_type == 'U' : 793 ztab [..., -1, : ] = sval 794 795 if cd_type == 'V' : 796 ztab [..., -1, : ] = sval 797 ztab [..., -2, jpi//2: ] = sval 798 799 if cd_type == 'F' : 800 ztab [..., -1, : ] = sval 801 ztab [..., -2, jpi//2+1:-1] = sval 1001 1002 if ax : 1003 # 1004 #> East-West boundary conditions 1005 # ------------------------------ 1006 if nperio in [1, 4, 6] : 1007 # ... cyclic 1008 ztab [..., :, 0] = ztab [..., :, -2] 1009 ztab [..., :, -1] = ztab [..., :, 1] 1010 1011 if ay : 1012 #> Masks south 1013 # ------------ 1014 if nperio in [4, 6] : 1015 ztab [..., 0, : ] = sval 1016 1017 # 1018 #> North-South boundary conditions 1019 # -------------------------------- 1020 if nperio in [3, 4] : # North fold T-point pivot 1021 if cd_type in [ 'T', 'W' ] : # T-, W-point 1022 ztab [..., -1, : ] = sval 1023 #ztab [..., -2, jpi//2: ] = sval 1024 ztab [..., -2, :jpi//2 ] = sval # Give better plots than above 1025 if cd_type == 'U' : 1026 ztab [..., -1, : ] = sval 1027 1028 if cd_type == 'V' : 1029 ztab [..., -2, : ] = sval 1030 ztab [..., -1, : ] = sval 1031 1032 if cd_type == 'F' : 1033 ztab [..., -2, : ] = sval 1034 ztab [..., -1, : ] = sval 1035 1036 if nperio in [4.2] : # North fold T-point pivot 1037 if cd_type in [ 'T', 'W' ] : # T-, W-point 1038 ztab [..., -1, jpi//2: ] = sval 1039 1040 if cd_type == 'U' : 1041 ztab [..., -1, jpi//2-1:-1] = sval 1042 1043 if cd_type == 'V' : 1044 ztab [..., -1, 1: ] = sval 1045 1046 if cd_type == 'F' : 1047 ztab [..., -1, 0:-1 ] = sval 1048 1049 if nperio in [5, 6] : # North fold F-point pivot 1050 if cd_type in ['T', 'W'] : 1051 ztab [..., -1, : ] = sval 1052 1053 if cd_type == 'U' : 1054 ztab [..., -1, : ] = sval 1055 1056 if cd_type == 'V' : 1057 ztab [..., -1, : ] = sval 1058 ztab [..., -2, jpi//2: ] = sval 1059 1060 if cd_type == 'F' : 1061 ztab [..., -1, : ] = sval 1062 ztab [..., -2, jpi//2+1:-1] = sval 802 1063 803 1064 return ztab 804 1065 805 #@numba.jit(forceobj=True) 806 def lbc_add (ptab, nperio=None, cd_type=None, psgn=1, sval=None) : 807 ''' 808 Handle NEMO domain changes between NEMO 4.0 to NEMO 4.2 809 Peridodicity halo has been removed 1066 def lbc_add (ptab, nperio=None, cd_type=None, psgn=1) : 1067 '''Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2 1068 1069 Periodicity and north fold halos has been removed in NEMO 4.2 810 1070 This routine adds the halos if needed 811 1071 812 ptab : Input array (works 1072 ptab : Input array (works 813 1073 rank 2 at least : ptab[...., lat, lon] 814 1074 nperio : Type of periodicity 815 1075 816 1076 See NEMO documentation for further details 817 1077 ''' 818 mmath = __mmath__ (ptab) 819 jpj, jpi, nperio = lbc_init (ptab, nperio) 1078 mmath = __mmath__ (ptab) 1079 nperio = lbc_init (ptab, nperio)[-1] 1080 lshape = get_shape (ptab) 1081 ix = __find_axis__ (ptab, 'x')[-1] 1082 jy = __find_axis__ (ptab, 'y')[-1] 820 1083 821 1084 t_shape = np.array (ptab.shape) 822 1085 823 if nperio == 4.2 or nperio == 6.2 : 824 825 ext_shape = t_shape 826 ext_shape[-1] = ext_shape[-1] + 2 827 ext_shape[-2] = ext_shape[-2] + 1 1086 if nperio in [4.2, 6.2] : 1087 1088 ext_shape = t_shape.copy() 1089 if 'X' in lshape : 1090 ext_shape[ix] = ext_shape[ix] + 2 1091 if 'Y' in lshape : 1092 ext_shape[jy] = ext_shape[jy] + 1 828 1093 829 1094 if mmath == xr : 830 ptab_ext = xr.DataArray (np.zeros (ext_shape), dims=ptab.dims) 831 ptab_ext.values[..., :-1, 1:-1] = ptab.values.copy () 1095 ptab_ext = xr.DataArray (np.zeros (ext_shape), dims=ptab.dims) 1096 if 'X' in lshape and 'Y' in lshape : 1097 ptab_ext.values[..., :-1, 1:-1] = ptab.values.copy () 1098 else : 1099 if 'X' in lshape : 1100 ptab_ext.values[..., 1:-1] = ptab.values.copy () 1101 if 'Y' in lshape : 1102 ptab_ext.values[..., :-1 ] = ptab.values.copy () 832 1103 else : 833 1104 ptab_ext = np.zeros (ext_shape) 834 ptab_ext[..., :-1, 1:-1] = ptab.copy () 835 836 #if sval != None : ptab_ext[..., 0, :] = sval 837 838 if nperio == 4.2 : ptab_ext = lbc (ptab_ext, nperio=4, cd_type=cd_type, psgn=psgn) 839 if nperio == 6.2 : ptab_ext = lbc (ptab_ext, nperio=6, cd_type=cd_type, psgn=psgn) 840 1105 if 'X' in lshape and 'Y' in lshape : 1106 ptab_ext [..., :-1, 1:-1] = ptab.copy () 1107 else : 1108 if 'X' in lshape : 1109 ptab_ext [..., 1:-1] = ptab.copy () 1110 if 'Y' in lshape : 1111 ptab_ext [..., :-1 ] = ptab.copy () 1112 1113 if nperio == 4.2 : 1114 ptab_ext = lbc (ptab_ext, nperio=4, cd_type=cd_type, psgn=psgn) 1115 if nperio == 6.2 : 1116 ptab_ext = lbc (ptab_ext, nperio=6, cd_type=cd_type, psgn=psgn) 1117 841 1118 if mmath == xr : 842 1119 ptab_ext.attrs = ptab.attrs 1120 az = __find_axis__ (ptab, 'z')[0] 1121 at = __find_axis__ (ptab, 't')[0] 1122 if az : 1123 ptab_ext = ptab_ext.assign_coords ( {az:ptab.coords[az]} ) 1124 if at : 1125 ptab_ext = ptab_ext.assign_coords ( {at:ptab.coords[at]} ) 843 1126 844 1127 else : ptab_ext = lbc (ptab, nperio=nperio, cd_type=cd_type, psgn=psgn) 845 1128 846 1129 return ptab_ext 847 1130 848 1131 def lbc_del (ptab, nperio=None, cd_type='T', psgn=1) : 849 ''' 850 Handle NEMO domain changes between NEMO 4.0 to NEMO 4.2 851 Periodicity halo has been removed1132 '''Handles NEMO domain changes between NEMO 4.0 to NEMO 4.2 1133 1134 Periodicity and north fold halos has been removed in NEMO 4.2 852 1135 This routine removes the halos if needed 853 1136 854 ptab : Input array (works 1137 ptab : Input array (works 855 1138 rank 2 at least : ptab[...., lat, lon] 856 1139 nperio : Type of periodicity 857 1140 858 1141 See NEMO documentation for further details 859 1142 ''' 860 861 jpj, jpi, nperio = lbc_init (ptab, nperio) 862 863 if nperio == 4.2 or nperio == 6.2 : 864 return lbc (ptab[..., :-1, 1:-1], nperio=nperio, cd_type=cd_type, psgn=psgn) 865 else : 866 return ptab 867 868 #@numba.jit(forceobj=True) 1143 nperio = lbc_init (ptab, nperio)[-1] 1144 #lshape = get_shape (ptab) 1145 ax = __find_axis__ (ptab, 'x')[0] 1146 ay = __find_axis__ (ptab, 'y')[0] 1147 1148 if nperio in [4.2, 6.2] : 1149 if ax or ay : 1150 if ax and ay : 1151 ztab = lbc (ptab[..., :-1, 1:-1], 1152 nperio=nperio, cd_type=cd_type, psgn=psgn) 1153 else : 1154 if ax : 1155 ztab = lbc (ptab[..., 1:-1], 1156 nperio=nperio, cd_type=cd_type, psgn=psgn) 1157 if ay : 1158 ztab = lbc (ptab[..., -1], 1159 nperio=nperio, cd_type=cd_type, psgn=psgn) 1160 else : 1161 ztab = ptab 1162 else : 1163 ztab = ptab 1164 1165 return ztab 1166 869 1167 def lbc_index (jj, ii, jpj, jpi, nperio=None, cd_type='T') : 870 ''' 871 For indexes of a NEMO point, give the corresponding point inside the util domain 1168 '''For indexes of a NEMO point, give the corresponding point 1169 inside the domain (i.e. not in the halo) 1170 872 1171 jj, ii : indexes 873 1172 jpi, jpi : size of domain 874 1173 nperio : type of periodicity 875 1174 cd_type : grid specification : T, U, V or F 876 1175 877 1176 See NEMO documentation for further details 878 1177 ''' 879 1178 880 if nperio == None : nperio = __guessNperio__ (jpj, jpi, nperio) 881 882 ## For the sake of simplicity, switch to the convention of original lbc Fortran routine from NEMO 883 ## : starts indexes at 1 884 jy = jj + 1 ; ix = ii + 1 1179 if nperio is None : 1180 nperio = __guess_nperio__ (jpj, jpi, nperio) 1181 1182 ## For the sake of simplicity, switch to the convention of original 1183 ## lbc Fortran routine from NEMO : starts indexes at 1 1184 jy = jj + 1 1185 ix = ii + 1 885 1186 886 1187 mmath = __mmath__ (jj) 887 if mmath == None : mmath=np 1188 if mmath is None : 1189 mmath=np 888 1190 889 1191 # … … 896 1198 897 1199 # 898 def mod IJ(cond, jy_new, ix_new) :1200 def mod_ij (cond, jy_new, ix_new) : 899 1201 jy_r = mmath.where (cond, jy_new, jy) 900 1202 ix_r = mmath.where (cond, ix_new, ix) … … 905 1207 if nperio in [ 3 , 4 ] : 906 1208 if cd_type in [ 'T' , 'W' ] : 907 (jy, ix) = modIJ (np.logical_and (jy==jpj , ix>=2 ), jpj-2, jpi-ix+2) 908 (jy, ix) = modIJ (np.logical_and (jy==jpj , ix==1 ), jpj-1, 3 ) 909 (jy, ix) = modIJ (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy , jpi-ix+2) 1209 jy, ix = mod_ij (np.logical_and (jy==jpj , ix>=2 ), jpj-2, jpi-ix+2) 1210 jy, ix = mod_ij (np.logical_and (jy==jpj , ix==1 ), jpj-1, 3 ) 1211 jy, ix = mod_ij (np.logical_and (jy==jpj-1, ix>=jpi//2+1), 1212 jy , jpi-ix+2) 910 1213 911 1214 if cd_type in [ 'U' ] : 912 (jy, ix) = modIJ (np.logical_and (jy==jpj , np.logical_and (ix>=1, ix <= jpi-1) ), jy , jpi-ix+1) 913 (jy, ix) = modIJ (np.logical_and (jy==jpj , ix==1 ) , jpj-2, 2 ) 914 (jy, ix) = modIJ (np.logical_and (jy==jpj , ix==jpi) , jpj-2, jpi-1 ) 915 (jy, ix) = modIJ (np.logical_and (jy==jpj-1, np.logical_and (ix>=jpi//2, ix<=jpi-1)), jy , jpi-ix+1) 916 1215 jy, ix = mod_ij (np.logical_and ( 1216 jy==jpj , 1217 np.logical_and (ix>=1, ix <= jpi-1) ), 1218 jy , jpi-ix+1) 1219 jy, ix = mod_ij (np.logical_and (jy==jpj , ix==1 ) , jpj-2, 2 ) 1220 jy, ix = mod_ij (np.logical_and (jy==jpj , ix==jpi) , jpj-2, jpi-1 ) 1221 jy, ix = mod_ij (np.logical_and (jy==jpj-1, 1222 np.logical_and (ix>=jpi//2, ix<=jpi-1)), jy , jpi-ix+1) 1223 917 1224 if cd_type in [ 'V' ] : 918 (jy, ix) = modIJ(np.logical_and (jy==jpj-1, ix>=2 ), jpj-2, jpi-ix+2)919 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix>=2 ), jpj-3, jpi-ix+2)920 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix==1 ), jpj-3, 3 )921 1225 jy, ix = mod_ij (np.logical_and (jy==jpj-1, ix>=2 ), jpj-2, jpi-ix+2) 1226 jy, ix = mod_ij (np.logical_and (jy==jpj , ix>=2 ), jpj-3, jpi-ix+2) 1227 jy, ix = mod_ij (np.logical_and (jy==jpj , ix==1 ), jpj-3, 3 ) 1228 922 1229 if cd_type in [ 'F' ] : 923 (jy, ix) = modIJ(np.logical_and (jy==jpj-1, ix<=jpi-1), jpj-2, jpi-ix+1)924 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix<=jpi-1), jpj-3, jpi-ix+1)925 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix==1 ), jpj-3, 2 )926 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix==jpi ), jpj-3, jpi-1 )1230 jy, ix = mod_ij (np.logical_and (jy==jpj-1, ix<=jpi-1), jpj-2, jpi-ix+1) 1231 jy, ix = mod_ij (np.logical_and (jy==jpj , ix<=jpi-1), jpj-3, jpi-ix+1) 1232 jy, ix = mod_ij (np.logical_and (jy==jpj , ix==1 ), jpj-3, 2 ) 1233 jy, ix = mod_ij (np.logical_and (jy==jpj , ix==jpi ), jpj-3, jpi-1 ) 927 1234 928 1235 if nperio in [ 5 , 6 ] : 929 1236 if cd_type in [ 'T' , 'W' ] : # T-, W-point 930 (jy, ix) = modIJ(jy==jpj, jpj-1, jpi-ix+1)931 1237 jy, ix = mod_ij (jy==jpj, jpj-1, jpi-ix+1) 1238 932 1239 if cd_type in [ 'U' ] : # U-point 933 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix<=jpi-1 ), jpj-1, jpi-ix )934 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix==jpi ), jpi-1, 1 )935 1240 jy, ix = mod_ij (np.logical_and (jy==jpj , ix<=jpi-1 ), jpj-1, jpi-ix ) 1241 jy, ix = mod_ij (np.logical_and (jy==jpj , ix==jpi ), jpi-1, 1 ) 1242 936 1243 if cd_type in [ 'V' ] : # V-point 937 (jy, ix) = modIJ(jy==jpj , jy , jpi-ix+1)938 (jy, ix) = modIJ(np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy , jpi-ix+1)939 1244 jy, ix = mod_ij (jy==jpj , jy , jpi-ix+1) 1245 jy, ix = mod_ij (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy , jpi-ix+1) 1246 940 1247 if cd_type in [ 'F' ] : # F-point 941 (jy, ix) = modIJ(np.logical_and (jy==jpj , ix<=jpi-1 ), jpj-2, jpi-ix )942 (jy, ix) = modIJ(np.logical_and (ix==jpj , ix==jpi ), jpj-2, 1 )943 (jy, ix) = modIJ(np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy , jpi-ix )1248 jy, ix = mod_ij (np.logical_and (jy==jpj , ix<=jpi-1 ), jpj-2, jpi-ix ) 1249 jy, ix = mod_ij (np.logical_and (ix==jpj , ix==jpi ), jpj-2, 1 ) 1250 jy, ix = mod_ij (np.logical_and (jy==jpj-1, ix>=jpi//2+1), jy , jpi-ix ) 944 1251 945 1252 ## Restore convention to Python/C : indexes start at 0 946 jy += -1 ; ix += -1 947 948 if isinstance (jj, int) : jy = jy.item () 949 if isinstance (ii, int) : ix = ix.item () 1253 jy += -1 1254 ix += -1 1255 1256 if isinstance (jj, int) : 1257 jy = jy.item () 1258 if isinstance (ii, int) : 1259 ix = ix.item () 950 1260 951 1261 return jy, ix 952 953 def geo2en (pxx, pyy, pzz, glam, gphi) : 954 ''' 955 Change vector from geocentric to east/north 1262 1263 def find_ji (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False, out=None) : 1264 ''' 1265 Description: seeks J,I indices of the grid point which is the closest 1266 of a given point 1267 1268 Usage: go FindJI <data latitude> <data longitude> <grid latitudes> <grid longitudes> [mask] 1269 <grid latitudes><grid longitudes> are 2D fields on J/I (Y/X) dimensions 1270 mask : if given, seek only non masked grid points (i.e with mask=1) 1271 1272 Example : findIJ (40, -20, nav_lat, nav_lon, mask=1.0) 1273 1274 Note : all longitudes and latitudes in degrees 1275 1276 Note : may work with 1D lon/lat (?) 1277 ''' 1278 # Get grid dimensions 1279 if len (lon_grid.shape) == 2 : 1280 jpi = lon_grid.shape[-1] 1281 else : 1282 jpi = len(lon_grid) 1283 1284 #mmath = __mmath__ (lat_grid) 1285 1286 # Compute distance from the point to all grid points (in RADian) 1287 arg = ( np.sin (RAD*lat_data) * np.sin (RAD*lat_grid) 1288 + np.cos (RAD*lat_data) * np.cos (RAD*lat_grid) * 1289 np.cos(RAD*(lon_data-lon_grid)) ) 1290 # Send masked points to 'infinite' 1291 distance = np.arccos (arg) + 4.0*RPI*(1.0-mask) 1292 1293 # Truncates to alleviate some precision problem with some grids 1294 prec = int (1E7) 1295 distance = (distance*prec).astype(int) / prec 1296 1297 # Compute minimum of distance, and index of minimum 1298 # 1299 #distance_min = distance.min () 1300 jimin = int (distance.argmin ()) 1301 1302 # Compute 2D indices (Python/C flavor : starting at 0) 1303 jmin = jimin // jpi 1304 imin = jimin - jmin*jpi 1305 1306 # Result 1307 if verbose : 1308 # Compute distance achieved 1309 #mindist = distance [jmin, imin] 1310 1311 # Compute azimuth 1312 dlon = lon_data-lon_grid[jmin,imin] 1313 arg = np.sin (RAD*dlon) / ( 1314 np.cos(RAD*lat_data)*np.tan(RAD*lat_grid[jmin,imin]) 1315 - np.sin(RAD*lat_data)*np.cos(RAD*dlon)) 1316 azimuth = DAR*np.arctan (arg) 1317 print ( f'I={imin:d} J={jmin:d} - Data:{lat_data:5.1f}°N {lon_data:5.1f}°E' \ 1318 + f'- Grid:{lat_grid[jmin,imin]:4.1f}°N ' \ 1319 + f'{lon_grid[jmin,imin]:4.1f}°E - Dist: {RA*distance[jmin,imin]:6.1f}km' \ 1320 + f' {DAR*distance[jmin,imin]:5.2f}° ' \ 1321 + f'- Azimuth: {RAD*azimuth:3.2f}RAD - {azimuth:5.1f}°' ) 1322 1323 if out=='dict' : 1324 return {'x':imin, 'y':jmin} 1325 elif out in ['array', 'numpy', 'np'] : 1326 return np.array ( [jmin, imin] ) 1327 elif out in ['xarray', 'xr'] : 1328 return xr.DataArray ( [jmin, imin] ) 1329 elif out=='list' : 1330 return [jmin, imin] 1331 elif out=='tuple' : 1332 return jmin, imin 1333 else : 1334 return jmin, imin 1335 1336 def curl (tx, ty, e1f, e2f, nperio=None) : 1337 '''Returns curl of a vector field 1338 ''' 1339 ax = __find_axis__ (tx, 'x')[0] 1340 ay = __find_axis__ (ty, 'y')[0] 1341 1342 tx_0 = lbc_add (tx , nperio=nperio, cd_type='U', psgn=-1) 1343 ty_0 = lbc_add (ty , nperio=nperio, cd_type='V', psgn=-1) 1344 e1f_0 = lbc_add (e1f, nperio=nperio, cd_type='U', psgn=-1) 1345 e2f_0 = lbc_add (e2f, nperio=nperio, cd_type='V', psgn=-1) 1346 1347 tx_1 = tx_0.roll ( {ay:-1} ) 1348 ty_1 = ty_0.roll ( {ax:-1} ) 1349 tx_1 = lbc (tx_1, nperio=nperio, cd_type='U', psgn=-1) 1350 ty_1 = lbc (ty_1, nperio=nperio, cd_type='V', psgn=-1) 1351 1352 zcurl = (ty_1 - ty_0)/e1f_0 - (tx_1 - tx_0)/e2f_0 1353 1354 mask = np.logical_or ( 1355 np.logical_or ( np.isnan(tx_0), np.isnan(tx_1)), 1356 np.logical_or ( np.isnan(ty_0), np.isnan(ty_1)) ) 1357 1358 zcurl = zcurl.where (np.logical_not (mask), np.nan) 1359 1360 zcurl = lbc_del (zcurl, nperio=nperio, cd_type='F', psgn=1) 1361 zcurl = lbc (zcurl, nperio=nperio, cd_type='F', psgn=1) 1362 1363 return zcurl 1364 1365 def div (ux, uy, e1t, e2t, nperio=None) : 1366 '''Returns divergence of a vector field 1367 ''' 1368 ax = __find_axis__ (ux, 'x')[0] 1369 ay = __find_axis__ (ux, 'y')[0] 1370 1371 ux_0 = lbc_add (ux , nperio=nperio, cd_type='U', psgn=-1) 1372 uy_0 = lbc_add (uy , nperio=nperio, cd_type='V', psgn=-1) 1373 e1t_0 = lbc_add (e1t, nperio=nperio, cd_type='U', psgn=-1) 1374 e2t_0 = lbc_add (e2t, nperio=nperio, cd_type='V', psgn=-1) 1375 1376 ux_1 = ux_0.roll ( {ay:1} ) 1377 uy_1 = uy_0.roll ( {ax:1} ) 1378 ux_1 = lbc (ux_1, nperio=nperio, cd_type='U', psgn=-1) 1379 uy_1 = lbc (uy_1, nperio=nperio, cd_type='V', psgn=-1) 1380 1381 zdiv = (ux_0 - ux_1)/e2t_0 + (uy_0 - uy_1)/e1t_0 1382 1383 mask = np.logical_or ( 1384 np.logical_or ( np.isnan(ux_0), np.isnan(ux_1)), 1385 np.logical_or ( np.isnan(uy_0), np.isnan(uy_1)) ) 1386 1387 zdiv = zdiv.where (np.logical_not (mask), np.nan) 1388 1389 zdiv = lbc_del (zdiv, nperio=nperio, cd_type='T', psgn=1) 1390 zdiv = lbc (zdiv, nperio=nperio, cd_type='T', psgn=1) 1391 1392 return zdiv 1393 1394 def geo2en (pxx, pyy, pzz, glam, gphi) : 1395 '''Change vector from geocentric to east/north 956 1396 957 1397 Inputs : … … 960 1400 ''' 961 1401 962 gsinlon = np.sin ( rad* glam)963 gcoslon = np.cos ( rad* glam)964 gsinlat = np.sin ( rad* gphi)965 gcoslat = np.cos ( rad* gphi)966 1402 gsinlon = np.sin (RAD * glam) 1403 gcoslon = np.cos (RAD * glam) 1404 gsinlat = np.sin (RAD * gphi) 1405 gcoslat = np.cos (RAD * gphi) 1406 967 1407 pte = - pxx * gsinlon + pyy * gcoslon 968 1408 ptn = - pxx * gcoslon * gsinlat - pyy * gsinlon * gsinlat + pzz * gcoslat … … 971 1411 972 1412 def en2geo (pte, ptn, glam, gphi) : 973 ''' 974 Change vector from east/north to geocentric 975 976 Inputs : 1413 '''Change vector from east/north to geocentric 1414 1415 Inputs : 977 1416 pte, ptn : eastward/northward components 978 1417 glam, gphi : longitude and latitude of the points 979 1418 ''' 980 981 gsinlon = np.sin ( rad* glam)982 gcoslon = np.cos ( rad* glam)983 gsinlat = np.sin ( rad* gphi)984 gcoslat = np.cos ( rad* gphi)1419 1420 gsinlon = np.sin (RAD * glam) 1421 gcoslon = np.cos (RAD * glam) 1422 gsinlat = np.sin (RAD * gphi) 1423 gcoslat = np.cos (RAD * gphi) 985 1424 986 1425 pxx = - pte * gsinlon - ptn * gcoslon * gsinlat 987 1426 pyy = pte * gcoslon - ptn * gsinlon * gsinlat 988 1427 pzz = ptn * gcoslat 989 1428 990 1429 return pxx, pyy, pzz 991 1430 992 def findJI (lat_data, lon_data, lat_grid, lon_grid, mask=1.0, verbose=False) : 993 ''' 994 Description: seeks J,I indices of the grid point which is the closest of a given point 995 Usage: go FindJI <data latitude> <data longitude> <grid latitudes> <grid longitudes> [mask] 996 <longitude fields> <latitude field> are 2D fields on J/I (Y/X) dimensions 997 mask : if given, seek only non masked grid points (i.e with mask=1) 998 999 Example : findIJ (40, -20, nav_lat, nav_lon, mask=1.0) 1000 1001 Note : all longitudes and latitudes in degrees 1002 1003 Note : may work with 1D lon/lat (?) 1004 ''' 1005 # Get grid dimensions 1006 if len (lon_grid.shape) == 2 : (jpj, jpi) = lon_grid.shape 1007 else : jpj = len(lat_grid) ; jpi=len(lon_grid) 1008 1009 mmath = __mmath__ (lat_grid) 1010 1011 # Compute distance from the point to all grid points (in radian) 1012 arg = np.sin (rad*lat_data) * np.sin (rad*lat_grid) \ 1013 + np.cos (rad*lat_data) * np.cos (rad*lat_grid) * np.cos(rad*(lon_data-lon_grid)) 1014 distance = np.arccos (arg) + 4.0*rpi*(1.0-mask) # Send masked points to 'infinite' 1015 1016 # Truncates to alleviate some precision problem with some grids 1017 prec = int (1E7) 1018 distance = (distance*prec).astype(int) / prec 1019 1020 # Compute minimum of distance, and index of minimum 1021 # 1022 distance_min = distance.min () 1023 jimin = int (distance.argmin ()) 1024 1025 # Compute 2D indices 1026 jmin = jimin // jpi ; imin = jimin - jmin*jpi 1027 1028 # Compute distance achieved 1029 mindist = distance[jmin, imin] 1030 1031 # Compute azimuth 1032 dlon = lon_data-lon_grid[jmin,imin] 1033 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)) 1034 azimuth = dar*np.arctan (arg) 1035 1036 # Result 1037 if verbose : 1038 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}°' 1039 .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)) 1040 1041 return jmin, imin 1042 1043 def clo_lon (lon, lon0) : 1044 '''Choose closest to lon0 longitude, adding or substacting 360° if needed''' 1431 1432 def clo_lon (lon, lon0=0., rad=False, deg=True) : 1433 '''Choose closest to lon0 longitude, adding/substacting 360° 1434 if needed 1435 ''' 1045 1436 mmath = __mmath__ (lon, np) 1046 1047 clo_lon = lon 1048 clo_lon = mmath.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon) 1049 clo_lon = mmath.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon) 1050 clo_lon = mmath.where (clo_lon > lon0 + 180., clo_lon-360., clo_lon) 1051 clo_lon = mmath.where (clo_lon < lon0 - 180., clo_lon+360., clo_lon) 1052 if clo_lon.shape == () : clo_lon = clo_lon.item () 1053 return clo_lon 1054 1055 def angle_full (glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif, nperio=None) : 1056 '''Compute sinus and cosinus of model line direction with respect to east''' 1437 if rad : 1438 lon_range = 2.*np.pi 1439 if deg : 1440 lon_range = 360. 1441 c_lon = lon 1442 c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, 1443 c_lon-lon_range, clo_lon) 1444 c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, 1445 c_lon+lon_range, clo_lon) 1446 c_lon = mmath.where (c_lon > lon0 + lon_range*0.5, 1447 c_lon-lon_range, clo_lon) 1448 c_lon = mmath.where (c_lon < lon0 - lon_range*0.5, 1449 c_lon+lon_range, clo_lon) 1450 if c_lon.shape == () : 1451 c_lon = c_lon.item () 1452 if mmath == xr : 1453 if lon.attrs : 1454 c_lon.attrs.update ( lon.attrs ) 1455 return c_lon 1456 1457 def index2depth (pk, gdept_0) : 1458 '''From index (real, continuous), get depth 1459 1460 Needed to use transforms in Matplotlib 1461 ''' 1462 jpk = gdept_0.shape[0] 1463 kk = xr.DataArray(pk) 1464 k = np.maximum (0, np.minimum (jpk-1, kk )) 1465 k0 = np.floor (k).astype (int) 1466 k1 = np.maximum (0, np.minimum (jpk-1, k0+1)) 1467 zz = k - k0 1468 gz = (1.0-zz)*gdept_0[k0]+ zz*gdept_0[k1] 1469 return gz.values 1470 1471 def depth2index (pz, gdept_0) : 1472 '''From depth, get index (real, continuous) 1473 1474 Needed to use transforms in Matplotlib 1475 ''' 1476 jpk = gdept_0.shape[0] 1477 if isinstance (pz, xr.core.dataarray.DataArray ) : 1478 zz = xr.DataArray (pz.values, dims=('zz',)) 1479 elif isinstance (pz, np.ndarray) : 1480 zz = xr.DataArray (pz.ravel(), dims=('zz',)) 1481 else : 1482 zz = xr.DataArray (np.array([pz]).ravel(), dims=('zz',)) 1483 zz = np.minimum (gdept_0[-1], np.maximum (0, zz)) 1484 1485 idk1 = np.minimum ( (gdept_0-zz), 0.).argmax (axis=0).astype(int) 1486 idk1 = np.maximum (0, np.minimum (jpk-1, idk1 )) 1487 idk2 = np.maximum (0, np.minimum (jpk-1, idk1-1)) 1488 1489 zff = (zz - gdept_0[idk2])/(gdept_0[idk1]-gdept_0[idk2]) 1490 idk = zff*idk1 + (1.0-zff)*idk2 1491 idk = xr.where ( np.isnan(idk), idk1, idk) 1492 return idk.values 1493 1494 def index2depth_panels (pk, gdept_0, depth0, fact) : 1495 '''From index (real, continuous), get depth, with bottom part compressed 1496 1497 Needed to use transforms in Matplotlib 1498 ''' 1499 jpk = gdept_0.shape[0] 1500 kk = xr.DataArray (pk) 1501 k = np.maximum (0, np.minimum (jpk-1, kk )) 1502 k0 = np.floor (k).astype (int) 1503 k1 = np.maximum (0, np.minimum (jpk-1, k0+1)) 1504 zz = k - k0 1505 gz = (1.0-zz)*gdept_0[k0]+ zz*gdept_0[k1] 1506 gz = xr.where ( gz<depth0, gz, depth0 + (gz-depth0)*fact) 1507 return gz.values 1508 1509 def depth2index_panels (pz, gdept_0, depth0, fact) : 1510 '''From index (real, continuous), get depth, with bottom part compressed 1511 1512 Needed to use transforms in Matplotlib 1513 ''' 1514 jpk = gdept_0.shape[0] 1515 if isinstance (pz, xr.core.dataarray.DataArray) : 1516 zz = xr.DataArray (pz.values , dims=('zz',)) 1517 elif isinstance (pz, np.ndarray) : 1518 zz = xr.DataArray (pz.ravel(), dims=('zz',)) 1519 else : 1520 zz = xr.DataArray (np.array([pz]).ravel(), dims=('zz',)) 1521 zz = np.minimum (gdept_0[-1], np.maximum (0, zz)) 1522 gdept_comp = xr.where ( gdept_0>depth0, 1523 (gdept_0-depth0)*fact+depth0, gdept_0) 1524 zz_comp = xr.where ( zz >depth0, (zz -depth0)*fact+depth0, 1525 zz ) 1526 zz_comp = np.minimum (gdept_comp[-1], np.maximum (0, zz_comp)) 1527 1528 idk1 = np.minimum ( (gdept_0-zz_comp), 0.).argmax (axis=0).astype(int) 1529 idk1 = np.maximum (0, np.minimum (jpk-1, idk1 )) 1530 idk2 = np.maximum (0, np.minimum (jpk-1, idk1-1)) 1531 1532 zff = (zz_comp - gdept_0[idk2])/(gdept_0[idk1]-gdept_0[idk2]) 1533 idk = zff*idk1 + (1.0-zff)*idk2 1534 idk = xr.where ( np.isnan(idk), idk1, idk) 1535 return idk.values 1536 1537 def depth2comp (pz, depth0, fact ) : 1538 '''Form depth, get compressed depth, with bottom part compressed 1539 1540 Needed to use transforms in Matplotlib 1541 ''' 1542 #print ('start depth2comp') 1543 if isinstance (pz, xr.core.dataarray.DataArray) : 1544 zz = pz.values 1545 elif isinstance(pz, list) : 1546 zz = np.array (pz) 1547 else : 1548 zz = pz 1549 gz = np.where ( zz>depth0, (zz-depth0)*fact+depth0, zz) 1550 #print ( f'depth2comp : {gz=}' ) 1551 if type (pz) in [int, float] : 1552 return gz.item() 1553 else : 1554 return gz 1555 1556 def comp2depth (pz, depth0, fact ) : 1557 '''Form compressed depth, get depth, with bottom part compressed 1558 1559 Needed to use transforms in Matplotlib 1560 ''' 1561 if isinstance (pz, xr.core.dataarray.DataArray) : 1562 zz = pz.values 1563 elif isinstance (pz, list) : 1564 zz = np.array (pz) 1565 else : 1566 zz = pz 1567 gz = np.where ( zz>depth0, (zz-depth0)/fact+depth0, zz) 1568 if type (pz) in [int, float] : 1569 gz = gz.item() 1570 1571 return gz 1572 1573 def index2lon (pi, plon_1d) : 1574 '''From index (real, continuous), get longitude 1575 1576 Needed to use transforms in Matplotlib 1577 ''' 1578 jpi = plon_1d.shape[0] 1579 ii = xr.DataArray (pi) 1580 i = np.maximum (0, np.minimum (jpi-1, ii )) 1581 i0 = np.floor (i).astype (int) 1582 i1 = np.maximum (0, np.minimum (jpi-1, i0+1)) 1583 xx = i - i0 1584 gx = (1.0-xx)*plon_1d[i0]+ xx*plon_1d[i1] 1585 return gx.values 1586 1587 def lon2index (px, plon_1d) : 1588 '''From longitude, get index (real, continuous) 1589 1590 Needed to use transforms in Matplotlib 1591 ''' 1592 jpi = plon_1d.shape[0] 1593 if isinstance (px, xr.core.dataarray.DataArray) : 1594 xx = xr.DataArray (px.values , dims=('xx',)) 1595 elif isinstance (px, np.ndarray) : 1596 xx = xr.DataArray (px.ravel(), dims=('xx',)) 1597 else : 1598 xx = xr.DataArray (np.array([px]).ravel(), dims=('xx',)) 1599 xx = xr.where ( xx>plon_1d.max(), xx-360.0, xx) 1600 xx = xr.where ( xx<plon_1d.min(), xx+360.0, xx) 1601 xx = np.minimum (plon_1d.max(), np.maximum(xx, plon_1d.min() )) 1602 idi1 = np.minimum ( (plon_1d-xx), 0.).argmax (axis=0).astype(int) 1603 idi1 = np.maximum (0, np.minimum (jpi-1, idi1 )) 1604 idi2 = np.maximum (0, np.minimum (jpi-1, idi1-1)) 1605 1606 zff = (xx - plon_1d[idi2])/(plon_1d[idi1]-plon_1d[idi2]) 1607 idi = zff*idi1 + (1.0-zff)*idi2 1608 idi = xr.where ( np.isnan(idi), idi1, idi) 1609 return idi.values 1610 1611 def index2lat (pj, plat_1d) : 1612 '''From index (real, continuous), get latitude 1613 1614 Needed to use transforms in Matplotlib 1615 ''' 1616 jpj = plat_1d.shape[0] 1617 jj = xr.DataArray (pj) 1618 j = np.maximum (0, np.minimum (jpj-1, jj )) 1619 j0 = np.floor (j).astype (int) 1620 j1 = np.maximum (0, np.minimum (jpj-1, j0+1)) 1621 yy = j - j0 1622 gy = (1.0-yy)*plat_1d[j0]+ yy*plat_1d[j1] 1623 return gy.values 1624 1625 def lat2index (py, plat_1d) : 1626 '''From latitude, get index (real, continuous) 1627 1628 Needed to use transforms in Matplotlib 1629 ''' 1630 jpj = plat_1d.shape[0] 1631 if isinstance (py, xr.core.dataarray.DataArray) : 1632 yy = xr.DataArray (py.values , dims=('yy',)) 1633 elif isinstance (py, np.ndarray) : 1634 yy = xr.DataArray (py.ravel(), dims=('yy',)) 1635 else : 1636 yy = xr.DataArray (np.array([py]).ravel(), dims=('yy',)) 1637 yy = np.minimum (plat_1d.max(), np.minimum(yy, plat_1d.max() )) 1638 idj1 = np.minimum ( (plat_1d-yy), 0.).argmax (axis=0).astype(int) 1639 idj1 = np.maximum (0, np.minimum (jpj-1, idj1 )) 1640 idj2 = np.maximum (0, np.minimum (jpj-1, idj1-1)) 1641 1642 zff = (yy - plat_1d[idj2])/(plat_1d[idj1]-plat_1d[idj2]) 1643 idj = zff*idj1 + (1.0-zff)*idj2 1644 idj = xr.where ( np.isnan(idj), idj1, idj) 1645 return idj.values 1646 1647 def angle_full (glamt, gphit, glamu, gphiu, glamv, gphiv, 1648 glamf, gphif, nperio=None) : 1649 '''Computes sinus and cosinus of model line direction with 1650 respect to east 1651 ''' 1057 1652 mmath = __mmath__ (glamt) 1058 1653 … … 1065 1660 zlamf = lbc_add (glamf, nperio, 'F', 1.) 1066 1661 zphif = lbc_add (gphif, nperio, 'F', 1.) 1067 1662 1068 1663 # north pole direction & modulous (at T-point) 1069 zxnpt = 0. - 2.0 * np.cos ( rad*zlamt) * np.tan (rpi/4.0 - rad*zphit/2.0)1070 zynpt = 0. - 2.0 * np.sin ( rad*zlamt) * np.tan (rpi/4.0 - rad*zphit/2.0)1664 zxnpt = 0. - 2.0 * np.cos (RAD*zlamt) * np.tan (RPI/4.0 - RAD*zphit/2.0) 1665 zynpt = 0. - 2.0 * np.sin (RAD*zlamt) * np.tan (RPI/4.0 - RAD*zphit/2.0) 1071 1666 znnpt = zxnpt*zxnpt + zynpt*zynpt 1072 1667 1073 1668 # north pole direction & modulous (at U-point) 1074 zxnpu = 0. - 2.0 * np.cos ( rad*zlamu) * np.tan (rpi/4.0 - rad*zphiu/2.0)1075 zynpu = 0. - 2.0 * np.sin ( rad*zlamu) * np.tan (rpi/4.0 - rad*zphiu/2.0)1669 zxnpu = 0. - 2.0 * np.cos (RAD*zlamu) * np.tan (RPI/4.0 - RAD*zphiu/2.0) 1670 zynpu = 0. - 2.0 * np.sin (RAD*zlamu) * np.tan (RPI/4.0 - RAD*zphiu/2.0) 1076 1671 znnpu = zxnpu*zxnpu + zynpu*zynpu 1077 1672 1078 1673 # north pole direction & modulous (at V-point) 1079 zxnpv = 0. - 2.0 * np.cos ( rad*zlamv) * np.tan (rpi/4.0 - rad*zphiv/2.0)1080 zynpv = 0. - 2.0 * np.sin ( rad*zlamv) * np.tan (rpi/4.0 - rad*zphiv/2.0)1674 zxnpv = 0. - 2.0 * np.cos (RAD*zlamv) * np.tan (RPI/4.0 - RAD*zphiv/2.0) 1675 zynpv = 0. - 2.0 * np.sin (RAD*zlamv) * np.tan (RPI/4.0 - RAD*zphiv/2.0) 1081 1676 znnpv = zxnpv*zxnpv + zynpv*zynpv 1082 1677 1083 1678 # north pole direction & modulous (at F-point) 1084 zxnpf = 0. - 2.0 * np.cos( rad*zlamf ) * np.tan ( rpi/4. - rad*zphif/2. )1085 zynpf = 0. - 2.0 * np.sin( rad*zlamf ) * np.tan ( rpi/4. - rad*zphif/2. )1679 zxnpf = 0. - 2.0 * np.cos( RAD*zlamf ) * np.tan ( RPI/4. - RAD*zphif/2. ) 1680 zynpf = 0. - 2.0 * np.sin( RAD*zlamf ) * np.tan ( RPI/4. - RAD*zphif/2. ) 1086 1681 znnpf = zxnpf*zxnpf + zynpf*zynpf 1087 1682 1088 1683 # j-direction: v-point segment direction (around T-point) 1089 zlam = zlamv 1684 zlam = zlamv 1090 1685 zphi = zphiv 1091 1686 zlan = np.roll ( zlamv, axis=-2, shift=1) # glamv (ji,jj-1) 1092 1687 zphh = np.roll ( zphiv, axis=-2, shift=1) # gphiv (ji,jj-1) 1093 zxvvt = 2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1094 - 2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1095 zyvvt = 2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1096 - 2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1688 zxvvt = 2.0 * np.cos ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1689 - 2.0 * np.cos ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1690 zyvvt = 2.0 * np.sin ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1691 - 2.0 * np.sin ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1097 1692 znvvt = np.sqrt ( znnpt * ( zxvvt*zxvvt + zyvvt*zyvvt ) ) 1098 1693 … … 1102 1697 zlan = np.roll (zlamf, axis=-2, shift=1) # glamf (ji,jj-1) 1103 1698 zphh = np.roll (zphif, axis=-2, shift=1) # gphif (ji,jj-1) 1104 zxffu = 2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1105 - 2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1106 zyffu = 2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1107 - 2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1699 zxffu = 2.0 * np.cos ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1700 - 2.0 * np.cos ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1701 zyffu = 2.0 * np.sin ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1702 - 2.0 * np.sin ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1108 1703 znffu = np.sqrt ( znnpu * ( zxffu*zxffu + zyffu*zyffu ) ) 1109 1704 1110 1705 # i-direction: f-point segment direction (around v-point) 1111 zlam = zlamf 1706 zlam = zlamf 1112 1707 zphi = zphif 1113 1708 zlan = np.roll (zlamf, axis=-1, shift=1) # glamf (ji-1,jj) 1114 1709 zphh = np.roll (zphif, axis=-1, shift=1) # gphif (ji-1,jj) 1115 zxffv = 2.0 * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1116 - 2.0 * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1117 zyffv = 2.0 * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1118 - 2.0 * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1710 zxffv = 2.0 * np.cos ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1711 - 2.0 * np.cos ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1712 zyffv = 2.0 * np.sin ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1713 - 2.0 * np.sin ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1119 1714 znffv = np.sqrt ( znnpv * ( zxffv*zxffv + zyffv*zyffv ) ) 1120 1715 1121 1716 # j-direction: u-point segment direction (around f-point) 1122 zlam = np.roll (zlamu, axis=-2, shift=-1) # glamu (ji,jj+1) 1717 zlam = np.roll (zlamu, axis=-2, shift=-1) # glamu (ji,jj+1) 1123 1718 zphi = np.roll (zphiu, axis=-2, shift=-1) # gphiu (ji,jj+1) 1124 1719 zlan = zlamu 1125 1720 zphh = zphiu 1126 zxuuf = 2. * np.cos ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1127 - 2. * np.cos ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1128 zyuuf = 2. * np.sin ( rad*zlam ) * np.tan ( rpi/4. - rad*zphi/2. )\1129 - 2. * np.sin ( rad*zlan ) * np.tan ( rpi/4. - rad*zphh/2. )1721 zxuuf = 2. * np.cos ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1722 - 2. * np.cos ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1723 zyuuf = 2. * np.sin ( RAD*zlam ) * np.tan ( RPI/4. - RAD*zphi/2. ) \ 1724 - 2. * np.sin ( RAD*zlan ) * np.tan ( RPI/4. - RAD*zphh/2. ) 1130 1725 znuuf = np.sqrt ( znnpf * ( zxuuf*zxuuf + zyuuf*zyuuf ) ) 1131 1726 1132 1727 1133 1728 # cosinus and sinus using scalar and vectorial products 1134 1729 gsint = ( zxnpt*zyvvt - zynpt*zxvvt ) / znvvt 1135 1730 gcost = ( zxnpt*zxvvt + zynpt*zyvvt ) / znvvt 1136 1731 1137 1732 gsinu = ( zxnpu*zyffu - zynpu*zxffu ) / znffu 1138 1733 gcosu = ( zxnpu*zxffu + zynpu*zyffu ) / znffu 1139 1734 1140 1735 gsinf = ( zxnpf*zyuuf - zynpf*zxuuf ) / znuuf 1141 1736 gcosf = ( zxnpf*zxuuf + zynpf*zyuuf ) / znuuf 1142 1737 1143 1738 gsinv = ( zxnpv*zxffv + zynpv*zyffv ) / znffv 1144 gcosv =-( zxnpv*zyffv - zynpv*zxffv ) / znffv # (caution, rotation of 90 degres) 1145 1146 #gsint = lbc (gsint, cd_type='T', nperio=nperio, psgn=-1.) 1147 #gcost = lbc (gcost, cd_type='T', nperio=nperio, psgn=-1.) 1148 #gsinu = lbc (gsinu, cd_type='U', nperio=nperio, psgn=-1.) 1149 #gcosu = lbc (gcosu, cd_type='U', nperio=nperio, psgn=-1.) 1150 #gsinv = lbc (gsinv, cd_type='V', nperio=nperio, psgn=-1.) 1151 #gcosv = lbc (gcosv, cd_type='V', nperio=nperio, psgn=-1.) 1152 #gsinf = lbc (gsinf, cd_type='F', nperio=nperio, psgn=-1.) 1153 #gcosf = lbc (gcosf, cd_type='F', nperio=nperio, psgn=-1.) 1739 # (caution, rotation of 90 degres) 1740 gcosv =-( zxnpv*zyffv - zynpv*zxffv ) / znffv 1154 1741 1155 1742 gsint = lbc_del (gsint, cd_type='T', nperio=nperio, psgn=-1.) … … 1175 1762 1176 1763 def angle (glam, gphi, nperio, cd_type='T') : 1177 '''Compute sinus and cosinus of model line direction with respect to east''' 1764 '''Computes sinus and cosinus of model line direction with 1765 respect to east 1766 ''' 1178 1767 mmath = __mmath__ (glam) 1179 1768 1180 1769 zlam = lbc_add (glam, nperio, cd_type, 1.) 1181 1770 zphi = lbc_add (gphi, nperio, cd_type, 1.) 1182 1771 1183 1772 # north pole direction & modulous 1184 zxnp = 0. - 2.0 * np.cos ( rad*zlam) * np.tan (rpi/4.0 - rad*zphi/2.0)1185 zynp = 0. - 2.0 * np.sin ( rad*zlam) * np.tan (rpi/4.0 - rad*zphi/2.0)1773 zxnp = 0. - 2.0 * np.cos (RAD*zlam) * np.tan (RPI/4.0 - RAD*zphi/2.0) 1774 zynp = 0. - 2.0 * np.sin (RAD*zlam) * np.tan (RPI/4.0 - RAD*zphi/2.0) 1186 1775 znnp = zxnp*zxnp + zynp*zynp 1187 1776 … … 1191 1780 zlan_s = np.roll (zlam, axis=-2, shift= 1) # glam [jj-1, ji] 1192 1781 zphh_s = np.roll (zphi, axis=-2, shift= 1) # gphi [jj-1, ji] 1193 1194 zxff = 2.0 * np.cos ( rad*zlan_n) * np.tan (rpi/4.0 - rad*zphh_n/2.0) \1195 - 2.0 * np.cos ( rad*zlan_s) * np.tan (rpi/4.0 - rad*zphh_s/2.0)1196 zyff = 2.0 * np.sin ( rad*zlan_n) * np.tan (rpi/4.0 - rad*zphh_n/2.0) \1197 - 2.0 * np.sin ( rad*zlan_s) * np.tan (rpi/4.0 - rad*zphh_s/2.0)1782 1783 zxff = 2.0 * np.cos (RAD*zlan_n) * np.tan (RPI/4.0 - RAD*zphh_n/2.0) \ 1784 - 2.0 * np.cos (RAD*zlan_s) * np.tan (RPI/4.0 - RAD*zphh_s/2.0) 1785 zyff = 2.0 * np.sin (RAD*zlan_n) * np.tan (RPI/4.0 - RAD*zphh_n/2.0) \ 1786 - 2.0 * np.sin (RAD*zlan_s) * np.tan (RPI/4.0 - RAD*zphh_s/2.0) 1198 1787 znff = np.sqrt (znnp * (zxff*zxff + zyff*zyff) ) 1199 1788 1200 1789 gsin = (zxnp*zyff - zynp*zxff) / znff 1201 1790 gcos = (zxnp*zxff + zynp*zyff) / znff … … 1207 1796 gsin = gsin.assign_coords ( glam.coords ) 1208 1797 gcos = gcos.assign_coords ( glam.coords ) 1209 1798 1210 1799 return gsin, gcos 1211 1800 1212 1801 def rot_en2ij ( u_e, v_n, gsin, gcos, nperio, cd_type ) : 1213 ''' 1214 ** Purpose : Rotate the Repere: Change vector componantes between 1802 '''Rotates the Repere: Change vector componantes between 1215 1803 geographic grid --> stretched coordinates grid. 1216 All components are on the same grid (T, U, V or F) 1804 1805 All components are on the same grid (T, U, V or F) 1217 1806 ''' 1218 1807 1219 1808 u_i = + u_e * gcos + v_n * gsin 1220 1809 v_j = - u_e * gsin + v_n * gcos 1221 1810 1222 1811 u_i = lbc (u_i, nperio=nperio, cd_type=cd_type, psgn=-1.0) 1223 1812 v_j = lbc (v_j, nperio=nperio, cd_type=cd_type, psgn=-1.0) 1224 1813 1225 1814 return u_i, v_j 1226 1815 1227 1816 def rot_ij2en ( u_i, v_j, gsin, gcos, nperio, cd_type='T' ) : 1228 ''' 1229 ** Purpose : Rotate the Repere: Change vector componantes from 1817 '''Rotates the Repere: Change vector componantes from 1230 1818 stretched coordinates grid --> geographic grid 1231 All components are on the same grid (T, U, V or F) 1819 1820 All components are on the same grid (T, U, V or F) 1232 1821 ''' 1233 1822 u_e = + u_i * gcos - v_j * gsin 1234 1823 v_n = + u_i * gsin + v_j * gcos 1824 1825 u_e = lbc (u_e, nperio=nperio, cd_type=cd_type, psgn=1.0) 1826 v_n = lbc (v_n, nperio=nperio, cd_type=cd_type, psgn=1.0) 1827 1828 return u_e, v_n 1829 1830 def rot_uv2en ( uo, vo, gsint, gcost, nperio, zdim=None ) : 1831 '''Rotate the Repere: Change vector componantes from 1832 stretched coordinates grid --> geographic grid 1833 1834 uo : velocity along i at the U grid point 1835 vo : valocity along j at the V grid point 1235 1836 1236 u_e = lbc (u_e, nperio=nperio, cd_type=cd_type, psgn= 1.0) 1237 v_n = lbc (v_n, nperio=nperio, cd_type=cd_type, psgn= 1.0) 1238 1239 return u_e, v_n 1240 1241 def rot_uv2en ( uo, vo, gsint, gcost, nperio, zdim='deptht' ) : 1242 ''' 1243 ** Purpose : Rotate the Repere: Change vector componantes from 1244 stretched coordinates grid --> geographic grid 1245 uo is on the U grid point, vo is on the V grid point 1246 east-north components on the T grid point 1247 ''' 1248 mmath = __mmath__ (uo) 1249 1250 ut = U2T (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 1251 vt = V2T (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 1252 1837 Returns east-north components on the T grid point 1838 ''' 1839 ut = u2t (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 1840 vt = v2t (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 1841 1253 1842 u_e = + ut * gcost - vt * gsint 1254 1843 v_n = + ut * gsint + vt * gcost … … 1256 1845 u_e = lbc (u_e, nperio=nperio, cd_type='T', psgn=1.0) 1257 1846 v_n = lbc (v_n, nperio=nperio, cd_type='T', psgn=1.0) 1847 1848 return u_e, v_n 1849 1850 def rot_uv2enf ( uo, vo, gsinf, gcosf, nperio, zdim=None ) : 1851 '''Rotates the Repere: Change vector componantes from 1852 stretched coordinates grid --> geographic grid 1853 1854 uo : velocity along i at the U grid point 1855 vo : valocity along j at the V grid point 1258 1856 1259 return u_e, v_n 1260 1261 def rot_uv2enF ( uo, vo, gsinf, gcosf, nperio, zdim='deptht' ) : 1262 ''' 1263 ** Purpose : Rotate the Repere: Change vector componantes from 1264 stretched coordinates grid --> geographic grid 1265 uo is on the U grid point, vo is on the V grid point 1266 east-north components on the T grid point 1267 ''' 1268 mmath = __mmath__ (uo) 1269 1270 uf = U2F (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 1271 vf = V2F (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 1272 1857 Returns east-north components on the F grid point 1858 ''' 1859 uf = u2f (uo, nperio=nperio, psgn=-1.0, zdim=zdim) 1860 vf = v2f (vo, nperio=nperio, psgn=-1.0, zdim=zdim) 1861 1273 1862 u_e = + uf * gcosf - vf * gsinf 1274 1863 v_n = + uf * gsinf + vf * gcosf … … 1276 1865 u_e = lbc (u_e, nperio=nperio, cd_type='F', psgn= 1.0) 1277 1866 v_n = lbc (v_n, nperio=nperio, cd_type='F', psgn= 1.0) 1278 1867 1279 1868 return u_e, v_n 1280 1869 1281 #@numba.jit(forceobj=True) 1282 def U2T (utab, nperio=None, psgn=-1.0, zdim='deptht', action='ave') : 1283 ''' Interpolate an array from U grid to T grid i-mean)'''1870 def u2t (utab, nperio=None, psgn=-1.0, zdim=None, action='ave') : 1871 '''Interpolates an array from U grid to T grid (i-mean) 1872 ''' 1284 1873 mmath = __mmath__ (utab) 1285 1874 utab_0 = mmath.where ( np.isnan(utab), 0., utab) 1286 lperio, aperio = lbc_diag (nperio)1875 #lperio, aperio = lbc_diag (nperio) 1287 1876 utab_0 = lbc_add (utab_0, nperio=nperio, cd_type='U', psgn=psgn) 1288 ix, ax = __findAxis__ (utab_0, 'x') 1289 iz, az = __findAxis__ (utab_0, 'z') 1290 if action == 'ave' : ttab = 0.5 * (utab_0 + np.roll (utab_0, axis=ix, shift=1)) 1291 if action == 'min' : ttab = np.minimum (utab_0 , np.roll (utab_0, axis=ix, shift=1)) 1292 if action == 'max' : ttab = np.maximum (utab_0 , np.roll (utab_0, axis=ix, shift=1)) 1293 if action == 'mult' : ttab = utab_0 * np.roll (utab_0, axis=ix, shift=1) 1294 ttab = lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn) 1295 1877 ax, ix = __find_axis__ (utab_0, 'x') 1878 az = __find_axis__ (utab_0, 'z')[0] 1879 1880 if ax : 1881 if action == 'ave' : 1882 ttab = 0.5 * (utab_0 + np.roll (utab_0, axis=ix, shift=1)) 1883 if action == 'min' : 1884 ttab = np.minimum (utab_0 , np.roll (utab_0, axis=ix, shift=1)) 1885 if action == 'max' : 1886 ttab = np.maximum (utab_0 , np.roll (utab_0, axis=ix, shift=1)) 1887 if action == 'mult': 1888 ttab = utab_0 * np.roll (utab_0, axis=ix, shift=1) 1889 ttab = lbc_del (ttab , nperio=nperio, cd_type='T', psgn=psgn) 1890 else : 1891 ttab = lbc_del (utab_0, nperio=nperio, cd_type='T', psgn=psgn) 1892 1296 1893 if mmath == xr : 1297 if ax != None:1894 if ax : 1298 1895 ttab = ttab.assign_coords({ax:np.arange (ttab.shape[ix])+1.}) 1299 if zdim != None and iz != None and az != 'olevel' : 1300 ttab = ttab.rename( {az:zdim}) 1896 if zdim and az : 1897 if az != zdim : 1898 ttab = ttab.rename( {az:zdim}) 1301 1899 return ttab 1302 1900 1303 #@numba.jit(forceobj=True) 1304 def V2T (vtab, nperio=None, psgn=-1.0, zdim='deptht', action='ave') : 1305 ''' Interpolate an array from V grid to T grid (j-mean)'''1901 def v2t (vtab, nperio=None, psgn=-1.0, zdim=None, action='ave') : 1902 '''Interpolates an array from V grid to T grid (j-mean) 1903 ''' 1306 1904 mmath = __mmath__ (vtab) 1307 lperio, aperio = lbc_diag (nperio)1905 #lperio, aperio = lbc_diag (nperio) 1308 1906 vtab_0 = mmath.where ( np.isnan(vtab), 0., vtab) 1309 1907 vtab_0 = lbc_add (vtab_0, nperio=nperio, cd_type='V', psgn=psgn) 1310 iy, ay = __findAxis__ (vtab_0, 'y') 1311 iz, az = __findAxis__ (vtab_0, 'z') 1312 if action == 'ave' : ttab = 0.5 * (vtab_0 + np.roll (vtab_0, axis=iy, shift=1)) 1313 if action == 'min' : ttab = np.minimum (vtab_0 , np.roll (vtab_0, axis=iy, shift=1)) 1314 if action == 'max' : ttab = np.maximum (vtab_0 , np.roll (vtab_0, axis=iy, shift=1)) 1315 if action == 'mult' : ttab = vtab_0 * np.roll (vtab_0, axis=iy, shift=1) 1316 ttab = lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn) 1908 ay, jy = __find_axis__ (vtab_0, 'y') 1909 az = __find_axis__ (vtab_0, 'z')[0] 1910 if ay : 1911 if action == 'ave' : 1912 ttab = 0.5 * (vtab_0 + np.roll (vtab_0, axis=jy, shift=1)) 1913 if action == 'min' : 1914 ttab = np.minimum (vtab_0 , np.roll (vtab_0, axis=jy, shift=1)) 1915 if action == 'max' : 1916 ttab = np.maximum (vtab_0 , np.roll (vtab_0, axis=jy, shift=1)) 1917 if action == 'mult' : 1918 ttab = vtab_0 * np.roll (vtab_0, axis=jy, shift=1) 1919 ttab = lbc_del (ttab , nperio=nperio, cd_type='T', psgn=psgn) 1920 else : 1921 ttab = lbc_del (vtab_0, nperio=nperio, cd_type='T', psgn=psgn) 1922 1317 1923 if mmath == xr : 1318 if ay !=None : 1319 ttab = ttab.assign_coords({ay:np.arange(ttab.shape[iy])+1.}) 1320 if zdim != None and iz != None and az != 'olevel' : 1321 ttab = ttab.rename( {az:zdim}) 1924 if ay : 1925 ttab = ttab.assign_coords({ay:np.arange(ttab.shape[jy])+1.}) 1926 if zdim and az : 1927 if az != zdim : 1928 ttab = ttab.rename( {az:zdim}) 1322 1929 return ttab 1323 1930 1324 #@numba.jit(forceobj=True) 1325 def F2T (ftab, nperio=None, psgn=1.0, zdim='depthf', action='ave') : 1326 ''' Interpolate an array from F grid to T grid (i- and j- means)'''1931 def f2t (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') : 1932 '''Interpolates an array from F grid to T grid (i- and j- means) 1933 ''' 1327 1934 mmath = __mmath__ (ftab) 1328 1935 ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab) 1329 1936 ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 1330 ttab = V2T(F2V(ftab_0, nperio=nperio, psgn=psgn, zdim=zdim, action=action), nperio=nperio, psgn=psgn, zdim=zdim, action=action) 1937 ttab = v2t (f2v (ftab_0, nperio=nperio, psgn=psgn, zdim=zdim, action=action), 1938 nperio=nperio, psgn=psgn, zdim=zdim, action=action) 1331 1939 return lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn) 1332 1940 1333 #@numba.jit(forceobj=True) 1334 def T2U (ttab, nperio=None, psgn=1.0, zdim='depthu', action='ave') : 1335 ''' Interpolate an array from T grid to U grid (i-mean)'''1941 def t2u (ttab, nperio=None, psgn=1.0, zdim=None, action='ave') : 1942 '''Interpolates an array from T grid to U grid (i-mean) 1943 ''' 1336 1944 mmath = __mmath__ (ttab) 1337 1945 ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 1338 1946 ttab_0 = lbc_add (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn) 1339 ix, ax = __findAxis__ (ttab_0, 'x') 1340 iz, az = __findAxis__ (ttab_0, 'z') 1341 if action == 'ave' : utab = 0.5 * (ttab_0 + np.roll (ttab_0, axis=ix, shift=-1)) 1342 if action == 'min' : utab = np.minimum (ttab_0 , np.roll (ttab_0, axis=ix, shift=-1)) 1343 if action == 'max' : utab = np.maximum (ttab_0 , np.roll (ttab_0, axis=ix, shift=-1)) 1344 if action == 'mult' : utab = ttab_0 * np.roll (ttab_0, axis=ix, shift=-1) 1345 utab = lbc_del (utab, nperio=nperio, cd_type='U', psgn=psgn) 1346 1347 if mmath == xr : 1348 if ax != None : 1947 ax, ix = __find_axis__ (ttab_0, 'x')[0] 1948 az = __find_axis__ (ttab_0, 'z') 1949 if ix : 1950 if action == 'ave' : 1951 utab = 0.5 * (ttab_0 + np.roll (ttab_0, axis=ix, shift=-1)) 1952 if action == 'min' : 1953 utab = np.minimum (ttab_0 , np.roll (ttab_0, axis=ix, shift=-1)) 1954 if action == 'max' : 1955 utab = np.maximum (ttab_0 , np.roll (ttab_0, axis=ix, shift=-1)) 1956 if action == 'mult' : 1957 utab = ttab_0 * np.roll (ttab_0, axis=ix, shift=-1) 1958 utab = lbc_del (utab , nperio=nperio, cd_type='U', psgn=psgn) 1959 else : 1960 utab = lbc_del (ttab_0, nperio=nperio, cd_type='U', psgn=psgn) 1961 1962 if mmath == xr : 1963 if ax : 1349 1964 utab = ttab.assign_coords({ax:np.arange(utab.shape[ix])+1.}) 1350 if zdim != None and iz != None and az != 'olevel' : 1351 utab = utab.rename( {az:zdim}) 1965 if zdim and az : 1966 if az != zdim : 1967 utab = utab.rename( {az:zdim}) 1352 1968 return utab 1353 1969 1354 #@numba.jit(forceobj=True) 1355 def T2V (ttab, nperio=None, psgn=1.0, zdim='depthv', action='ave') : 1356 ''' Interpolate an array from T grid to V grid (j-mean)'''1970 def t2v (ttab, nperio=None, psgn=1.0, zdim=None, action='ave') : 1971 '''Interpolates an array from T grid to V grid (j-mean) 1972 ''' 1357 1973 mmath = __mmath__ (ttab) 1358 1974 ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 1359 1975 ttab_0 = lbc_add (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn) 1360 iy, ay = __findAxis__ (ttab_0, 'y') 1361 iz, az = __findAxis__ (ttab_0, 'z') 1362 if action == 'ave' : vtab = 0.5 * (ttab_0 + np.roll (ttab_0, axis=iy, shift=-1)) 1363 if action == 'min' : vtab = np.minimum (ttab_0 , np.roll (ttab_0, axis=iy, shift=-1)) 1364 if action == 'max' : vtab = np.maximum (ttab_0 , np.roll (ttab_0, axis=iy, shift=-1)) 1365 if action == 'mult' : vtab = ttab_0 * np.roll (ttab_0, axis=iy, shift=-1) 1366 1367 vtab = lbc_del (vtab, nperio=nperio, cd_type='V', psgn=psgn) 1976 ay, jy = __find_axis__ (ttab_0, 'y') 1977 az = __find_axis__ (ttab_0, 'z')[0] 1978 if jy : 1979 if action == 'ave' : 1980 vtab = 0.5 * (ttab_0 + np.roll (ttab_0, axis=jy, shift=-1)) 1981 if action == 'min' : 1982 vtab = np.minimum (ttab_0 , np.roll (ttab_0, axis=jy, shift=-1)) 1983 if action == 'max' : 1984 vtab = np.maximum (ttab_0 , np.roll (ttab_0, axis=jy, shift=-1)) 1985 if action == 'mult' : 1986 vtab = ttab_0 * np.roll (ttab_0, axis=jy, shift=-1) 1987 vtab = lbc_del (vtab , nperio=nperio, cd_type='V', psgn=psgn) 1988 else : 1989 vtab = lbc_del (ttab_0, nperio=nperio, cd_type='V', psgn=psgn) 1990 1368 1991 if mmath == xr : 1369 if ay != None : 1370 vtab = vtab.assign_coords({ay:np.arange(vtab.shape[iy])+1.}) 1371 if zdim != None and iz != None and az != 'olevel' : 1372 vtab = vtab.rename( {az:zdim}) 1992 if ay : 1993 vtab = vtab.assign_coords({ay:np.arange(vtab.shape[jy])+1.}) 1994 if zdim and az : 1995 if az != zdim : 1996 vtab = vtab.rename( {az:zdim}) 1373 1997 return vtab 1374 1998 1375 #@numba.jit(forceobj=True) 1376 def V2F (vtab, nperio=None, psgn=-1.0, zdim='depthf', action='ave') : 1377 ''' Interpolate an array from V grid to F grid (i-mean)'''1999 def v2f (vtab, nperio=None, psgn=-1.0, zdim=None, action='ave') : 2000 '''Interpolates an array from V grid to F grid (i-mean) 2001 ''' 1378 2002 mmath = __mmath__ (vtab) 1379 2003 vtab_0 = mmath.where ( np.isnan(vtab), 0., vtab) 1380 2004 vtab_0 = lbc_add (vtab_0 , nperio=nperio, cd_type='V', psgn=psgn) 1381 ix, ax = __findAxis__ (vtab_0, 'x') 1382 iz, az = __findAxis__ (vtab_0, 'z') 1383 if action == 'ave' : 0.5 * (vtab_0 + np.roll (vtab_0, axis=ix, shift=-1)) 1384 if action == 'min' : np.minimum (vtab_0 , np.roll (vtab_0, axis=ix, shift=-1)) 1385 if action == 'max' : np.maximum (vtab_0 , np.roll (vtab_0, axis=ix, shift=-1)) 1386 if action == 'mult' : vtab_0 * np.roll (vtab_0, axis=ix, shift=-1) 1387 ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 1388 2005 ax, ix = __find_axis__ (vtab_0, 'x') 2006 az = __find_axis__ (vtab_0, 'z')[0] 2007 if ix : 2008 if action == 'ave' : 2009 ftab = 0.5 * (vtab_0 + np.roll (vtab_0, axis=ix, shift=-1)) 2010 if action == 'min' : 2011 ftab = np.minimum (vtab_0 , np.roll (vtab_0, axis=ix, shift=-1)) 2012 if action == 'max' : 2013 ftab = np.maximum (vtab_0 , np.roll (vtab_0, axis=ix, shift=-1)) 2014 if action == 'mult' : 2015 ftab = vtab_0 * np.roll (vtab_0, axis=ix, shift=-1) 2016 ftab = lbc_del (ftab , nperio=nperio, cd_type='F', psgn=psgn) 2017 else : 2018 ftab = lbc_del (vtab_0, nperio=nperio, cd_type='F', psgn=psgn) 2019 1389 2020 if mmath == xr : 1390 if ax != None :2021 if ax : 1391 2022 ftab = ftab.assign_coords({ax:np.arange(ftab.shape[ix])+1.}) 1392 if zdim != None and iz != None and az != 'olevel' : 1393 ftab = ftab.rename( {az:zdim}) 2023 if zdim and az : 2024 if az != zdim : 2025 ftab = ftab.rename( {az:zdim}) 1394 2026 return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 1395 2027 1396 #@numba.jit(forceobj=True) 1397 def U2F (utab, nperio=None, psgn=-1.0, zdim='depthf', action='ave') : 1398 ''' Interpolate an array from U grid to F grid i-mean)'''2028 def u2f (utab, nperio=None, psgn=-1.0, zdim=None, action='ave') : 2029 '''Interpolates an array from U grid to F grid i-mean) 2030 ''' 1399 2031 mmath = __mmath__ (utab) 1400 2032 utab_0 = mmath.where ( np.isnan(utab), 0., utab) 1401 2033 utab_0 = lbc_add (utab_0 , nperio=nperio, cd_type='U', psgn=psgn) 1402 iy, ay = __findAxis__ (utab_0, 'y') 1403 iz, az = __findAxis__ (utab_0, 'z') 1404 if action == 'ave' : ftab = 0.5 * (utab_0 + np.roll (utab_0, axis=iy, shift=-1)) 1405 if action == 'min' : ftab = np.minimum (utab_0 , np.roll (utab_0, axis=iy, shift=-1)) 1406 if action == 'max' : ftab = np.maximum (utab_0 , np.roll (utab_0, axis=iy, shift=-1)) 1407 if action == 'mult' : ftab = utab_0 * np.roll (utab_0, axis=iy, shift=-1) 1408 ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 2034 ay, jy = __find_axis__ (utab_0, 'y') 2035 az = __find_axis__ (utab_0, 'z')[0] 2036 if jy : 2037 if action == 'ave' : 2038 ftab = 0.5 * (utab_0 + np.roll (utab_0, axis=jy, shift=-1)) 2039 if action == 'min' : 2040 ftab = np.minimum (utab_0 , np.roll (utab_0, axis=jy, shift=-1)) 2041 if action == 'max' : 2042 ftab = np.maximum (utab_0 , np.roll (utab_0, axis=jy, shift=-1)) 2043 if action == 'mult' : 2044 ftab = utab_0 * np.roll (utab_0, axis=jy, shift=-1) 2045 ftab = lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 2046 else : 2047 ftab = lbc_del (utab_0, nperio=nperio, cd_type='F', psgn=psgn) 1409 2048 1410 2049 if mmath == xr : 1411 if ay != None : 1412 ftab = ftab.assign_coords({'y':np.arange(ftab.shape[iy])+1.}) 1413 if zdim != None and iz != None and az != 'olevel' : 1414 ftab = ftab.rename( {az:zdim}) 2050 if ay : 2051 ftab = ftab.assign_coords({'y':np.arange(ftab.shape[jy])+1.}) 2052 if zdim and az : 2053 if az != zdim : 2054 ftab = ftab.rename( {az:zdim}) 1415 2055 return ftab 1416 2056 1417 #@numba.jit(forceobj=True) 1418 def F2T (ftab, nperio=None, psgn=1.0, zdim='deptht', action='ave') : 1419 '''Interpolate an array on F grid to T grid (i- and j- means)''' 1420 mmath = __mmath__ (ftab) 1421 ftab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 1422 ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 1423 ttab = U2T(F2U(ftab_0, nperio=nperio, psgn=psgn, zdim=zdim, action=action), nperio=nperio, psgn=psgn, zdim=zdim, action=action) 1424 return lbc_del (ttab, nperio=nperio, cd_type='T', psgn=psgn) 1425 1426 #@numba.jit(forceobj=True) 1427 def T2F (ttab, nperio=None, psgn=1.0, zdim='deptht', action='mean') : 1428 '''Interpolate an array on T grid to F grid (i- and j- means)''' 2057 def t2f (ttab, nperio=None, psgn=1.0, zdim=None, action='mean') : 2058 '''Interpolates an array on T grid to F grid (i- and j- means) 2059 ''' 1429 2060 mmath = __mmath__ (ttab) 1430 2061 ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 1431 2062 ttab_0 = lbc_add (ttab_0 , nperio=nperio, cd_type='T', psgn=psgn) 1432 ftab = T2U(U2F(ttab, nperio=nperio, psgn=psgn, zdim=zdim, action=action), nperio=nperio, psgn=psgn, zdim=zdim, action=action) 1433 2063 ftab = t2u (u2f (ttab, nperio=nperio, psgn=psgn, zdim=zdim, action=action), 2064 nperio=nperio, psgn=psgn, zdim=zdim, action=action) 2065 1434 2066 return lbc_del (ftab, nperio=nperio, cd_type='F', psgn=psgn) 1435 2067 1436 #@numba.jit(forceobj=True) 1437 def F2U (ftab, nperio=None, psgn=1.0, zdim='depthu', action='ave') : 1438 ''' Interpolate an array on F grid to FUgrid (i-mean)'''2068 def f2u (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') : 2069 '''Interpolates an array on F grid to U grid (j-mean) 2070 ''' 1439 2071 mmath = __mmath__ (ftab) 1440 2072 ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab) 1441 2073 ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 1442 iy, ay = __findAxis__ (ftab_0, 'y') 1443 iz, az = __findAxis__ (ftab_0, 'z') 1444 if action == 'ave' : utab = 0.5 * (ftab_0 + np.roll (ftab_0, axis=iy, shift=-1)) 1445 if action == 'min' : utab = np.minimum (ftab_0 , np.roll (ftab_0, axis=iy, shift=-1)) 1446 if action == 'max' : utab = np.maximum (ftab_0 , np.roll (ftab_0, axis=iy, shift=-1)) 1447 if action == 'mult' : utab = ftab_0 * np.roll (ftab_0, axis=iy, shift=-1) 1448 1449 utab = lbc_del (utab, nperio=nperio, cd_type='U', psgn=psgn) 1450 2074 ay, jy = __find_axis__ (ftab_0, 'y') 2075 az = __find_axis__ (ftab_0, 'z')[0] 2076 if jy : 2077 if action == 'ave' : 2078 utab = 0.5 * (ftab_0 + np.roll (ftab_0, axis=jy, shift=-1)) 2079 if action == 'min' : 2080 utab = np.minimum (ftab_0 , np.roll (ftab_0, axis=jy, shift=-1)) 2081 if action == 'max' : 2082 utab = np.maximum (ftab_0 , np.roll (ftab_0, axis=jy, shift=-1)) 2083 if action == 'mult' : 2084 utab = ftab_0 * np.roll (ftab_0, axis=jy, shift=-1) 2085 utab = lbc_del (utab , nperio=nperio, cd_type='U', psgn=psgn) 2086 else : 2087 utab = lbc_del (ftab_0, nperio=nperio, cd_type='U', psgn=psgn) 2088 1451 2089 if mmath == xr : 1452 utab = utab.assign_coords({ay:np.arange(ftab.shape[ iy])+1.})1453 if zdim != None and iz != None and az != 'olevel':1454 utab = utab.rename( {az:zdim}) 2090 utab = utab.assign_coords({ay:np.arange(ftab.shape[jy])+1.}) 2091 if zdim and az and az != zdim : 2092 utab = utab.rename( {az:zdim}) 1455 2093 return utab 1456 2094 1457 #@numba.jit(forceobj=True) 1458 def F2V (ftab, nperio=None, psgn=1.0, zdim='depthv', action='ave') : 1459 ''' Interpolate an array from F grid to V grid (i-mean)'''2095 def f2v (ftab, nperio=None, psgn=1.0, zdim=None, action='ave') : 2096 '''Interpolates an array from F grid to V grid (i-mean) 2097 ''' 1460 2098 mmath = __mmath__ (ftab) 1461 2099 ftab_0 = mmath.where ( np.isnan(ftab), 0., ftab) 1462 2100 ftab_0 = lbc_add (ftab_0 , nperio=nperio, cd_type='F', psgn=psgn) 1463 ix, ax = __findAxis__ (ftab_0, 'x') 1464 iz, az = __findAxis__ (ftab_0, 'z') 1465 if action == 'ave' : vtab = 0.5 * (ftab_0 + np.roll (ftab_0, axis=ix, shift=-1)) 1466 if action == 'min' : vtab = np.minimum (ftab_0 , np.roll (ftab_0, axis=ix, shift=-1)) 1467 if action == 'max' : vtab = np.maximum (ftab_0 , np.roll (ftab_0, axis=ix, shift=-1)) 1468 if action == 'mult' : vtab = ftab_0 * np.roll (ftab_0, axis=ix, shift=-1) 1469 1470 vtab = lbc_del (vtab, nperio=nperio, cd_type='V', psgn=psgn) 2101 ax, ix = __find_axis__ (ftab_0, 'x') 2102 az = __find_axis__ (ftab_0, 'z')[0] 2103 if ix : 2104 if action == 'ave' : 2105 vtab = 0.5 * (ftab_0 + np.roll (ftab_0, axis=ix, shift=-1)) 2106 if action == 'min' : 2107 vtab = np.minimum (ftab_0 , np.roll (ftab_0, axis=ix, shift=-1)) 2108 if action == 'max' : 2109 vtab = np.maximum (ftab_0 , np.roll (ftab_0, axis=ix, shift=-1)) 2110 if action == 'mult' : 2111 vtab = ftab_0 * np.roll (ftab_0, axis=ix, shift=-1) 2112 vtab = lbc_del (vtab , nperio=nperio, cd_type='V', psgn=psgn) 2113 else : 2114 vtab = lbc_del (ftab_0, nperio=nperio, cd_type='V', psgn=psgn) 2115 1471 2116 if mmath == xr : 1472 2117 vtab = vtab.assign_coords({ax:np.arange(ftab.shape[ix])+1.}) 1473 if zdim != None and iz != None and az != 'olevel' : 1474 vtab = vtab.rename( {az:zdim}) 2118 if zdim and az : 2119 if az != zdim : 2120 vtab = vtab.rename( {az:zdim}) 1475 2121 return vtab 1476 2122 1477 #@numba.jit(forceobj=True) 1478 def W2T (wtab, zcoord=None, zdim='deptht', sval=np.nan) : 1479 ''' 1480 Interpolate an array on W grid to T grid (k-mean) 2123 def w2t (wtab, zcoord=None, zdim=None, sval=np.nan) : 2124 '''Interpolates an array on W grid to T grid (k-mean) 2125 1481 2126 sval is the bottom value 1482 2127 ''' … … 1484 2129 wtab_0 = mmath.where ( np.isnan(wtab), 0., wtab) 1485 2130 1486 iz, az = __findAxis__ (wtab_0, 'z') 1487 1488 ttab = 0.5 * ( wtab_0 + np.roll (wtab_0, axis=iz, shift=-1) ) 1489 2131 az, kz = __find_axis__ (wtab_0, 'z') 2132 2133 if kz : 2134 ttab = 0.5 * ( wtab_0 + np.roll (wtab_0, axis=kz, shift=-1) ) 2135 else : 2136 ttab = wtab_0 2137 1490 2138 if mmath == xr : 1491 ttab[{az:iz}] = sval 1492 if zdim != None and iz != None and az != 'olevel' : 1493 ttab = ttab.rename ( {az:zdim} ) 1494 try : ttab = ttab.assign_coords ( {zdim:zcoord} ) 1495 except : pass 2139 ttab[{az:kz}] = sval 2140 if zdim and az : 2141 if az != zdim : 2142 ttab = ttab.rename ( {az:zdim} ) 2143 if zcoord is not None : 2144 ttab = ttab.assign_coords ( {zdim:zcoord} ) 1496 2145 else : 1497 2146 ttab[..., -1, :, :] = sval … … 1499 2148 return ttab 1500 2149 1501 #@numba.jit(forceobj=True) 1502 def T2W (ttab, zcoord=None, zdim='depthw', sval=np.nan, extrap_surf=False) : 1503 '''Interpolate an array from T grid to W grid (k-mean) 2150 def t2w (ttab, zcoord=None, zdim=None, sval=np.nan, extrap_surf=False) : 2151 '''Interpolates an array from T grid to W grid (k-mean) 2152 1504 2153 sval is the surface value 1505 2154 if extrap_surf==True, surface value is taken from 1st level value. … … 1507 2156 mmath = __mmath__ (ttab) 1508 2157 ttab_0 = mmath.where ( np.isnan(ttab), 0., ttab) 1509 iz, az = __findAxis__ (ttab_0, 'z')1510 wtab = 0.5 * ( ttab_0 + np.roll (ttab_0, axis= iz, shift=1) )2158 az, kz = __find_axis__ (ttab_0, 'z') 2159 wtab = 0.5 * ( ttab_0 + np.roll (ttab_0, axis=kz, shift=1) ) 1511 2160 1512 2161 if mmath == xr : 1513 if extrap_surf : wtab[{az:0}] = ttabb[{az:0}] 1514 else : wtab[{az:0}] = sval 1515 else : 1516 if extrap_surf : wtab[..., 0, :, :] = ttab[..., 0, :, :] 1517 else : wtab[..., 0, :, :] = sval 2162 if extrap_surf : 2163 wtab[{az:0}] = ttab[{az:0}] 2164 else : 2165 wtab[{az:0}] = sval 2166 else : 2167 if extrap_surf : 2168 wtab[..., 0, :, :] = ttab[..., 0, :, :] 2169 else : 2170 wtab[..., 0, :, :] = sval 1518 2171 1519 2172 if mmath == xr : 1520 if zdim != None and iz != None and az != 'olevel' : 1521 wtab = wtab.rename ( {az:zdim}) 1522 if zcoord != None : wtab = wtab.assign_coords ( {zdim:zcoord}) 1523 else : ztab = wtab.assign_coords ( {zdim:np.arange(ttab.shape[iz])+1.} ) 2173 if zdim and az and az != zdim : 2174 wtab = wtab.rename ( {az:zdim}) 2175 if zcoord is not None : 2176 wtab = wtab.assign_coords ( {zdim:zcoord}) 2177 else : 2178 wtab = wtab.assign_coords ( {zdim:np.arange(ttab.shape[kz])+1.} ) 1524 2179 return wtab 1525 2180 1526 #@numba.jit(forceobj=True) 1527 def fill (ptab, nperio, cd_type='T', npass=1, sval=0.) : 1528 ''' 1529 Fill sval values with mean of neighbours 1530 2181 def fill (ptab, nperio, cd_type='T', npass=1, sval=np.nan) : 2182 '''Fills np.nan values with mean of neighbours 2183 1531 2184 Inputs : 1532 2185 ptab : input field to fill 1533 2186 nperio, cd_type : periodicity characteristics 1534 ''' 2187 ''' 1535 2188 1536 2189 mmath = __mmath__ (ptab) 1537 2190 1538 DoPerio = False ; lperio = nperio 2191 do_perio = False 2192 lperio = nperio 1539 2193 if nperio == 4.2 : 1540 DoPerio = True ; lperio =42194 do_perio, lperio = True, 4 1541 2195 if nperio == 6.2 : 1542 DoPerio = True ; lperio =61543 1544 if DoPerio :1545 ztab = lbc_add (ptab, nperio=nperio , sval=sval)1546 else : 2196 do_perio, lperio = True, 6 2197 2198 if do_perio : 2199 ztab = lbc_add (ptab, nperio=nperio) 2200 else : 1547 2201 ztab = ptab 1548 1549 if np.isnan (sval) : 2202 2203 if np.isnan (sval) : 1550 2204 ztab = mmath.where (np.isnan(ztab), np.nan, ztab) 1551 2205 else : 1552 2206 ztab = mmath.where (ztab==sval , np.nan, ztab) 1553 1554 for nn in np.arange (npass) :2207 2208 for _ in np.arange (npass) : 1555 2209 zmask = mmath.where ( np.isnan(ztab), 0., 1. ) 1556 2210 ztab0 = mmath.where ( np.isnan(ztab), 0., ztab ) … … 1576 2230 zcount = lbc (zcount, nperio=lperio, cd_type=cd_type) 1577 2231 znew = lbc (znew , nperio=lperio, cd_type=cd_type) 1578 2232 1579 2233 ztab = mmath.where (np.logical_and (zmask==0., zcount>0), znew/zcount, ztab) 1580 2234 1581 2235 ztab = mmath.where (zcount==0, sval, ztab) 1582 if DoPerio : ztab = lbc_del (ztab, nperio=lperio) 2236 if do_perio : 2237 ztab = lbc_del (ztab, nperio=lperio) 1583 2238 1584 2239 return ztab 1585 2240 1586 #@numba.jit(forceobj=True)1587 2241 def correct_uv (u, v, lat) : 1588 2242 ''' 1589 Correct a Cartopy bug in Orthographic projection2243 Corrects a Cartopy bug in orthographic projection 1590 2244 1591 2245 See https://github.com/SciTools/cartopy/issues/1179 … … 1593 2247 The correction is needed with cartopy <= 0.20 1594 2248 It seems that version 0.21 will correct the bug (https://github.com/SciTools/cartopy/pull/1926) 2249 Later note : the bug is still present in Cartopy 0.22 1595 2250 1596 2251 Inputs : 1597 u, v : eastward/no thward components2252 u, v : eastward/northward components 1598 2253 lat : latitude of the point (degrees north) 1599 2254 1600 Outputs : 2255 Outputs : 1601 2256 modified eastward/nothward components to have correct polar projections in cartopy 1602 2257 ''' 1603 2258 uv = np.sqrt (u*u + v*v) # Original modulus 1604 2259 zu = u 1605 zv = v * np.cos ( rad*lat)2260 zv = v * np.cos (RAD*lat) 1606 2261 zz = np.sqrt ( zu*zu + zv*zv ) # Corrected modulus 1607 uc = zu*uv/zz ; vc = zv*uv/zz # Final corrected values 2262 uc = zu*uv/zz 2263 vc = zv*uv/zz # Final corrected values 1608 2264 return uc, vc 1609 2265 1610 def msf (v_e1v_e3v, lat1d, depthw) : 1611 ''' 1612 Computes the meridonal stream function 1613 First input is meridional_velocity*e1v*e3v 1614 ''' 1615 @numba.jit(forceobj=True) 1616 def iin (tab, dim) : 1617 ''' 1618 Integrate from the bottom 1619 ''' 1620 result = tab * 0.0 1621 nlen = len(tab.coords[dim]) 1622 for jn in np.arange (nlen-2, 0, -1) : 1623 result [{dim:jn}] = result [{dim:jn+1}] - tab [{dim:jn}] 1624 result = result.where (result !=0, np.nan) 1625 return result 1626 1627 zomsf = iin ((v_e1v_e3v).sum (dim='x', keep_attrs=True)*1E-6, dim='depthv') 1628 zomsf = zomsf.assign_coords ( {'depthv':depthw.values, 'y':lat1d}) 1629 zomsf = zomsf.rename ( {'depthv':'depthw', 'y':'lat'}) 2266 def norm_uv (u, v) : 2267 '''Returns norm of a 2 components vector 2268 ''' 2269 return np.sqrt (u*u + v*v) 2270 2271 def normalize_uv (u, v) : 2272 '''Normalizes 2 components vector 2273 ''' 2274 uv = norm_uv (u, v) 2275 return u/uv, v/uv 2276 2277 def msf (vv, e1v_e3v, plat1d, depthw) : 2278 '''Computes the meridonal stream function 2279 2280 vv : meridional_velocity 2281 e1v_e3v : prodcut of scale factors e1v*e3v 2282 ''' 2283 2284 v_e1v_e3v = vv * e1v_e3v 2285 v_e1v_e3v.attrs = vv.attrs 2286 2287 ax = __find_axis__ (v_e1v_e3v, 'x')[0] 2288 az = __find_axis__ (v_e1v_e3v, 'z')[0] 2289 if az == 'olevel' : 2290 new_az = 'olevel' 2291 else : 2292 new_az = 'depthw' 2293 2294 zomsf = -v_e1v_e3v.cumsum ( dim=az, keep_attrs=True).sum (dim=ax, keep_attrs=True)*1.E-6 2295 zomsf = zomsf - zomsf.isel ( { az:-1} ) 2296 2297 ay = __find_axis__ (zomsf, 'y' )[0] 2298 zomsf = zomsf.assign_coords ( {az:depthw.values, ay:plat1d.values}) 2299 2300 zomsf = zomsf.rename ( {ay:'lat'}) 2301 if az != new_az : 2302 zomsf = zomsf.rename ( {az:new_az} ) 2303 zomsf.attrs['standard_name'] = 'Meridional stream function' 1630 2304 zomsf.attrs['long_name'] = 'Meridional stream function' 1631 1632 2305 zomsf.attrs['units'] = 'Sv' 1633 zomsf .depthw.attrs=depthw.attrs1634 zomsf.lat.attrs= lat1d.attrs1635 2306 zomsf[new_az].attrs = depthw.attrs 2307 zomsf.lat.attrs=plat1d.attrs 2308 1636 2309 return zomsf 1637 2310 1638 def bsf (u_e2u_e3u, mask, nperio=None, bsf0=None ) : 1639 ''' 1640 Computes the barotropic stream function 1641 First input is zonal_velocity*e2u*e3u 1642 bsf0 is the point with bsf=0 (ex: bsf0={'x':5, 'y':120} ) 1643 ''' 1644 @numba.jit(forceobj=True) 1645 def iin (tab, dim) : 1646 ''' 1647 Integrate from the south 1648 ''' 1649 result = tab * 0.0 1650 nlen = len(tab.coords[dim]) 1651 for jn in np.arange (3, nlen) : 1652 result [{dim:jn}] = result [{dim:jn-1}] + tab [{dim:jn}] 1653 return result 1654 1655 bsf = iin ((u_e2u_e3u).sum(dim='depthu', keep_attrs=True)*1E-6, dim='y') 1656 bsf.attrs = u_e2u_e3u.attrs 1657 if bsf0 != None : 1658 bsf = bsf - bsf.isel (bsf0) 1659 1660 bsf = bsf.where (mask !=0, np.nan) 1661 bsf.attrs['long_name'] = 'Barotropic stream function' 1662 bsf.attrs['units'] = 'Sv' 1663 bsf = lbc (bsf, nperio=nperio, cd_type='F') 1664 1665 return bsf 1666 1667 def namelist_read (ref=None, cfg=None, out='dict', flat=False, verbose=False) : 1668 ''' 1669 Read NEMO namelist(s) and return either a dictionnary or an xarray dataset 1670 1671 ref : file with reference namelist, or a f90nml.namelist.Namelist object 1672 cfg : file with config namelist, or a f90nml.namelist.Namelist object 1673 At least one namelist neaded 1674 1675 out: 2311 def bsf (uu, e2u_e3u, mask, nperio=None, bsf0=None ) : 2312 '''Computes the barotropic stream function 2313 2314 uu : zonal_velocity 2315 e2u_e3u : product of scales factor e2u*e3u 2316 bsf0 : the point with bsf=0 2317 (ex: bsf0={'x':3, 'y':120} for orca2, 2318 bsf0={'x':5, 'y':300} for eORCA1 2319 ''' 2320 u_e2u_e3u = uu * e2u_e3u 2321 u_e2u_e3u.attrs = uu.attrs 2322 2323 ay = __find_axis__ (u_e2u_e3u, 'y')[0] 2324 az = __find_axis__ (u_e2u_e3u, 'z')[0] 2325 2326 zbsf = -u_e2u_e3u.cumsum ( dim=ay, keep_attrs=True ) 2327 zbsf = zbsf.sum (dim=az, keep_attrs=True)*1.E-6 2328 2329 if bsf0 : 2330 zbsf = zbsf - zbsf.isel (bsf0) 2331 2332 zbsf = zbsf.where (mask !=0, np.nan) 2333 zbsf.attrs.update (uu.attrs) 2334 zbsf.attrs['standard_name'] = 'barotropic_stream_function' 2335 zbsf.attrs['long_name'] = 'Barotropic stream function' 2336 zbsf.attrs['units'] = 'Sv' 2337 zbsf = lbc (zbsf, nperio=nperio, cd_type='F') 2338 2339 return zbsf 2340 2341 if f90nml : 2342 def namelist_read (ref=None, cfg=None, out='dict', flat=False, verbose=False) : 2343 '''Read NEMO namelist(s) and return either a dictionnary or an xarray dataset 2344 2345 ref : file with reference namelist, or a f90nml.namelist.Namelist object 2346 cfg : file with config namelist, or a f90nml.namelist.Namelist object 2347 At least one namelist neaded 2348 2349 out: 1676 2350 'dict' to return a dictonnary 1677 2351 'xr' to return an xarray dataset 1678 flat : only for dict output. Output a flat dictionnary with all values. 1679 1680 ''' 1681 1682 if ref != None : 1683 if isinstance (ref, str) : nml_ref = f90nml.read (ref) 1684 if isinstance (ref, f90nml.namelist.Namelist) : nml_ref = ref 1685 1686 if cfg != None : 1687 if isinstance (cfg, str) : nml_cfg = f90nml.read (cfg) 1688 if isinstance (cfg, f90nml.namelist.Namelist) : nml_cfg = cfg 1689 1690 if out == 'dict' : dict_namelist = {} 1691 if out == 'xr' : xr_namelist = xr.Dataset () 1692 1693 list_nml = [] ; list_comment = [] 1694 1695 if ref != None : 1696 list_nml.append (nml_ref) ; list_comment.append ('ref') 1697 if cfg != None : 1698 list_nml.append (nml_cfg) ; list_comment.append ('cfg') 1699 1700 for nml, comment in zip (list_nml, list_comment) : 1701 if verbose : print (comment) 1702 if flat and out =='dict' : 1703 for nam in nml.keys () : 1704 if verbose : print (nam) 1705 for value in nml[nam] : 1706 if out == 'dict' : dict_namelist[value] = nml[nam][value] 1707 if verbose : print (nam, ':', value, ':', nml[nam][value]) 1708 else : 1709 for nam in nml.keys () : 1710 if verbose : print (nam) 1711 if out == 'dict' : 1712 if nam not in dict_namelist.keys () : dict_namelist[nam] = {} 1713 for value in nml[nam] : 1714 if out == 'dict' : dict_namelist[nam][value] = nml[nam][value] 1715 if out == 'xr' : xr_namelist[value] = nml[nam][value] 1716 if verbose : print (nam, ':', value, ':', nml[nam][value]) 1717 1718 if out == 'dict' : return dict_namelist 1719 if out == 'xr' : return xr_namelist 1720 2352 flat : only for dict output. Output a flat dictionary with all values. 2353 2354 ''' 2355 if ref : 2356 if isinstance (ref, str) : 2357 nml_ref = f90nml.read (ref) 2358 if isinstance (ref, f90nml.namelist.Namelist) : 2359 nml_ref = ref 2360 2361 if cfg : 2362 if isinstance (cfg, str) : 2363 nml_cfg = f90nml.read (cfg) 2364 if isinstance (cfg, f90nml.namelist.Namelist) : 2365 nml_cfg = cfg 2366 2367 if out == 'dict' : 2368 dict_namelist = {} 2369 if out == 'xr' : 2370 xr_namelist = xr.Dataset () 2371 2372 list_nml = [] 2373 list_comment = [] 2374 2375 if ref : 2376 list_nml.append (nml_ref) 2377 list_comment.append ('ref') 2378 if cfg : 2379 list_nml.append (nml_cfg) 2380 list_comment.append ('cfg') 2381 2382 for nml, comment in zip (list_nml, list_comment) : 2383 if verbose : 2384 print (comment) 2385 if flat and out =='dict' : 2386 for nam in nml.keys () : 2387 if verbose : 2388 print (nam) 2389 for value in nml[nam] : 2390 if out == 'dict' : 2391 dict_namelist[value] = nml[nam][value] 2392 if verbose : 2393 print (nam, ':', value, ':', nml[nam][value]) 2394 else : 2395 for nam in nml.keys () : 2396 if verbose : 2397 print (nam) 2398 if out == 'dict' : 2399 if nam not in dict_namelist.keys () : 2400 dict_namelist[nam] = {} 2401 for value in nml[nam] : 2402 if out == 'dict' : 2403 dict_namelist[nam][value] = nml[nam][value] 2404 if out == 'xr' : 2405 xr_namelist[value] = nml[nam][value] 2406 if verbose : 2407 print (nam, ':', value, ':', nml[nam][value]) 2408 2409 if out == 'dict' : 2410 return dict_namelist 2411 if out == 'xr' : 2412 return xr_namelist 2413 else : 2414 def namelist_read (ref=None, cfg=None, out='dict', flat=False, verbose=False) : 2415 '''Shadow version of namelist read, when f90nm module was not found 2416 2417 namelist_read : 2418 Read NEMO namelist(s) and return either a dictionnary or an xarray dataset 2419 ''' 2420 print ( 'Error : module f90nml not found' ) 2421 print ( 'Cannot call namelist_read' ) 2422 print ( 'Call parameters where : ') 2423 print ( f'{err=} {ref=} {cfg=} {out=} {flat=} {verbose=}' ) 1721 2424 1722 2425 def fill_closed_seas (imask, nperio=None, cd_type='T') : 1723 2426 '''Fill closed seas with image processing library 2427 1724 2428 imask : mask, 1 on ocean, 0 on land 1725 2429 ''' … … 1730 2434 1731 2435 return imask_filled 2436 2437 # ====================================================== 2438 # Sea water state function parameters from NEMO code 2439 2440 RDELTAS = 32. 2441 R1_S0 = 0.875/35.16504 2442 R1_T0 = 1./40. 2443 R1_Z0 = 1.e-4 2444 2445 EOS000 = 8.0189615746e+02 2446 EOS100 = 8.6672408165e+02 2447 EOS200 = -1.7864682637e+03 2448 EOS300 = 2.0375295546e+03 2449 EOS400 = -1.2849161071e+03 2450 EOS500 = 4.3227585684e+02 2451 EOS600 = -6.0579916612e+01 2452 EOS010 = 2.6010145068e+01 2453 EOS110 = -6.5281885265e+01 2454 EOS210 = 8.1770425108e+01 2455 EOS310 = -5.6888046321e+01 2456 EOS410 = 1.7681814114e+01 2457 EOS510 = -1.9193502195 2458 EOS020 = -3.7074170417e+01 2459 EOS120 = 6.1548258127e+01 2460 EOS220 = -6.0362551501e+01 2461 EOS320 = 2.9130021253e+01 2462 EOS420 = -5.4723692739 2463 EOS030 = 2.1661789529e+01 2464 EOS130 = -3.3449108469e+01 2465 EOS230 = 1.9717078466e+01 2466 EOS330 = -3.1742946532 2467 EOS040 = -8.3627885467 2468 EOS140 = 1.1311538584e+01 2469 EOS240 = -5.3563304045 2470 EOS050 = 5.4048723791e-01 2471 EOS150 = 4.8169980163e-01 2472 EOS060 = -1.9083568888e-01 2473 EOS001 = 1.9681925209e+01 2474 EOS101 = -4.2549998214e+01 2475 EOS201 = 5.0774768218e+01 2476 EOS301 = -3.0938076334e+01 2477 EOS401 = 6.6051753097 2478 EOS011 = -1.3336301113e+01 2479 EOS111 = -4.4870114575 2480 EOS211 = 5.0042598061 2481 EOS311 = -6.5399043664e-01 2482 EOS021 = 6.7080479603 2483 EOS121 = 3.5063081279 2484 EOS221 = -1.8795372996 2485 EOS031 = -2.4649669534 2486 EOS131 = -5.5077101279e-01 2487 EOS041 = 5.5927935970e-01 2488 EOS002 = 2.0660924175 2489 EOS102 = -4.9527603989 2490 EOS202 = 2.5019633244 2491 EOS012 = 2.0564311499 2492 EOS112 = -2.1311365518e-01 2493 EOS022 = -1.2419983026 2494 EOS003 = -2.3342758797e-02 2495 EOS103 = -1.8507636718e-02 2496 EOS013 = 3.7969820455e-01 2497 2498 def rhop ( ptemp, psal ) : 2499 '''Returns potential density referenced to surface 2500 2501 Computation from NEMO code 2502 ''' 2503 zt = ptemp * R1_T0 # Temperature (°C) 2504 zs = np.sqrt ( np.abs( psal + RDELTAS ) * R1_S0 ) # Square root of salinity (PSS) 2505 # 2506 prhop = ( 2507 (((((EOS060*zt 2508 + EOS150*zs + EOS050)*zt 2509 + (EOS240*zs + EOS140)*zs + EOS040)*zt 2510 + ((EOS330*zs + EOS230)*zs + EOS130)*zs + EOS030)*zt 2511 + (((EOS420*zs + EOS320)*zs + EOS220)*zs + EOS120)*zs + EOS020)*zt 2512 + ((((EOS510*zs + EOS410)*zs + EOS310)*zs + EOS210)*zs + EOS110)*zs + EOS010)*zt 2513 + (((((EOS600*zs+ EOS500)*zs + EOS400)*zs + EOS300)*zs + EOS200)*zs + EOS100)*zs + EOS000 ) 2514 # 2515 return prhop 2516 2517 def rho ( pdep, ptemp, psal ) : 2518 '''Returns in situ density 2519 2520 Computation from NEMO code 2521 ''' 2522 zh = pdep * R1_Z0 # Depth (m) 2523 zt = ptemp * R1_T0 # Temperature (°C) 2524 zs = np.sqrt ( np.abs( psal + RDELTAS ) * R1_S0 ) # Square root salinity (PSS) 2525 # 2526 zn3 = EOS013*zt + EOS103*zs+EOS003 2527 # 2528 zn2 = (EOS022*zt + EOS112*zs+EOS012)*zt + (EOS202*zs+EOS102)*zs+EOS002 2529 # 2530 zn1 = ( 2531 (((EOS041*zt 2532 + EOS131*zs + EOS031)*zt 2533 + (EOS221*zs + EOS121)*zs + EOS021)*zt 2534 + ((EOS311*zs + EOS211)*zs + EOS111)*zs + EOS011)*zt 2535 + (((EOS401*zs + EOS301)*zs + EOS201)*zs + EOS101)*zs + EOS001 ) 2536 # 2537 zn0 = ( 2538 (((((EOS060*zt 2539 + EOS150*zs + EOS050)*zt 2540 + (EOS240*zs + EOS140)*zs + EOS040)*zt 2541 + ((EOS330*zs + EOS230)*zs + EOS130)*zs + EOS030)*zt 2542 + (((EOS420*zs + EOS320)*zs + EOS220)*zs + EOS120)*zs + EOS020)*zt 2543 + ((((EOS510*zs + EOS410)*zs + EOS310)*zs + EOS210)*zs + EOS110)*zs + EOS010)*zt 2544 + (((((EOS600*zs + EOS500)*zs + EOS400)*zs + EOS300)*zs + 2545 EOS200)*zs + EOS100)*zs + EOS000 ) 2546 # 2547 prho = ( ( zn3 * zh + zn2 ) * zh + zn1 ) * zh + zn0 2548 # 2549 return prho 1732 2550 1733 2551 ## =========================================================================== … … 1737 2555 ## =========================================================================== 1738 2556 1739 def __is_orca_north_fold__ ( Xtest, cname_long='T' ) :1740 '''1741 Ported (pirated !!?) from Sosie1742 1743 Tell if there is a 2/point band overlaping folding at the north pole typical of the ORCA grid1744 1745 0 => not an orca grid (or unknown one)1746 4 => North fold T-point pivot (ex: ORCA2)1747 6 => North fold F-point pivot (ex: ORCA1)1748 1749 We need all this 'cname_long' stuff because with our method, there is a1750 confusion between "Grid_U with T-fold" and "Grid_V with F-fold"1751 => so knowing the name of the longitude array (as in namelist, and hence as1752 in netcdf file) might help taking the righ decision !!! UGLY!!!1753 => not implemented yet1754 '''1755 1756 ifld_nord = 0 ; cgrd_type = 'X'1757 ny, nx = Xtest.shape[-2:]1758 1759 if ny > 3 : # (case if called with a 1D array, ignoring...)1760 if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. :1761 ifld_nord = 4 ; cgrd_type = 'T' # T-pivot, grid_T 1762 1763 if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2 :-1] ).sum() == 0. :1764 if cnlon == 'U' : ifld_nord = 4 ; cgrd_type = 'U' # T-pivot, grid_T1765 ## LOLO: PROBLEM == 6, V !!!1766 1767 if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. :1768 ifld_nord = 4 ; cgrd_type = 'V' # T-pivot, grid_V1769 1770 if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-2, nx-1-1:nx-nx//2:-1] ).sum() == 0. :1771 ifld_nord = 6 ; cgrd_type = 'T'# F-pivot, grid_T1772 1773 if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-1, nx-1:nx-nx//2-1:-1] ).sum() == 0. :1774 ifld_nord = 6 ; cgrd_type = 'U' # F-pivot, grid_U1775 1776 if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2 :-1] ).sum() == 0. :1777 if cnlon == 'V' : ifld_nord = 6 ; cgrd_type = 'V' # F-pivot, grid_V1778 ## LOLO: PROBLEM == 4, U !!!1779 1780 return ifld_nord, cgrd_type2557 # def __is_orca_north_fold__ ( Xtest, cname_long='T' ) : 2558 # ''' 2559 # Ported (pirated !!?) from Sosie 2560 2561 # Tell if there is a 2/point band overlaping folding at the north pole typical of the ORCA grid 2562 2563 # 0 => not an orca grid (or unknown one) 2564 # 4 => North fold T-point pivot (ex: ORCA2) 2565 # 6 => North fold F-point pivot (ex: ORCA1) 2566 2567 # We need all this 'cname_long' stuff because with our method, there is a 2568 # confusion between "Grid_U with T-fold" and "Grid_V with F-fold" 2569 # => so knowing the name of the longitude array (as in namelist, and hence as 2570 # in netcdf file) might help taking the righ decision !!! UGLY!!! 2571 # => not implemented yet 2572 # ''' 2573 2574 # ifld_nord = 0 ; cgrd_type = 'X' 2575 # ny, nx = Xtest.shape[-2:] 2576 2577 # if ny > 3 : # (case if called with a 1D array, ignoring...) 2578 # if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. : 2579 # ifld_nord = 4 ; cgrd_type = 'T' # T-pivot, grid_T 2580 2581 # if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2 :-1] ).sum() == 0. : 2582 # if cnlon == 'U' : ifld_nord = 4 ; cgrd_type = 'U' # T-pivot, grid_T 2583 # ## LOLO: PROBLEM == 6, V !!! 2584 2585 # if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-1:nx-nx//2+1:-1] ).sum() == 0. : 2586 # ifld_nord = 4 ; cgrd_type = 'V' # T-pivot, grid_V 2587 2588 # if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-2, nx-1-1:nx-nx//2:-1] ).sum() == 0. : 2589 # ifld_nord = 6 ; cgrd_type = 'T'# F-pivot, grid_T 2590 2591 # if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-1, nx-1:nx-nx//2-1:-1] ).sum() == 0. : 2592 # ifld_nord = 6 ; cgrd_type = 'U' # F-pivot, grid_U 2593 2594 # if ( Xtest [ny-1, 1:nx//2-1] - Xtest [ny-3, nx-2:nx-nx//2 :-1] ).sum() == 0. : 2595 # if cnlon == 'V' : ifld_nord = 6 ; cgrd_type = 'V' # F-pivot, grid_V 2596 # ## LOLO: PROBLEM == 4, U !!! 2597 2598 # return ifld_nord, cgrd_type -
TOOLS/MOSAIX/update_xml.py
r6190 r6666 26 26 ## software by incorrectly or partially configured 27 27 ## 28 ## 28 ''' 29 Python script used to perform on the fly 30 editing of xml files. Used by `CreateWeights.bash`. More 31 information with `python update_xml.py -h`. 32 29 33 ## SVN information 30 __Author__ = "$Author$" 31 __Date__ = "$Date$" 32 __Revision__ = "$Revision$" 33 __Id__ = "$Id$" 34 __HeadURL = "$HeadURL$" 34 Author = "$Author$" 35 Date = "$Date$" 36 Revision = "$Revision$" 37 Id = "$Id$" 38 HeadURL = "$HeadURL$" 39 ''' 40 41 ## SVN information 42 __SVN__ = ({ 43 'Author' : "$Author$", 44 'Date' : "$Date$", 45 'Revision' : "$Revision$", 46 'Id' : "$Id$", 47 'HeadURL' : "$HeadURL$", 48 }) 35 49 36 50 # python update_xml.py -i ~/Unix/TOOLS/MOSAIX/iodef_atm_to_oce.xml -o essai.xml -n 'context[@id="interpol_read"]/file_definition/file[@id="file_src"]/field[@id="mask_src"]' -k name=Bidon … … 41 55 # 42 56 import xml.etree.ElementTree 43 import argparse, sys, textwrap,shlex44 57 import argparse, sys, shlex 58 45 59 # Check version of Python 46 60 Version = sys.version_info 47 if Version < (2,7,0) :48 sys.stderr.write ( "You need python 2.7 or later to run this script\n" )49 sys.stderr.write ( "Present version is: " + str(Version[0]) + "." + str(Version[1]) + "." + str(Version[2]) + "\n" )50 sys.exit (1)51 61 52 62 ## ============================================================================ … … 56 66 '''Concatenate some elements of the list of strings when needed''' 57 67 zlist = list_str.copy () ; list_new = [] 58 while ( len (zlist) > 0 ):59 arg = zlist.pop (0) 68 while len (zlist) > 0 : 69 arg = zlist.pop (0) 60 70 if arg[0] == '"' : 61 71 for arg2 in zlist.copy () : … … 66 76 list_new.append (arg) 67 77 return list_new 68 69 def UpdateNode ( iodef, Node, Text=None,Key=None) :78 79 def UpdateNode (piodef, pNode, pText=None, pKey=None) : 70 80 '''Update an xml node''' 71 81 # Remove whitespaces at both ends 72 Node =Node.rstrip().lstrip()82 zNode = pNode.rstrip().lstrip() 73 83 74 84 ## Find node 75 nodeList = iodef.findall (Node)85 nodeList = piodef.findall (zNode) 76 86 77 87 ## Check that one and only one node is found 78 88 if len (nodeList) == 0 : 79 89 print ( "Error : node not found" ) 80 print ( "Node : ", Node )90 print ( "Node : ", zNode ) 81 91 sys.exit (1) 82 92 83 93 if len (nodeList) > 1 : 84 94 print ( "Error : " + len (nodeList)+" occurences of node found in file" ) 85 print ( "Node : ", Node )95 print ( "Node : ", zNode ) 86 96 sys.exit (2) 87 97 … … 89 99 elem = nodeList[0] 90 100 91 if Text !=None :92 if Verbose : print ( 'Node:', Node, ' -- Text:',Text )101 if pText is not None : 102 if Verbose : print ( 'Node:', zNode, ' -- Text:', pText ) 93 103 if Debug : 94 104 print ( 'Attributes of node: ' + str (elem.attrib) ) 95 105 print ( 'Text : ' + str (elem.text) ) 96 elem.text = Text97 98 if Key !=None :99 106 elem.text = pText 107 108 if pKey is not None : 109 100 110 # Check the syntax 101 if not '=' in Key :111 if not '=' in pKey : 102 112 print ( 'Key syntax error. Correct syntax is -k Key=Value' ) 103 113 sys.exit (-1) 104 114 else : 105 KeyName, Value = Key.split ('=')106 if Verbose : print ( 'Node:', Node, ' -- Key:',Key )115 KeyName, Value = pKey.split ('=') 116 if Verbose : print ( 'Node:', zNode, ' -- Key:', pKey ) 107 117 # To do : check that KeyName exist (it is added if not : do we want that ?) 108 118 if Debug : … … 110 120 elem.attrib.update ( { KeyName:Value } ) 111 121 112 return iodef122 return piodef 113 123 114 124 ## ============================================================================ … … 118 128 # The first step in using the argparse is creating an ArgumentParser object: 119 129 parser = argparse.ArgumentParser (description = """ 120 Examples with the modification on the command line : 130 Examples with the modification on the command line : 121 131 python %(prog)s -i iodef.xml -n 'context[@id="interpol_run"]/file_definition/file[@id="file_src"]/field[@id="mask_source"]' -k name -v maskutil_T 122 132 python %(prog)s -i iodef.xml -n 'context[@id="interpol_run"]/file_definition/file[@id="dia"]/variable[@name="dest_grid"]' -t dstDomainType 123 133 124 Usage with a command file : 134 Usage with a command file : 125 135 python %(prog)s -i iodef.xml -c commands.txt 126 127 Syntax in the command file : removes the quote around the node description : 136 137 Syntax in the command file : removes the quote around the node description : 128 138 -n context[@id="interpol_run"]/file_definition/file[@id="dia"]/variable[@name="dest_grid"] -t dstDomainType 129 139 130 """ + "SVN : " + __ Revision__, formatter_class=argparse.RawDescriptionHelpFormatter, epilog='-------- This is the end of the help message --------')140 """ + "SVN : " + __SVN__['Revision'], formatter_class=argparse.RawDescriptionHelpFormatter, epilog='-------- This is the end of the help message --------') 131 141 132 142 # Adding arguments … … 155 165 Key = myargs.key 156 166 Text = myargs.text 157 158 if FileCommand != None :159 if ( Node != None or Key != None or Text !=None ) :167 168 if FileCommand is not None : 169 if ( Node is not None or Key is not None or Text is not None ) : 160 170 print ('Error : when a command file is specified, options -k|--key, -n|--node are unused' ) 161 171 exit (-2) 162 163 if FileOut ==None : FileOut = FileIn164 172 173 if FileOut is None : FileOut = FileIn 174 165 175 ## Get XML tree from input file 166 176 iodef = xml.etree.ElementTree.parse ( FileIn ) 167 177 168 if FileCommand ==None :178 if FileCommand is None : 169 179 ## Only one node to modify 170 180 iodef = UpdateNode (iodef, Node, Key, Text) … … 172 182 else : 173 183 ## Read a list of modification commands in a command file 174 fic = open (FileCommand, 'r' )175 lignes = fic.readlines ()184 fic = open (FileCommand, 'r', encoding='utf-8') 185 lignes = fic.readlines () 176 186 177 187 for nn, ligne in enumerate (lignes) : 178 188 ligne = ligne.strip ().split ('#')[0] # Remove leading and trailing blanks, and trailing comments 179 189 if Debug : print (nn+1, ':', type (ligne) , ':', len (ligne) , ':', ligne, ':') 180 if ligne == '' or ligne ==None or len(ligne) == 0 or ligne[0] == '#' :190 if ligne == '' or ligne is None or len(ligne) == 0 or ligne[0] == '#' : 181 191 if Debug : print ('Skips blank or comment line') 182 else : 192 else : 183 193 list_args = shlex.split (ligne) 184 194 185 195 if Debug : 186 196 print ( '{:3d} : '.format(nn+1), end='') … … 193 203 Text = myargs_ligne.text 194 204 195 UpdateNode (iodef, Node, Text, Key)196 205 iodef = UpdateNode (iodef, Node, Text, Key) 206 197 207 ## Writes XML tree to file 198 208 iodef.write ( FileOut ) … … 201 211 if Debug : print ('This is the end') 202 212 #sys.exit (0) 203 213 204 214 ### =========================================================================== 205 215 ### … … 207 217 ### 208 218 ### =========================================================================== 209
Note: See TracChangeset
for help on using the changeset viewer.