Changeset 6066 for TOOLS/MOSAIX
- Timestamp:
- 02/24/22 15:17:45 (2 years ago)
- Location:
- TOOLS/MOSAIX
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
TOOLS/MOSAIX/CalvingWeights.py
r6064 r6066 14 14 ## personal. 15 15 ## 16 import n etCDF4, numpy as np16 import numpy as np, xarray as xr 17 17 import sys, os, platform, argparse, textwrap, time 18 18 from scipy import ndimage … … 36 36 37 37 # Adding arguments 38 parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025', 'eORCA025.1'] ) 38 parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', 39 choices=['ORCA2.3', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 39 40 parser.add_argument ('--atm' , help='atm model name (ICO* or LMD*)', type=str, default='ICO40' ) 40 41 parser.add_argument ('--type' , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full' ) 41 parser.add_argument ('--dir' , help='Directory of input file', type=str, default='./' )42 ## parser.add_argument ('--dir' , help='Directory of input file', type=str, default='./' ) 42 43 parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' ) 43 44 parser.add_argument ('--repartition_var' , help='Variable name for iceshelf' , type=str, default=None) 45 parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) 46 parser.add_argument ('--areas' , help='masks file name', default='areas.nc' ) 47 parser.add_argument ('--masks' , help='areas file name', default='masks.nc' ) 48 parser.add_argument ('--o2a' , help='o2a file name' , default='o2a.nc' ) 44 49 parser.add_argument ('--output' , help='output rmp file name', default='rmp_tlmd_to_torc_calving.nc' ) 45 parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 46 parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=int, default=0 ) 50 parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', 51 choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 52 parser.add_argument ('--ocePerio', help='periodicity of ocean grid', type=int, default=0 ) 47 53 48 54 # Parse command line … … 89 95 90 96 # Reading domains characteristics 91 areas = myargs.dir + '/areas_' + dst_Name + 'x' + src_Name + '.nc' 92 masks = myargs.dir + '/masks_' + dst_Name + 'x' + src_Name + '.nc' 93 grids = myargs.dir + '/grids_' + dst_Name + 'x' + src_Name + '.nc' 97 grids = myargs.grids 98 areas = myargs.areas 99 masks = myargs.masks 100 o2a = myargs.o2a 94 101 95 102 print ('Opening ' + areas) 96 f_areas = netCDF4.Dataset ( areas, "r")103 f_areas = xr.open_dataset ( areas ) 97 104 print ('Opening ' + masks) 98 f_masks = netCDF4.Dataset ( masks, "r")105 f_masks = xr.open_dataset ( masks ) 99 106 print ('Opening ' + grids) 100 f_grids = netCDF4.Dataset ( grids, "r")101 102 src_lon = f_grids .variables ['t' + src_name + '.lon'][:]103 src_lat = f_grids .variables ['t' + src_name + '.lat'][:]104 dst_lon = f_grids .variables ['t' + dst_name + '.lon'][:]105 dst_lat = f_grids .variables ['t' + dst_name + '.lat'][:]106 107 src_msk = f_masks .variables ['t' + src_name + '.msk'][:]108 dst_msk = f_masks .variables ['t' + dst_name + '.msk'][:]107 f_grids = xr.open_dataset ( grids ) 108 109 src_lon = f_grids ['t' + src_name + '.lon'] 110 src_lat = f_grids ['t' + src_name + '.lat'] 111 dst_lon = f_grids ['t' + dst_name + '.lon'] 112 dst_lat = f_grids ['t' + dst_name + '.lat'] 113 114 src_msk = f_masks ['t' + src_name + '.msk'] 115 dst_msk = f_masks ['t' + dst_name + '.msk'] 109 116 dst_mskutil = 1-dst_msk # Reversing the convention : 0 on continent, 1 on ocean 110 print ('dst_msk : ' + str(np.sum(dst_msk)))111 print ('dst_mskutil : ' + str(np.sum(dst_mskutil)))117 print ('dst_msk : ', dst_msk.sum().values ) 118 print ('dst_mskutil : ', dst_mskutil.sum().values ) 112 119 113 120 # Ocean grid periodicity … … 116 123 # Periodicity masking for NEMO 117 124 if nperio_dst == 0 : 118 if dst_Name == 'ORCA2.3' : nperio_dst = 4 119 if dst_Name == 'eORCA1.2' : nperio_dst = 6 120 if dst_Name == 'ORCA025' : nperio_dst = 6 121 if dst_Name == 'eORCA025' : nperio_dst = 6 122 if dst_Name == 'eORCA025.1' : nperio_dst = 6 125 if dst_Name in ['ORCA2.3', 'ORCA2.4'] : nperio_dst = 4 126 if dst_Name in ['eORCA1.2', 'eORCA1.3'] : nperio_dst = 6 127 if dst_Name in ['ORCA025', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] : nperio_dst = 6 128 123 129 print ('nperio_dst: ' + str(nperio_dst) ) 124 130 dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T', sval=0 ) … … 144 150 dst_mskutil[997:1012,907] = 0 145 151 146 if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4' ] :152 if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4', 'eORCA1.4.2' ] : 147 153 print ('Filling some seas for ' + dst_Name) 148 154 # Set Gibraltar strait to 0 to fill Mediterrean sea … … 184 190 dst_srfutil_sum = np.sum( dst_srfutil) 185 191 186 src_clo = f_grids .variables ['t' + src_name + '.clo'][:]187 src_cla = f_grids .variables ['t' + src_name + '.cla'][:]188 dst_clo = f_grids .variables ['t' + dst_name + '.clo'][:]189 dst_cla = f_grids .variables ['t' + dst_name + '.cla'][:]192 src_clo = f_grids ['t' + src_name + '.clo'] 193 src_cla = f_grids ['t' + src_name + '.cla'] 194 dst_clo = f_grids ['t' + dst_name + '.clo'] 195 dst_cla = f_grids ['t' + dst_name + '.cla'] 190 196 191 197 # Indices … … 198 204 # Reading data file for calving or iceberg geometry around Antarctica 199 205 print ( 'Opening ' + myargs.repartition_file) 200 f_repartition = netCDF4.Dataset ( myargs.repartition_file, "r")201 repartition = np.sum ( f_repartition.variables [repartitionVar] [:], axis=0 )206 f_repartition = xr.open_dataset ( myargs.repartition_file ) 207 repartition = np.sum ( f_repartition.variables [repartitionVar], axis=0 ) 202 208 203 209 ## Before loop on basins … … 246 252 247 253 # Integral and normalisation 248 sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil ) 254 sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil ).values 249 255 key_repartition = key_repartition / sum_repartition 250 256 251 257 print ( 'Sum of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] )) ) 252 print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] * dst_srfutil ) ) )258 print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] * dst_srfutil ).values) ) 253 259 254 260 # Basin surface (modulated by repartition key) … … 257 263 # Weights and links 258 264 poids = 1.0 / basin_srfutil 259 matrix_local = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids , 0. )265 matrix_local = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids.values, 0. ) 260 266 matrix_local = matrix_local[matrix_local > 0.0] # Keep only non zero values 261 267 num_links = np.count_nonzero (matrix_local) … … 297 303 # Output file 298 304 calving = myargs.output 299 f_calving = netCDF4.Dataset ( calving, 'w', format=FmtNetcdf )300 305 301 306 print ('Output file: ' + calving ) 302 307 303 f_calving.Conventions = "CF-1.6" 304 f_calving.source = "IPSL Earth system model" 305 f_calving.group = "ICMC IPSL Climate Modelling Center" 306 f_calving.Institution = "IPSL https.//www.ipsl.fr" 307 f_calving.Ocean = dst_Name + " https://www.nemo-ocean.eu" 308 f_calving.Atmosphere = src_Name + " http://lmdz.lmd.jussieu.fr" 309 if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.repartition_file 310 f_calving.associatedFiles = grids + " " + areas + " " + masks 311 f_calving.directory = os.getcwd () 312 f_calving.description = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX" 313 f_calving.title = calving 314 f_calving.Program = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) 315 f_calving.repartitionType = myargs.type 316 if myargs.type in [ 'iceberg', 'iceshelf' ] : 317 f_calving.repartitionFile = myargs.repartition_file 318 f_calving.repartitionVar = repartitionVar 319 f_calving.gridsFile = grids 320 f_calving.areasFile = areas 321 f_calving.masksFile = masks 322 f_calving.timeStamp = time.asctime() 323 f_calving.HOSTNAME = platform.node() 324 #f_calving.LOGNAME = os.getlogin() 325 f_calving.Python = "Python version " + platform.python_version() 326 f_calving.OS = platform.system() 327 f_calving.release = platform.release() 328 f_calving.hardware = platform.machine() 329 f_calving.conventions = "SCRIP" 330 if src_name == 'lmd' : f_calving.source_grid = "curvilinear" 331 if src_name == 'ico' : f_calving.source_grid = "unstructured" 332 f_calving.dest_grid = "curvilinear" 333 f_calving.Model = "IPSL CM6" 334 f_calving.SVN_Author = "$Author$" 335 f_calving.SVN_Date = "$Date$" 336 f_calving.SVN_Revision = "$Revision$" 337 f_calving.SVN_Id = "$Id$" 338 f_calving.SVN_HeadURL = "$HeadURL$" 339 340 d_nb_zone = f_calving.createDimension ('nb_zone' , nb_zone ) 341 d_num_links = f_calving.createDimension ('num_links' , num_links ) 342 d_num_wgts = f_calving.createDimension ('num_wgts' , 1 ) 343 344 d_src_grid_size = f_calving.createDimension ('src_grid_size' , src_grid_size ) 345 d_src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0] ) 346 d_src_grid_rank = f_calving.createDimension ('src_grid_rank' , 2 ) 347 348 d_dst_grid_size = f_calving.createDimension ('dst_grid_size' , dst_grid_size ) 349 d_dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] ) 350 d_dst_grid_rank = f_calving.createDimension ('dst_grid_rank' , 2 ) 351 352 v_remap_matrix = f_calving.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') ) 353 v_src_address = f_calving.createVariable ( 'src_address' , 'i4', ('num_links',) ) 354 v_src_address.convention = "Fortran style addressing, starting at 1" 355 v_dst_address = f_calving.createVariable ( 'dst_address' , 'i4', ('num_links',) ) 356 v_dst_address.convention = "Fortran style addressing, starting at 1" 357 358 v_remap_matrix[:] = remap_matrix 359 v_src_address [:] = src_address 360 v_dst_address [:] = dst_address 361 362 v_src_grid_dims = f_calving.createVariable ( 'src_grid_dims' , 'i4', ('src_grid_rank',) ) 363 v_src_grid_center_lon = f_calving.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) ) 364 v_src_grid_center_lat = f_calving.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) ) 365 v_src_grid_center_lon.units='degrees_east' ; v_src_grid_center_lon.long_name='Longitude' 366 v_src_grid_center_lon.long_name='longitude' ; v_src_grid_center_lon.bounds="src_grid_corner_lon" 367 v_src_grid_center_lat.units='degrees_north' ; v_src_grid_center_lat.long_name='Latitude' 368 v_src_grid_center_lat.long_name='latitude ' ; v_src_grid_center_lat.bounds="src_grid_corner_lat" 369 v_src_grid_corner_lon = f_calving.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners') ) 370 v_src_grid_corner_lat = f_calving.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners') ) 371 v_src_grid_corner_lon.units="degrees_east" 372 v_src_grid_corner_lat.units="degrees_north" 373 v_src_grid_area = f_calving.createVariable ( 'src_grid_area' , 'f8', ('src_grid_size',) ) 374 v_src_grid_area.long_name="Grid area" ; v_src_grid_area.standard_name="cell_area" ; v_src_grid_area.units="m2" 375 v_src_grid_imask = f_calving.createVariable ( 'src_grid_imask' , 'f8', ('src_grid_size',) ) 376 v_src_grid_imask.long_name="Land-sea mask" ; v_src_grid_imask.units="Land:1, Ocean:0" 377 378 v_src_grid_dims [:] = ( src_jpi, src_jpi ) 379 v_src_grid_center_lon[:] = src_lon[:].ravel() 380 v_src_grid_center_lat[:] = src_lat[:].ravel() 381 v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) 382 v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) 383 v_src_grid_area [:] = src_srf[:].ravel() 384 v_src_grid_imask [:] = src_msk[:].ravel() 308 remap_matrix = xr.DataArray (np.reshape(remap_matrix, (num_links, 1)), dims = ['num_links', 'num_wgts'] ) 309 src_address = xr.DataArray (src_address.astype(np.int32) , dims = ['num_links',] ) 310 src_address.attrs['convention'] = "Fortran style addressing, starting at 1" 311 dst_address = xr.DataArray (dst_address.astype(np.int32) , dims = ['num_links',] ) 312 dst_address.attrs['convention'] = "Fortran style addressing, starting at 1" 313 314 src_grid_dims = xr.DataArray ( np.array (( src_jpi, src_jpi ), dtype=np.int32), dims = ['src_grid_rank',] ) 315 src_grid_center_lon = xr.DataArray ( src_lon.values.ravel(), dims = ['src_grid_size',] ) 316 src_grid_center_lat = xr.DataArray ( src_lat.values.ravel(), dims = ['src_grid_size',] ) 317 src_grid_center_lon.attrs['units']='degrees_east' ; src_grid_center_lon.attrs['long_name']='Longitude' 318 src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" 319 src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' 320 src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" 321 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'] ) 322 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' ] ) 323 src_grid_corner_lon.attrs['units']="degrees_east" 324 src_grid_corner_lat.attrs['units']="degrees_north" 325 src_grid_area = xr.DataArray ( src_srf.values.ravel(), 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 src_grid_imask = xr.DataArray ( src_msk.values.ravel().astype(np.int32) , dims = ['src_grid_size',] ) 328 src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" 385 329 386 330 # -- 387 388 v_dst_grid_dims = f_calving.createVariable ( 'dst_grid_dims' , 'i4', ('dst_grid_rank',) ) 389 v_dst_grid_center_lon = f_calving.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) ) 390 v_dst_grid_center_lat = f_calving.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) ) 391 v_dst_grid_center_lon.units='degrees_east' ; v_dst_grid_center_lon.long_name='Longitude' 392 v_dst_grid_center_lon.long_name='longitude' ; v_dst_grid_center_lon.bounds="dst_grid_corner_lon" 393 v_dst_grid_center_lat.units='degrees_north' ; v_dst_grid_center_lat.long_name='Latitude' 394 v_dst_grid_center_lat.long_name='latitude' ; v_dst_grid_center_lat.bounds="dst_grid_corner_lat" 395 v_dst_grid_corner_lon = f_calving.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners') ) 396 v_dst_grid_corner_lat = f_calving.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners') ) 397 v_dst_grid_corner_lon.units="degrees_east" 398 v_dst_grid_corner_lat.units="degrees_north" 399 v_dst_grid_area = f_calving.createVariable ( 'dst_grid_area' , 'f8', ('dst_grid_size',) ) 400 v_dst_grid_area.long_name="Grid area" ; v_dst_grid_area.standard_name="cell_area" ; v_dst_grid_area.units="m2" 401 v_dst_grid_imask = f_calving.createVariable ( 'dst_grid_imask' , 'f8', ('dst_grid_size',) ) 402 v_dst_grid_imask.long_name="Land-sea mask" ; v_dst_grid_imask.units="Land:1, Ocean:0" 403 404 v_dst_bassin = f_calving.createVariable ( 'dst_bassin' , 'f8', ('nb_zone', 'dst_grid_size',) ) 405 v_dst_repartition = f_calving.createVariable ( 'dst_repartition' , 'f8', ('nb_zone', 'dst_grid_size',) ) 406 407 v_dst_southLimit = f_calving.createVariable ('dst_southLimit', 'f4', ('nb_zone',) ) 408 v_dst_northLimit = f_calving.createVariable ('dst_northLimit', 'f4', ('nb_zone',) ) 409 v_dst_southLimit[:] = np.min(limit_lat, axis=(1,) ) 410 v_dst_northLimit[:] = np.max(limit_lat, axis=(1,) ) 411 412 v_dst_grid_dims [:] = ( dst_jpi, dst_jpi ) 413 v_dst_grid_center_lon[:] = dst_lon[:].ravel() 414 v_dst_grid_center_lat[:] = dst_lat[:].ravel() 415 v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) 416 v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) 417 v_dst_grid_area [:] = dst_srf[:].ravel() 418 v_dst_grid_imask [:] = dst_msk[:].ravel() 419 420 v_dst_bassin[:] = basin_msk.reshape ( (nb_zone,dst_grid_size) ) 421 v_dst_repartition[:] = key_repartition.reshape( (nb_zone,dst_grid_size) ) 331 dst_grid_dims = xr.DataArray ( np.array(( dst_jpi, dst_jpi), dtype=np.int32) , dims = ['dst_grid_rank', ]) 332 dst_grid_center_lon = xr.DataArray ( dst_lon.values.ravel(), dims = ['dst_grid_size', ]) 333 dst_grid_center_lat = xr.DataArray ( dst_lat.values.ravel(), dims = ['dst_grid_size', ]) 334 dst_grid_center_lon.attrs['units']='degrees_east' ; dst_grid_center_lon.attrs['long_name']='Longitude' 335 dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" 336 dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' 337 dst_grid_center_lat.attrs['long_name']='latitude' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" 338 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'] ) 339 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'] ) 340 dst_grid_corner_lon.attrs['units']="degrees_east" 341 dst_grid_corner_lat.attrs['units']="degrees_north" 342 dst_grid_area = xr.DataArray ( dst_srf.values.ravel(), dims = ['dst_grid_size', ] ) 343 dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 344 dst_grid_imask = xr.DataArray ( dst_msk.values.ravel(), dims = ['dst_grid_size',] ) 345 dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" 346 347 dst_bassin = xr.DataArray ( basin_msk.reshape ( (nb_zone,dst_grid_size) ) , dims = ['nb_zone', 'dst_grid_size', ] ) 348 dst_repartition = xr.DataArray ( key_repartition.reshape( (nb_zone,dst_grid_size) ) , dims = ['nb_zone', 'dst_grid_size', ] ) 349 350 dst_southLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , dims = ['nb_zone', ] ) 351 dst_northLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , dims = ['nb_zone', ] ) 422 352 423 353 # Additionnal fields for debugging 424 354 # ================================ 425 v_dst_grid_imaskutil = f_calving.createVariable ( 'dst_grid_imaskutil' , 'f8', ('dst_grid_size',) ) 426 v_dst_grid_imaskutil.long_name="Land-sea mask" ; v_dst_grid_imaskutil.units="Land:1, Ocean:0" 427 v_dst_grid_imaskclose = f_calving.createVariable ( 'dst_grid_imaskclose' , 'f8', ('dst_grid_size',) ) 428 v_dst_grid_imaskclose.long_name="Land-sea mask" ; v_dst_grid_imaskclose.units="Land:1, Ocean:0" 429 v_dst_grid_imaskutil [:] = dst_mskutil[:].ravel() 430 v_dst_grid_imaskclose[:] = dst_closed[:].ravel() 431 432 v_dst_lon_addressed = f_calving.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) ) 433 v_dst_lat_addressed = f_calving.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) ) 434 v_dst_lon_addressed.long_name="Longitude" ; v_dst_lon_addressed.standard_name="longitude" ; v_dst_lon_addressed.units="degrees_east" 435 v_dst_lat_addressed.long_name="Latitude" ; v_dst_lat_addressed.standard_name="latitude" ; v_dst_lat_addressed.units="degrees_north" 436 v_dst_area_addressed = f_calving.createVariable ( 'dst_area_addressed', 'f8', ('num_links',) ) 437 v_dst_area_addressed.long_name="Grid area" ; v_dst_area_addressed.standard_name="cell_area" ; v_dst_grid_area.units="m2" 438 v_dst_imask_addressed = f_calving.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) ) 439 v_dst_imask_addressed.long_name="Land-sea mask" ; v_dst_imask_addressed.units="Land:1, Ocean:0" 440 v_dst_imaskutil_addressed = f_calving.createVariable ( 'dst_imaskutil_addressed', 'i4', ('num_links',) ) 441 v_dst_imaskutil_addressed.long_name="Land-sea mask" ; v_dst_imaskutil_addressed.units="Land:1, Ocean:0" 442 443 v_dst_lon_addressed [:] = dst_lon.ravel()[dst_address-1].ravel() 444 v_dst_lat_addressed [:] = dst_lat.ravel()[dst_address-1].ravel() 445 v_dst_area_addressed [:] = dst_srf[:].ravel()[dst_address-1].ravel() 446 v_dst_imask_addressed[:] = dst_msk[:].ravel()[dst_address-1].ravel() 447 v_dst_imaskutil_addressed[:] = dst_mskutil.ravel()[dst_address-1].ravel() 448 449 v_src_field = f_calving.createVariable ( 'src_field', 'f8', ('src_grid_size',) ) 450 v_dst_field = f_calving.createVariable ( 'dst_field', 'f8', ('dst_grid_size',) ) 451 v_src_field[:] = src_field 452 v_dst_field[:] = dst_field 453 454 f_calving.close () 455 456 457 ### ===== That's all Folk's !! ============================================== 355 dst_grid_imaskutil = xr.DataArray ( dst_mskutil.ravel().astype(np.float32), dims = ['dst_grid_size',] ) 356 dst_grid_imaskutil.attrs['long_name']="Land-sea mask" ; dst_grid_imaskutil.attrs['units']="Land:1, Ocean:0" 357 dst_grid_imaskclose = xr.DataArray ( dst_closed.values.ravel(), dims = ['dst_grid_size',] ) 358 dst_grid_imaskclose.attrs['long_name']="Land-sea mask" ; dst_grid_imaskclose.attrs['units']="Land:1, Ocean:0" 359 360 dst_lon_addressed = xr.DataArray ( dst_lon.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 361 dst_lat_addressed = xr.DataArray ( dst_lat.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 362 dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" 363 dst_lat_addressed.attrs['long_name']="Latitude" ; dst_lat_addressed.attrs['standard_name']="latitude" ; dst_lat_addressed.attrs['units']="degrees_north" 364 dst_area_addressed = xr.DataArray ( dst_srf.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 365 dst_area_addressed.attrs['long_name']="Grid area" ; dst_area_addressed.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 366 dst_imask_addressed = xr.DataArray ( dst_msk.values.ravel()[dst_address-1].ravel(), dims = ['num_links', ] ) 367 dst_imask_addressed.attrs['long_name']="Land-sea mask" ; dst_imask_addressed.attrs['units']="Land:1, Ocean:0" 368 dst_imaskutil_addressed = xr.DataArray ( dst_mskutil.ravel()[dst_address-1].astype(np.float32), dims = ['num_links'] ) 369 dst_imaskutil_addressed.attrs['long_name']="Land-sea mask" ; dst_imaskutil_addressed.attrs['units']="Land:1, Ocean:0" 370 371 src_field = xr.DataArray ( src_field, dims = ['src_grid_size',] ) 372 dst_field = xr.DataArray ( dst_field, dims = ['dst_grid_size',] ) 373 374 f_calving = xr.Dataset ( { 375 'remap_matrix' : remap_matrix, 376 'src_address' : src_address, 377 'dst_address' : dst_address, 378 'src_grid_dims' : src_grid_dims, 379 'src_grid_center_lon' : src_grid_center_lat, 380 'src_grid_center_lat' : src_grid_center_lat, 381 'src_grid_corner_lon' : src_grid_corner_lon, 382 'src_grid_corner_lat' : src_grid_corner_lat, 383 'src_grid_area' : src_grid_area, 384 'src_grid_imask' : src_grid_imask, 385 'dst_grid_dims' : dst_grid_dims, 386 'dst_grid_center_lon' : dst_grid_center_lon, 387 'dst_grid_center_lat' : dst_grid_center_lat, 388 'dst_grid_corner_lon' : dst_grid_corner_lon, 389 'dst_grid_corner_lat' : dst_grid_corner_lat, 390 'dst_grid_area' : dst_grid_area, 391 'dst_grid_imask' : dst_grid_imask, 392 'dst_bassin' : dst_bassin, 393 'dst_repartition' : dst_repartition, 394 'dst_southLimit' : dst_southLimit, 395 'dst_northLimit' : dst_northLimit, 396 'dst_grid_imaskutil' : dst_grid_imaskutil, 397 'dst_grid_imaskclose' : dst_grid_imaskclose, 398 'dst_lon_addressed' : dst_lon_addressed, 399 'dst_lat_addressed' : dst_lat_addressed, 400 'dst_area_addressed' : dst_area_addressed, 401 'dst_imask_addressed' : dst_imask_addressed, 402 'dst_imaskutil_addressed' : dst_imaskutil_addressed, 403 'src_field' : src_field, 404 'dst_field' : dst_field, 405 } ) 406 407 f_calving.attrs['Conventions'] = "CF-1.6" 408 f_calving.attrs['source'] = "IPSL Earth system model" 409 f_calving.attrs['group'] = "ICMC IPSL Climate Modelling Center" 410 f_calving.attrs['Institution'] = "IPSL https.//www.ipsl.fr" 411 f_calving.attrs['Ocean'] = dst_Name + " https://www.nemo-ocean.eu" 412 f_calving.attrs['Atmosphere'] = src_Name + " http://lmdz.lmd.jussieu.fr" 413 if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.attrs['originalFiles'] = myargs.repartition_file 414 f_calving.attrs['associatedFiles'] = grids + " " + areas + " " + masks 415 f_calving.attrs['directory'] = os.getcwd () 416 f_calving.attrs['description'] = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX" 417 f_calving.attrs['title'] = calving 418 f_calving.attrs['Program'] = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) 419 f_calving.attrs['repartitionType'] = myargs.type 420 if myargs.type in [ 'iceberg', 'iceshelf' ] : 421 f_calving.attrs['repartitionFile'] = myargs.repartition_file 422 f_calving.attrs['repartitionVar'] = repartitionVar 423 f_calving.attrs['gridsFile'] = grids 424 f_calving.attrs['areasFile'] = areas 425 f_calving.attrs['masksFile'] = masks 426 f_calving.attrs['timeStamp'] = time.asctime() 427 f_calving.attrs['HOSTNAME'] = platform.node() 428 f_calving.attrs['LOGNAME'] = os.getlogin() 429 f_calving.attrs['Python'] = "Python version " + platform.python_version() 430 f_calving.attrs['OS'] = platform.system() 431 f_calving.attrs['release'] = platform.release() 432 f_calving.attrs['hardware'] = platform.machine() 433 f_calving.attrs['conventions'] = "SCRIP" 434 if src_name == 'lmd' : f_calving.attrs['source_grid'] = "curvilinear" 435 if src_name == 'ico' : f_calving.attrs['source_grid'] = "unstructured" 436 f_calving.attrs['dest_grid'] = "curvilinear" 437 f_calving.attrs['Model'] = "IPSL CM6" 438 f_calving.attrs['SVN_Author'] = "$Author$" 439 f_calving.attrs['SVN_Date'] = "$Date$" 440 f_calving.attrs['SVN_Revision'] = "$Revision$" 441 f_calving.attrs['SVN_Id'] = "$Id$" 442 f_calving.attrs['SVN_HeadURL'] = "$HeadURL$" 443 444 ## 445 f_calving.to_netcdf ( calving, mode='w', format=FmtNetcdf ) 446 f_calving.close() 447 448 print ('That''s all folks !') 449 ## ====================================================================================== -
TOOLS/MOSAIX/CreateWeightsMask.bash
r6064 r6066 69 69 # ================================= 70 70 #CplModel=ORCA2.3xLMD9695 71 CplModel=ORCA2.3xICO3072 #CplModel=ORCA2.3xICO4071 #CplModel=ORCA2.3xICO30 72 CplModel=ORCA2.3xICO40 73 73 #CplModel=eORCA1.2xLMD144142 74 74 #CplModel=eORCA1.2xLMD256256 … … 800 800 --atmQuantity=${runOff_atmQuantity} --oceQuantity=${runOff_oceQuantity} --ocePerio=${OcePerio} 801 801 fi 802 exit 802 803 803 ## 804 804 echo ${Titre}"Calving weights"${Norm} … … 809 809 cp /ccc/work/cont003/gencmip6/deshayej/eORCA_R025_runoff_v1.2.nc . 810 810 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_nosouth.nc --fmt=${FMT_XIOS} \ 811 --oce=${OCE} --atm=${ATM} --type=nosouth --dir=. 811 --oce=${OCE} --atm=${ATM} --type=nosouth \ 812 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 813 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 812 814 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceberg.nc --fmt=${FMT_XIOS} \ 813 --oce=${OCE} --atm=${ATM} --type=iceberg --dir=. --repartition_file=eORCA_R025_runoff_v1.2.nc --repartition_var=Icb_flux 815 --oce=${OCE} --atm=${ATM} --type=iceberg --repartition_file=eORCA_R025_runoff_v1.2.nc --repartition_var=Icb_flux \ 816 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 817 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 814 818 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceshelf.nc --fmt=${FMT_XIOS} \ 815 --oce=${OCE} --atm=${ATM} --type=iceshelf --dir=. --repartition_file=eORCA_R025_runoff_v1.2.nc --repartition_var=sornfisf 819 --oce=${OCE} --atm=${ATM} --type=iceshelf --repartition_file=eORCA_R025_runoff_v1.2.nc --repartition_var=sornfisf \ 820 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 821 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 816 822 ;; 817 823 818 ( eORCA1.2 )824 ( eORCA1.2 | eORCA1.4 ) 819 825 cp ${R_IN}/OCE/NEMO/ORCA1_LIM3_PISCES/v3.6_stable/runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc . 820 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_nosouth.nc --fmt=${FMT_XIOS}\ 821 --oce=${OCE} --atm=${ATM} --type=nosouth --dir=. 822 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceberg.nc --fmt=${FMT_XIOS}\ 823 --oce=${OCE} --atm=${ATM} --type=iceberg --dir=. --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=Icb_flux 824 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceshelf.nc --fmt=${FMT_XIOS}\ 825 --oce=${OCE} --atm=${ATM} --type=iceshelf --dir=. --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=sornfisf 826 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_nosouth.nc --fmt=${FMT_XIOS} \ 827 --oce=${OCE} --atm=${ATM} --type=nosouth \ 828 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 829 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 830 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceberg.nc --fmt=${FMT_XIOS} \ 831 --oce=${OCE} --atm=${ATM} --type=iceberg --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=Icb_flux \ 832 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 833 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 834 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_iceshelf.nc --fmt=${FMT_XIOS }\ 835 --oce=${OCE} --atm=${ATM} --type=iceshelf --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=sornfisf \ 836 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 837 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 826 838 ;; 827 839 828 840 ( * ) 829 841 ${PyRun} python3 -u CalvingWeights.py --output=rmp_t${atm}_to_t${oce}_calving_full.nc --fmt=${FMT_XIOS} \ 830 --oce=${OCE} --atm=${ATM} --type=full --dir=. --ocePerio=${OcePerio} 842 --oce=${OCE} --atm=${ATM} --type=full --ocePerio=${OcePerio} \ 843 --grids=grids_${CplModel}.nc --areas=areas_${CplModel}.nc --masks=masks_${CplModel}.nc \ 844 --o2a=${ATM}_grid_maskFrom_${OCE}.nc 831 845 ;; 832 846 -
TOOLS/MOSAIX/RunoffWeights.py
r6064 r6066 26 26 27 27 ## Modules 28 import netCDF4 29 import numpy as np 28 import numpy as np, xarray as xr 30 29 import nemo 31 30 from scipy import ndimage … … 55 54 56 55 # Adding arguments 57 parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025', 'eORCA025.1'] ) 56 parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', 57 choices=['ORCA2.3', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) 58 58 parser.add_argument ('--atm' , help='atm model name', type=str, default='LMD9695' ) 59 59 parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', type=int, default=1 ) 60 60 parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)' , type=int, default=2 ) 61 parser.add_argument ('--atmQuantity' , help='Quantity if atm provides quantities (m/s), Surfacic if atm provided flux (m/s/m2)' , type=str, default='Quantity', choices=['Quantity', 'Surfacic'] ) 62 parser.add_argument ('--oceQuantity' , help='Quantity if oce requires quantities (ks/s), Surfacic if oce requires flux (m/s/m2)' , type=str, default='Surfacic', choices=['Quantity', 'Surfacic'] ) 61 parser.add_argument ('--atmQuantity' , help='Quantity if atm provides quantities (m/s), Surfacic if atm provided flux (m/s/m2)' , type=str, 62 default='Quantity', choices=['Quantity', 'Surfacic'] ) 63 parser.add_argument ('--oceQuantity' , help='Quantity if oce requires quantities (ks/s), Surfacic if oce requires flux (m/s/m2)', type=str, 64 default='Surfacic', choices=['Quantity', 'Surfacic'] ) 63 65 parser.add_argument ('--searchRadius' , help='max distance to connect a land point to an ocean point (in km)', type=float, default=600.0 ) 64 66 parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) … … 67 69 parser.add_argument ('--o2a' , help='o2a file name' , default='o2a.nc' ) 68 70 parser.add_argument ('--output', help='output rmp file name', default='rmp_tlmd_to_torc_runoff.nc' ) 69 parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 71 parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', 72 choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 70 73 parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=int, default=0 ) 71 74 … … 104 107 105 108 # Ocean grid periodicity 106 oce_perio =myargs.ocePerio109 oce_perio = myargs.ocePerio 107 110 108 111 ### Read coordinates of all models 109 112 ### 110 113 111 diaFile = netCDF4.Dataset ( o2a )112 gridFile = netCDF4.Dataset ( grids )113 areaFile = netCDF4.Dataset ( areas )114 maskFile = netCDF4.Dataset ( masks )115 116 o2aFrac = diaFile ['OceFrac'] [:].squeeze()114 diaFile = xr.open_dataset ( o2a ) 115 gridFile = xr.open_dataset ( grids ) 116 areaFile = xr.open_dataset ( areas ) 117 maskFile = xr.open_dataset ( masks ) 118 119 o2aFrac = diaFile ['OceFrac'].squeeze() 117 120 o2aFrac = np.where ( np.abs(o2aFrac) < 1E10, o2aFrac, 0.0) 118 121 … … 121 124 atm_grid_rank = len(gridFile['t'+atm_n+'.lat'][:].shape) 122 125 123 atm_grid_center_lat = gridFile['t'+atm_n+'.lat'] [:].ravel()124 atm_grid_center_lon = gridFile['t'+atm_n+'.lon'] [:].ravel()125 atm_grid_corner_lat = np.reshape ( gridFile['t'+atm_n+'.cla'][:], (atm_nvertex, atm_grid_size))126 atm_grid_corner_lon = np.reshape ( gridFile['t'+atm_n+'.clo'][:], (atm_nvertex, atm_grid_size))127 atm_grid_area = areaFile['t'+atm_n+'.srf'] [:].ravel()128 atm_grid_imask = 1-maskFile['t'+atm_n+'.msk'][:].squeeze() .ravel()126 atm_grid_center_lat = gridFile['t'+atm_n+'.lat'].squeeze() 127 atm_grid_center_lon = gridFile['t'+atm_n+'.lon'].squeeze() 128 atm_grid_corner_lat = gridFile['t'+atm_n+'.cla'].squeeze() 129 atm_grid_corner_lon = gridFile['t'+atm_n+'.clo'].squeeze() 130 atm_grid_area = areaFile['t'+atm_n+'.srf'].squeeze() 131 atm_grid_imask = 1-maskFile['t'+atm_n+'.msk'][:].squeeze() 129 132 atm_grid_dims = gridFile['t'+atm_n+'.lat'][:].shape 133 134 if atmDomainType == 'unstructured' : 135 atm_grid_center_lat = atm_grid_center_lat.rename ({'ycell':'cell'}) 136 atm_grid_center_lon = atm_grid_center_lon.rename ({'ycell':'cell'}) 137 atm_grid_corner_lat = atm_grid_corner_lat.rename ({'ycell':'cell'}) 138 atm_grid_corner_lon = atm_grid_corner_lon.rename ({'ycell':'cell'}) 139 atm_grid_area = atm_grid_area.rename ({'ycell':'cell'}) 140 atm_grid_imask = atm_grid_imask.rename ({'ycell':'cell'}) 141 142 if atmDomainType == 'rectilinear' : 143 atm_grid_center_lat = atm_grid_center_lat.stack (cell=['y', 'x']) 144 atm_grid_center_lon = atm_grid_center_lon.stack (cell=['y', 'x']) 145 atm_grid_corner_lat = atm_grid_corner_lat.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'}) 146 atm_grid_corner_lon = atm_grid_corner_lon.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'}) 147 atm_grid_area = atm_grid_area.stack (cell=['y', 'x']) 148 atm_grid_imask = atm_grid_imask.stack (cell=['y', 'x']) 130 149 131 150 atm_perio = 0 … … 137 156 oce_grid_rank = len(gridFile['torc.lat'][:].shape) 138 157 139 oce_grid_center_lat = gridFile['torc.lat'] [:].ravel()140 oce_grid_center_lon = gridFile['torc.lon'] [:].ravel()141 oce_grid_corner_lat = np.reshape( gridFile['torc.cla'][:], (oce_nvertex, oce_grid_size))142 oce_grid_corner_lon = np.reshape( gridFile['torc.clo'][:], (oce_nvertex, oce_grid_size))143 oce_grid_area = areaFile['torc.srf'] [:].ravel()144 oce_grid_imask = 1-maskFile['torc.msk'] [:].ravel()158 oce_grid_center_lat = gridFile['torc.lat'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 159 oce_grid_center_lon = gridFile['torc.lon'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 160 oce_grid_corner_lat = gridFile['torc.cla'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 161 oce_grid_corner_lon = gridFile['torc.clo'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 162 oce_grid_area = areaFile['torc.srf'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 163 oce_grid_imask = 1-maskFile['torc.msk'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 145 164 oce_grid_dims = gridFile['torc.lat'][:].shape 165 146 166 if oce_perio == 0 : 147 167 if oce_jpi == 182 : oce_perio = 4 # ORCA 2 148 168 if oce_jpi == 362 : oce_perio = 6 # ORCA 1 149 169 if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 170 150 171 print ("Oce NPERIO parameter : {:d}".format(oce_perio)) 151 oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask , (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel()172 oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel() 152 173 oce_address = np.arange(oce_jpj*oce_jpi) 153 174 154 175 print ("Fill closed sea with image processing library") 155 176 oce_grid_imask2D = np.reshape(oce_grid_pmask,(oce_jpj,oce_jpi)) 156 oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), nperio=oce_perio, cd_type='T', sval=0 ) 177 oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), 178 nperio=oce_perio, cd_type='T', sval=0 ) 157 179 oce_grid_imask = oce_grid_imask2D.ravel() 158 180 ## 159 181 print ("Computes an ocean coastal band") 160 182 161 oceLand2D = np.reshape ( np.where (oce_grid_pmask [:]< 0.5, True, False), (oce_jpj, oce_jpi) )162 oceOcean2D = np.reshape ( np.where (oce_grid_pmask [:]> 0.5, True, False), (oce_jpj, oce_jpi) )183 oceLand2D = np.reshape ( np.where (oce_grid_pmask < 0.5, True, False), (oce_jpj, oce_jpi) ) 184 oceOcean2D = np.reshape ( np.where (oce_grid_pmask > 0.5, True, False), (oce_jpj, oce_jpi) ) 163 185 164 186 NNocean = 1+2*oceCoastWidth … … 168 190 169 191 oceOceanFiltered = oceOceanFiltered2D.ravel() 170 oceLand = oceLand2D.ravel ()192 oceLand = oceLand2D.ravel () 171 193 oceOcean = oceOcean2D.ravel() 172 194 oceCoast = oceCoast2D.ravel() … … 232 254 num_links = int(z_mask.sum()) 233 255 if num_links == 0 : continue 234 z_area = oceCoast_grid_area[z_mask].sum() 256 z_area = oceCoast_grid_area[z_mask].sum().values 235 257 poids = np.ones ((num_links),dtype=np.float64) / z_area 236 258 if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja] 237 259 if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] 238 260 if ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print 239 print ( 'ja:{:8d}, num_links:{:8d}, z_area:{:8.4e}, atm area:{:8.4e}, weights sum:{:8.4e} '.format(ja, num_links, z_area, atm_grid_area[ja], poids.sum() ) ) 261 print ( 'ja:{:8d}, num_links:{:8d}, z_area:{:8.4e}, atm area:{:8.4e}, weights sum:{:8.4e} ' 262 .format(ja, num_links, z_area, atm_grid_area[ja].values, poids.sum() ) ) 240 263 # 241 264 matrix_local = poids … … 254 277 print ("Write output file") 255 278 runoff = myargs.output 256 f_runoff = netCDF4.Dataset ( runoff, 'w', format=FmtNetcdf )257 279 print ('Output file: ' + runoff ) 258 280 259 f_runoff.Conventions = "CF-1.6" 260 f_runoff.source = "IPSL Earth system model" 261 f_runoff.group = "ICMC IPSL Climate Modelling Center" 262 f_runoff.Institution = "IPSL https.//www.ipsl.fr" 263 f_runoff.Ocean = oce_Name + " https://www.nemo-ocean.eu" 264 f_runoff.Atmosphere = atm_Name + " http://lmdz.lmd.jussieu.fr" 265 f_runoff.associatedFiles = grids + " " + areas + " " + masks 266 f_runoff.directory = os.getcwd () 267 f_runoff.description = "Generated with RunoffWeights.py" 268 f_runoff.title = runoff 269 f_runoff.Program = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) 270 f_runoff.atmCoastWidth = "{:d} grid points".format(atmCoastWidth) 271 f_runoff.oceCoastWidth = "{:d} grid points".format(oceCoastWidth) 272 f_runoff.searchRadius = "{:.0f} km".format(searchRadius/1000.) 273 f_runoff.atmQuantity = myargs.atmQuantity 274 f_runoff.oceQuantity = myargs.oceQuantity 275 f_runoff.gridsFile = grids 276 f_runoff.areasFile = areas 277 f_runoff.masksFile = masks 278 f_runoff.o2aFile = o2a 279 f_runoff.timeStamp = time.asctime() 280 f_runoff.HOSTNAME = platform.node() 281 #f_runoff.LOGNAME = os.getlogin() 282 f_runoff.Python = "Python version " + platform.python_version() 283 f_runoff.OS = platform.system() 284 f_runoff.release = platform.release() 285 f_runoff.hardware = platform.machine() 286 f_runoff.conventions = "SCRIP" 287 f_runoff.source_grid = "curvilinear" 288 f_runoff.dest_grid = "curvilinear" 289 f_runoff.Model = "IPSL CM6" 290 f_runoff.SVN_Author = "$Author$" 291 f_runoff.SVN_Date = "$Date$" 292 f_runoff.SVN_Revision = "$Revision$" 293 f_runoff.SVN_Id = "$Id$" 294 f_runoff.SVN_HeadURL = "$HeadURL$" 295 296 d_num_links = f_runoff.createDimension ('num_links' , num_links ) 297 d_num_wgts = f_runoff.createDimension ('num_wgts' , 1 ) 298 299 d_atm_grid_size = f_runoff.createDimension ('src_grid_size' , atm_grid_size ) 300 d_atm_grid_corners = f_runoff.createDimension ('src_grid_corners', atm_grid_corner_lon.shape[0] ) 301 d_atm_grid_rank = f_runoff.createDimension ('src_grid_rank' , atm_grid_rank ) 302 303 d_oce_grid_size = f_runoff.createDimension ('dst_grid_size' , oce_grid_size ) 304 d_oce_grid_corners = f_runoff.createDimension ('dst_grid_corners', oce_grid_corner_lon.shape[0] ) 305 d_oce_grid_rank = f_runoff.createDimension ('dst_grid_rank' , oce_grid_rank ) 306 307 v_remap_matrix = f_runoff.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') ) 308 309 v_atm_address = f_runoff.createVariable ( 'src_address' , 'i4', ('num_links',) ) 310 v_atm_address.convention = "Fortran style addressing, starting at 1" 311 v_oce_address = f_runoff.createVariable ( 'dst_address' , 'i4', ('num_links',) ) 312 v_oce_address.convention = "Fortran style addressing, starting at 1" 313 314 v_remap_matrix[:] = remap_matrix 315 v_atm_address [:] = atm_address + 1 # OASIS uses Fortran style indexing, starting at one 316 v_oce_address [:] = oce_address + 1 317 318 v_atm_grid_dims = f_runoff.createVariable ( 'src_grid_dims' , 'i4', ('src_grid_rank',) ) 319 v_atm_grid_center_lon = f_runoff.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) ) 320 v_atm_grid_center_lat = f_runoff.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) ) 321 v_atm_grid_center_lon.units='degrees_east' ; v_atm_grid_center_lon.long_name='Longitude' ; v_atm_grid_center_lon.long_name='longitude' ; v_atm_grid_center_lon.bounds="src_grid_corner_lon" 322 v_atm_grid_center_lat.units='degrees_north' ; v_atm_grid_center_lat.long_name='Latitude' ; v_atm_grid_center_lat.long_name='latitude ' ; v_atm_grid_center_lat.bounds="src_grid_corner_lat" 323 v_atm_grid_corner_lon = f_runoff.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners') ) 324 v_atm_grid_corner_lat = f_runoff.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners') ) 325 v_atm_grid_corner_lon.units="degrees_east" 326 v_atm_grid_corner_lat.units="degrees_north" 327 v_atm_grid_area = f_runoff.createVariable ( 'src_grid_area' , 'f8', ('src_grid_size',) ) 328 v_atm_grid_area.long_name="Grid area" ; v_atm_grid_area.standard_name="cell_area" ; v_atm_grid_area.units="m2" 329 v_atm_grid_imask = f_runoff.createVariable ( 'src_grid_imask' , 'i4', ('src_grid_size',) ) 330 v_atm_grid_imask.long_name="Land-sea mask" ; v_atm_grid_imask.units="Land:1, Ocean:0" 331 v_atm_grid_pmask = f_runoff.createVariable ( 'src_grid_pmask' , 'i4', ('src_grid_size',) ) 332 v_atm_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_atm_grid_pmask.units="Land:1, Ocean:0" 333 334 v_atm_grid_dims [:] = atm_grid_dims 335 v_atm_grid_center_lon[:] = atm_grid_center_lon[:] 336 v_atm_grid_center_lat[:] = atm_grid_center_lat[:] 337 v_atm_grid_corner_lon[:] = np.transpose(atm_grid_corner_lon) 338 v_atm_grid_corner_lat[:] = np.transpose(atm_grid_corner_lat) 339 v_atm_grid_area [:] = atm_grid_area[:] 340 v_atm_grid_imask [:] = atm_grid_imask[:] 341 v_atm_grid_pmask [:] = atm_grid_pmask[:] 281 282 remap_matrix = xr.DataArray ( np.reshape(remap_matrix, (num_links, 1)), dims = ['num_links', 'num_wgts'] ) 283 284 # OASIS uses Fortran style indexing, starting at one 285 src_address = xr.DataArray ( atm_address.astype(np.int32)+1, dims = ['num_links'], 286 attrs={"convention": "Fortran style addressing, starting at 1"}) 287 dst_address = xr.DataArray ( oce_address.astype(np.int32)+1, dims = ['num_links'], 288 attrs={"convention": "Fortran style addressing, starting at 1"}) 289 290 src_grid_dims = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32), dims = ['src_grid_rank',] ) 291 src_grid_center_lon = xr.DataArray (atm_grid_center_lon.values , dims = ['src_grid_size',] ) 292 src_grid_center_lat = xr.DataArray (atm_grid_center_lat.values , dims = ['src_grid_size',] ) 293 294 src_grid_center_lon.attrs['units']='degrees_east' ; src_grid_center_lon.attrs['long_name']='Longitude' 295 src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" 296 src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' 297 src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" 298 299 src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) 300 src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) 301 src_grid_corner_lon.attrs['units']="degrees_east" 302 src_grid_corner_lat.attrs['units']="degrees_north" 303 304 src_grid_area = xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] ) 305 src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" 306 307 src_grid_imask = xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] ) 308 src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" 309 310 src_grid_pmask = xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] ) 311 src_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; src_grid_pmask.attrs['units']="Land:1, Ocean:0" 342 312 343 313 # -- 344 345 v_oce_grid_dims = f_runoff.createVariable ( 'dst_grid_dims' , 'i4', ('dst_grid_rank',) ) 346 v_oce_grid_center_lon = f_runoff.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) ) 347 v_oce_grid_center_lat = f_runoff.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) ) 348 v_oce_grid_center_lon.units='degrees_east' ; v_oce_grid_center_lon.long_name='Longitude' ; v_oce_grid_center_lon.long_name='longitude' ; v_oce_grid_center_lon.bounds="dst_grid_corner_lon" 349 v_oce_grid_center_lat.units='degrees_north' ; v_oce_grid_center_lat.long_name='Latitude' ; v_oce_grid_center_lat.long_name='latitude' ; v_oce_grid_center_lat.bounds="dst_grid_corner_lat" 350 v_oce_grid_corner_lon = f_runoff.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners') ) 351 v_oce_grid_corner_lat = f_runoff.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners') ) 352 v_oce_grid_corner_lon.units="degrees_east" 353 v_oce_grid_corner_lat.units="degrees_north" 354 v_oce_grid_area = f_runoff.createVariable ( 'dst_grid_area' , 'f8', ('dst_grid_size',) ) 355 v_oce_grid_area.long_name="Grid area" ; v_oce_grid_area.standard_name="cell_area" ; v_oce_grid_area.units="m2" 356 v_oce_grid_imask = f_runoff.createVariable ( 'dst_grid_imask' , 'i4', ('dst_grid_size',) ) 357 v_oce_grid_imask.long_name="Land-sea mask" ; v_oce_grid_imask.units="Land:1, Ocean:0" 358 v_oce_grid_pmask = f_runoff.createVariable ( 'dst_grid_pmask' , 'i4', ('dst_grid_size',) ) 359 v_oce_grid_pmask.long_name="Land-sea mask (periodicity removed)" ; v_oce_grid_pmask.units="Land:1, Ocean:0" 360 361 v_oce_grid_dims [:] = oce_grid_dims 362 v_oce_grid_center_lon[:] = oce_grid_center_lon[:] 363 v_oce_grid_center_lat[:] = oce_grid_center_lat[:] 364 v_oce_grid_corner_lon[:] = np.transpose(oce_grid_corner_lon) 365 v_oce_grid_corner_lat[:] = np.transpose(oce_grid_corner_lat) 366 v_oce_grid_area [:] = oce_grid_area[:] 367 v_oce_grid_imask [:] = oce_grid_imask[:] 368 v_oce_grid_pmask [:] = oce_grid_pmask[:] 369 370 v_atm_lon_addressed = f_runoff.createVariable ( 'src_lon_addressed' , 'f8', ('num_links',) ) 371 v_atm_lat_addressed = f_runoff.createVariable ( 'src_lat_addressed' , 'f8', ('num_links',) ) 372 v_atm_area_addressed = f_runoff.createVariable ( 'src_area_addressed' , 'f8', ('num_links',) ) 373 v_atm_imask_addressed = f_runoff.createVariable ( 'src_imask_addressed', 'i4', ('num_links',) ) 374 v_atm_pmask_addressed = f_runoff.createVariable ( 'src_pmask_addressed', 'i4', ('num_links',) ) 375 376 v_oce_lon_addressed = f_runoff.createVariable ( 'dst_lon_addressed' , 'f8', ('num_links',) ) 377 v_oce_lat_addressed = f_runoff.createVariable ( 'dst_lat_addressed' , 'f8', ('num_links',) ) 378 v_oce_area_addressed = f_runoff.createVariable ( 'dst_area_addressed' , 'f8', ('num_links',) ) 379 v_oce_imask_addressed = f_runoff.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) ) 380 v_oce_pmask_addressed = f_runoff.createVariable ( 'dst_pmask_addressed', 'i4', ('num_links',) ) 381 382 v_atm_lon_addressed.long_name="Longitude" ; v_atm_lon_addressed.standard_name="longitude" ; v_atm_lon_addressed.units="degrees_east" 383 v_atm_lat_addressed.long_name="Latitude" ; v_atm_lat_addressed.standard_name="latitude" ; v_atm_lat_addressed.units="degrees_north" 384 v_atm_lon_addressed [:] = atm_grid_center_lon[atm_address] 385 v_atm_lat_addressed [:] = atm_grid_center_lat[atm_address] 386 v_atm_area_addressed [:] = atm_grid_area[atm_address] 387 v_atm_imask_addressed[:] = 1-atm_grid_imask[atm_address] 388 v_atm_pmask_addressed[:] = 1-atm_grid_pmask[atm_address] 389 390 v_oce_lon_addressed.long_name="Longitude" ; v_oce_lon_addressed.standard_name="longitude" ; v_oce_lon_addressed.units="degrees_east" 391 v_oce_lat_addressed.long_name="Latitude" ; v_oce_lat_addressed.standard_name="latitude" ; v_oce_lat_addressed.units="degrees_north" 392 v_oce_lon_addressed [:] = oce_grid_center_lon[oce_address] 393 v_oce_lat_addressed [:] = oce_grid_center_lat[oce_address]#.ravel() 394 v_oce_area_addressed [:] = oce_grid_area[oce_address] 395 v_oce_imask_addressed[:] = 1-oce_grid_imask[oce_address] 396 v_oce_pmask_addressed[:] = 1-oce_grid_pmask[oce_address] 314 dst_grid_dims = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32), dims = ['dst_grid_rank',] ) 315 dst_grid_center_lon = xr.DataArray (oce_grid_center_lon.values, dims = ['dst_grid_size',] ) 316 dst_grid_center_lat = xr.DataArray (oce_grid_center_lat.values, dims = ['dst_grid_size',] ) 317 318 dst_grid_center_lon.attrs['units']='degrees_east' ; dst_grid_center_lon.attrs['long_name']='Longitude' 319 dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" 320 dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' 321 dst_grid_center_lat.attrs['long_name']='latitude ' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" 322 323 dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 324 dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 325 dst_grid_corner_lon.attrs['units']="degrees_east" 326 dst_grid_corner_lat.attrs['units']="degrees_north" 327 328 dst_grid_area = xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] ) 329 dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 330 331 dst_grid_imask = xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] ) 332 dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" 333 334 dst_grid_pmask = xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] ) 335 dst_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; dst_grid_pmask.attrs['units']="Land:1, Ocean:0" 336 337 src_lon_addressed = xr.DataArray (atm_grid_center_lon.values[atm_address] , dims = ['num_links'] ) 338 src_lat_addressed = xr.DataArray (atm_grid_center_lat.values[atm_address] , dims = ['num_links'] ) 339 src_area_addressed = xr.DataArray (atm_grid_area .values[atm_address] , dims = ['num_links'] ) 340 src_imask_addressed = xr.DataArray (1-atm_grid_imask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) 341 src_pmask_addressed = xr.DataArray (1-atm_grid_pmask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) 342 343 dst_lon_addressed = xr.DataArray (oce_grid_center_lon.values[atm_address], dims = ['num_links'] ) 344 dst_lat_addressed = xr.DataArray (oce_grid_center_lat.values[oce_address], dims = ['num_links'] ) 345 dst_area_addressed = xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32) , dims = ['num_links'] ) 346 dst_imask_addressed = xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32) , dims = ['num_links'] ) 347 dst_pmask_addressed = xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32) , dims = ['num_links'] ) 348 349 src_lon_addressed.attrs['long_name']="Longitude" ; src_lon_addressed.attrs['standard_name']="longitude" ; src_lon_addressed.attrs['units']="degrees_east" 350 src_lat_addressed.attrs['long_name']="Latitude" ; src_lat_addressed.attrs['standard_name']="latitude" ; src_lat_addressed.attrs['units']="degrees_north" 351 352 dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" 353 dst_lat_addressed.attrs['long_name']="Latitude" ; dst_lat_addressed.attrs['standard_name']="latitude" ; dst_lat_addressed.attrs['units']="degrees_north" 397 354 398 355 if atmDomainType == 'rectilinear' : 399 v_atmLand = f_runoff.createVariable ( 'atmLand' , 'i4', ('src_grid_size',) ) 400 v_atmLandFiltered = f_runoff.createVariable ( 'atmLandFiltered', 'f4', ('src_grid_size',) ) 401 v_atmLandFrac = f_runoff.createVariable ( 'atmLandFrac' , 'i4', ('src_grid_size',) ) 402 v_atmFrac = f_runoff.createVariable ( 'atmFrac' , 'i4', ('src_grid_size',) ) 403 v_atmOcean = f_runoff.createVariable ( 'atmOcean' , 'i4', ('src_grid_size',) ) 404 v_atmOceanFrac = f_runoff.createVariable ( 'atmOceanFrac' , 'i4', ('src_grid_size',) ) 405 406 v_atmLand[:] = atmLand.ravel() 407 v_atmLandFrac[:] = atmLandFrac.ravel() 408 v_atmLandFiltered[:] = atmLandFiltered.ravel() 409 v_atmFrac[:] = atmFrac.ravel() 410 v_atmOcean[:] = atmOcean.ravel() 411 v_atmOceanFrac[:] = atmOceanFrac.ravel() 412 413 v_atmCoast = f_runoff.createVariable ( 'atmCoast' , 'i4', ('src_grid_size',) ) 414 v_atmCoast[:] = atmCoast 415 416 v_oceLand = f_runoff.createVariable ( 'oceLand' , 'i4', ('dst_grid_size',) ) 417 v_oceOcean = f_runoff.createVariable ( 'oceOcean' , 'i4', ('dst_grid_size',) ) 418 v_oceOceanFiltered = f_runoff.createVariable ( 'oceOceanFiltered', 'f4', ('dst_grid_size',) ) 419 v_oceCoast = f_runoff.createVariable ( 'oceCoast' , 'i4', ('dst_grid_size',) ) 420 421 v_oceLand[:] = oceLand 422 v_oceOcean[:] = oceOcean 423 v_oceOceanFiltered[:] = oceOceanFiltered 424 v_oceCoast[:] = oceCoast 425 426 f_runoff.sync () 427 428 ## 356 atmLand = xr.DataArray ( atmLand.ravel() , dims = ['src_grid_size',] ) 357 atmLandFiltered = xr.DataArray ( atmLandFrac.ravel() , dims = ['src_grid_size',] ) 358 atmLandFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) 359 atmFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) 360 atmOcean = xr.DataArray ( atmOcean.ravel() , dims = ['src_grid_size',] ) 361 atmOceanFrac = xr.DataArray ( atmOceanFrac.ravel(), dims = ['src_grid_size',] ) 362 363 atmCoast = xr.DataArray (atmCoast.astype(np.int32) , dims = ['src_grid_size',]) 364 oceLand = xr.DataArray (oceLand.astype(np.int32) , dims = ['dst_grid_size',]) 365 oceOcean = xr.DataArray (oceOcean.astype(np.int32) , dims = ['dst_grid_size',]) 366 oceOceanFiltered = xr.DataArray (oceOceanFiltered.astype(np.float32), dims = ['dst_grid_size',]) 367 oceCoast = xr.DataArray (oceCoast.astype(np.int32) , dims = ['dst_grid_size',]) 368 369 370 f_runoff = xr.Dataset ( { 371 'src_address' : src_address, 372 'dst_address' : dst_address, 373 'src_grid_dims' : src_grid_dims, 374 'src_grid_center_lon' : src_grid_center_lon, 375 'src_grid_center_lat' : src_grid_center_lat, 376 'src_grid_corner_lon' : src_grid_corner_lon, 377 'src_grid_corner_lat' : src_grid_corner_lat, 378 'src_grid_area' : src_grid_area, 379 'src_grid_area' : src_grid_area, 380 'src_grid_pmask' : src_grid_pmask, 381 'dst_grid_dims' : dst_grid_dims, 382 'dst_grid_center_lon' : dst_grid_center_lon, 383 'st_grid_center_lat' : dst_grid_center_lat, 384 'dst_grid_corner_lon' : dst_grid_corner_lon, 385 'dst_grid_corner_lat' : dst_grid_corner_lat, 386 'dst_grid_area' : dst_grid_area, 387 'dst_grid_imask' : dst_grid_imask, 388 'dst_grid_pmask' : dst_grid_pmask, 389 'src_lon_addressed' : src_lon_addressed, 390 'src_lat_addressed' : src_lat_addressed, 391 'src_area_addressed' : src_area_addressed, 392 'dst_lon_addressed' : dst_lon_addressed, 393 'dst_lat_addressed' : dst_lat_addressed, 394 'dst_area_addressed' : dst_area_addressed, 395 'dst_imask_addressed' : dst_imask_addressed, 396 'dst_pmask_addressed' : dst_pmask_addressed, 397 'atmCoast' : atmCoast, 398 'oceLand' : oceLand, 399 'oceOcean' : oceOcean, 400 'oceOceanFiltered' : oceOceanFiltered, 401 'oceCoast' : oceCoast 402 } ) 403 404 f_runoff.attrs['Conventions'] = "CF-1.6" 405 f_runoff.attrs['source'] = "IPSL Earth system model" 406 f_runoff.attrs['group'] = "ICMC IPSL Climate Modelling Center" 407 f_runoff.attrs['Institution'] = "IPSL https.//www.ipsl.fr" 408 f_runoff.attrs['Ocean'] = oce_Name + " https://www.nemo-ocean.eu" 409 f_runoff.attrs['Atmosphere'] = atm_Name + " http://lmdz.lmd.jussieu.fr" 410 f_runoff.attrs['associatedFiles'] = grids + " " + areas + " " + masks 411 f_runoff.attrs['directory'] = os.getcwd () 412 f_runoff.attrs['description'] = "Generated with RunoffWeights.py" 413 f_runoff.attrs['title'] = runoff 414 args = '' 415 for i in np.arange(1,len(sys.argv[1:])) : args += sys.argv[i] + ' ' 416 f_runoff.attrs['Program'] = "Generated by " + sys.argv[0] + " with flags " + args 417 f_runoff.attrs['atmCoastWidth'] = "{:d} grid points".format(atmCoastWidth) 418 f_runoff.attrs['oceCoastWidth'] = "{:d} grid points".format(oceCoastWidth) 419 f_runoff.attrs['searchRadius'] = "{:.0f} km".format(searchRadius/1000.) 420 f_runoff.attrs['atmQuantity'] = myargs.atmQuantity 421 f_runoff.attrs['oceQuantity'] = myargs.oceQuantity 422 f_runoff.attrs['gridsFile'] = grids 423 f_runoff.attrs['areasFile'] = areas 424 f_runoff.attrs['masksFile'] = masks 425 f_runoff.attrs['o2aFile'] = o2a 426 f_runoff.attrs['timeStamp'] = time.asctime() 427 f_runoff.attrs['HOSTNAME'] = platform.node() 428 f_runoff.attrs['LOGNAME'] = os.getlogin() 429 f_runoff.attrs['Python'] = "Python version " + platform.python_version() 430 f_runoff.attrs['OS'] = platform.system() 431 f_runoff.attrs['release'] = platform.release() 432 f_runoff.attrs['hardware'] = platform.machine() 433 f_runoff.attrs['conventions'] = "SCRIP" 434 f_runoff.attrs['source_grid'] = "curvilinear" 435 f_runoff.attrs['dest_grid'] = "curvilinear" 436 f_runoff.attrs['Model'] = "IPSL CM6" 437 f_runoff.attrs['SVN_Author'] = "$Author$" 438 f_runoff.attrs['SVN_Date'] = "$Date$" 439 f_runoff.attrs['SVN_Revision'] = "$Revision$" 440 f_runoff.attrs['SVN_Id'] = "$Id$" 441 f_runoff.attrs['SVN_HeadURL'] = "$HeadURL$" 442 443 f_runoff.to_netcdf ( runoff, mode='w', format=FmtNetcdf ) 429 444 f_runoff.close () 430 445 446 ## 447 431 448 print ('That''s all folks !') 432 433 449 ## ======================================================================================
Note: See TracChangeset
for help on using the changeset viewer.