### =========================================================================== ### ### Compute calving weights. ### ### =========================================================================== ## ## Warning, to install, configure, run, use any of Olivier Marti's ## software or to read the associated documentation you'll need at least ## one (1) brain in a reasonably working order. Lack of this implement ## will void any warranties (either express or implied). ## O. Marti assumes no responsability for errors, omissions, ## data loss, or any other consequences caused directly or indirectly by ## the usage of his software by incorrectly or partially configured ## personal. ## import netCDF4, numpy as np import sys, os, platform, argparse, textwrap, time from scipy import ndimage import nemo ## SVN information __Author__ = "$Author$" __Date__ = "$Date$" __Revision__ = "$Revision$" __Id__ = "$Id$" __HeadURL__ = "$HeadURL$" __SVN_Date__ = "$SVN_Date: $" ### ### ===== Handling command line parameters ================================================== # Creating a parser parser = argparse.ArgumentParser ( description = """Compute calving weights""", epilog='-------- End of the help message --------') # Adding arguments parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025', 'eORCA025.1'] ) parser.add_argument ('--atm' , help='atm model name (ICO* or LMD*)', type=str, default='ICO40' ) parser.add_argument ('--type' , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full' ) parser.add_argument ('--dir' , help='Directory of input file', type=str, default='./' ) parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' ) parser.add_argument ('--repartition_var' , help='Variable name for iceshelf' , type=str, default=None) parser.add_argument ('--output' , help='output rmp file name', default='rmp_tlmd_to_torc_calving.nc' ) parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=int, default=0 ) # Parse command line myargs = parser.parse_args() # Model Names src_Name = myargs.atm ; dst_Name = myargs.oce # Default vars if myargs.repartition_var == None : # Runoff from Antarctica iceshelves (Depoorter, 2013) if myargs.type == 'iceshelf' : repartitionVar = 'sornfisf' # Runoff from Antarctica iceberg (Depoorter, 2013) if myargs.type == 'iceberg' : repartitionVar = 'Icb_Flux' else : repartitionVar = myargs.repartition_var # Latitude limits of each calving zone limit_lat = ( (40.0, 90.0), (-50.0, 40.0), ( -90.0, -50.0) ) nb_zone = len(limit_lat) if myargs.fmt == 'classic' : FmtNetcdf = 'CLASSIC' if myargs.fmt == 'netcdf3' : FmtNetcdf = 'CLASSIC' if myargs.fmt == '64bit' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' if myargs.fmt == '64bit_data' : FmtNetcdf = 'NETCDF3_64BIT_DATA' if myargs.fmt == '64bit_offset' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' if myargs.fmt == 'netcdf4' : FmtNetcdf = 'NETCDF4' if myargs.fmt == 'netcdf4_classic' : FmtNetcdf = 'NETCDF4_CLASSIC' ### ========================================================================================== # Model short names src_name = None ; dst_name = None if src_Name.count('ICO') != 0 : src_name = 'ico' if src_Name.count('LMD') != 0 : src_name = 'lmd' if dst_Name.count('ORCA') != 0 : dst_name = 'orc' CplModel = dst_Name + 'x' + src_Name print ('Atm names : ' + src_name + ' ' + src_Name ) print ('Oce names : ' + dst_name + ' ' + dst_Name ) print (' ') # Reading domains characteristics areas = myargs.dir + '/areas_' + dst_Name + 'x' + src_Name + '.nc' masks = myargs.dir + '/masks_' + dst_Name + 'x' + src_Name + '.nc' grids = myargs.dir + '/grids_' + dst_Name + 'x' + src_Name + '.nc' print ('Opening ' + areas) f_areas = netCDF4.Dataset ( areas, "r" ) print ('Opening ' + masks) f_masks = netCDF4.Dataset ( masks, "r" ) print ('Opening ' + grids) f_grids = netCDF4.Dataset ( grids, "r" ) src_lon = f_grids.variables ['t' + src_name + '.lon'][:] src_lat = f_grids.variables ['t' + src_name + '.lat'][:] dst_lon = f_grids.variables ['t' + dst_name + '.lon'][:] dst_lat = f_grids.variables ['t' + dst_name + '.lat'][:] src_msk = f_masks.variables ['t' + src_name + '.msk'][:] dst_msk = f_masks.variables ['t' + dst_name + '.msk'][:] dst_mskutil = 1-dst_msk # Reversing the convention : 0 on continent, 1 on ocean print ('dst_msk : ' + str(np.sum(dst_msk))) print ('dst_mskutil : ' + str(np.sum(dst_mskutil))) # Ocean grid periodicity nperio_dst=myargs.ocePerio # Periodicity masking for NEMO if nperio_dst == 0 : if dst_Name == 'ORCA2.3' : nperio_dst = 4 if dst_Name == 'eORCA1.2' : nperio_dst = 6 if dst_Name == 'ORCA025' : nperio_dst = 6 if dst_Name == 'eORCA025' : nperio_dst = 6 if dst_Name == 'eORCA025.1' : nperio_dst = 6 print ('nperio_dst: ' + str(nperio_dst) ) dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T', sval=0 ) print ('dst_mskutil : ' + str(np.sum(dst_mskutil))) ## ## Fill Closed seas and other ## ================================================== # Preparation by closing some straits # ----------------------------------- if dst_Name in [ 'eORCA025', 'eORCA025.1', 'eORCA025.4' ] : print ('Filling some seas for ' + dst_Name ) # Set Gibraltar strait to 0 to fill Mediterranean sea dst_mskutil[838:841,1125] = 0 # Set Bal-El-Mandeb to 0 to fill Red Sea dst_mskutil[736,1321:1324] = 0 # Set Stagerak to 0 to fill Baltic Sea dst_mskutil[961,1179:1182] = 0 # Set Ormuz Strait to 0 to fill Arabo-Persian Gulf dst_mskutil[794:798,1374] = 0 # Set Hudson Strait to 0 to fill Hudson Bay dst_mskutil[997:1012,907] = 0 if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4' ] : print ('Filling some seas for ' + dst_Name) # Set Gibraltar strait to 0 to fill Mediterrean sea dst_mskutil[240, 283] = 0 # Set Bal-El-Mandeb to 0 to fill Red Sea dst_mskutil[211:214, 332] = 0 # Set Stagerak to 0 to fill Baltic Sea dst_mskutil[272:276, 293] = 0 # Set Ormuz Strait to 0 to fill Arabo-Persian Gulf dst_mskutil[227:230, 345] = 0 # Set Hudson Strait to 0 to fill Hudson Bay dst_mskutil[284,222:225] = 0 if dst_Name in [ 'ORCA2.3', 'ORCA2.4' ] : print ('Filling some seas for ' + dst_Name) # Set Gibraltar strait to 0 to fill Mediterrean sea dst_mskutil[101,139] = 0 # Set Black Sea to zero. At the edge of the domain : binary_fill_holes fails dst_mskutil[ 99:111, 0: 5] = 0 dst_mskutil[106:111, 173:182] = 0 # Set Stagerak to 0 to fill Baltic Sea dst_mskutil[115,144] = 0 # Set Hudson Strait to 0 to fill Hudson Bay dst_mskutil[120:123,110] = 0 # Set Bal-El-Mandeb to 0 to fill Red Sea dst_mskutil[87:89,166] = 0 dst_closed = dst_mskutil # Fill closed seas with image processing library # ---------------------------------------------- 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) # Surfaces src_srf = f_areas.variables ['t' + src_name + '.srf'] dst_srf = f_areas.variables ['t' + dst_name + '.srf'] dst_srfutil = dst_srf * np.float64 (dst_mskutil) dst_srfutil_sum = np.sum( dst_srfutil) src_clo = f_grids.variables ['t' + src_name + '.clo'][:] src_cla = f_grids.variables ['t' + src_name + '.cla'][:] dst_clo = f_grids.variables ['t' + dst_name + '.clo'][:] dst_cla = f_grids.variables ['t' + dst_name + '.cla'][:] # Indices ( src_jpj, src_jpi) = src_lat.shape ; src_grid_size = src_jpj*src_jpi ( dst_jpj, dst_jpi) = dst_lat.shape ; dst_grid_size = dst_jpj*dst_jpi orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 # Fortran indices (starting at one) ### ===== Reading needed data ================================================== if myargs.type in ['iceberg', 'iceshelf' ]: # Reading data file for calving or iceberg geometry around Antarctica print ( 'Opening ' + myargs.repartition_file) f_repartition = netCDF4.Dataset ( myargs.repartition_file, "r" ) repartition = np.sum ( f_repartition.variables [repartitionVar][:], axis=0 ) ## Before loop on basins remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) src_address = np.empty ( shape=(0), dtype=np.int32 ) dst_address = np.empty ( shape=(0), dtype=np.int32 ) print (' ') ### ===== Starting loop on basins============================================== # Initialise some fields remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) src_address = np.empty ( shape=(0), dtype=np.int32 ) dst_address = np.empty ( shape=(0), dtype=np.int32 ) basin_msk = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float32) key_repartition = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float64) ## Loop for n_bas in np.arange ( nb_zone ) : south = False ; ok_run = False lat_south = np.min(limit_lat[n_bas]) ; lat_north = np.max(limit_lat[n_bas]) if lat_south <= -60.0 : south = True print ( 'basin: {:2d} -- Latitudes: {:+.1f} {:+.1f} --'.format(n_bas, lat_south, lat_north) ) ## if myargs.type == 'iceberg' and south : ok_run = True ; print ('Applying iceberg to south' ) if myargs.type == 'iceshelf' and south : ok_run = True ; print ('Applying iceshelf to south' ) if myargs.type == 'iceberg' and not south : ok_run = False ; print ('Skipping iceberg: not south ') if myargs.type == 'iceshelf' and not south : ok_run = False ; print ('Skipping iceshelf: not south ') if myargs.type == 'nosouth' and south : ok_run = False ; print ('Skipping south: nosouth case' ) if myargs.type == 'nosouth' and not south : ok_run = True ; print ('Running: not in south, uniform repartition') if myargs.type == 'full' : ok_run = True ; print ('Running general trivial case, uniform repartition' ) if ok_run : # Calving integral send to one point per latitude bands index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1 # Basin mask basin_msk[n_bas] = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 ) # Repartition pattern if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition[n_bas] = repartition * basin_msk[n_bas] else : key_repartition[n_bas] = basin_msk[n_bas] # Integral and normalisation sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil ) key_repartition = key_repartition / sum_repartition print ( 'Sum of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] )) ) print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] * dst_srfutil )) ) # Basin surface (modulated by repartition key) basin_srfutil = np.sum ( key_repartition[n_bas] * dst_srfutil ) # Weights and links poids = 1.0 / basin_srfutil matrix_local = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids, 0. ) matrix_local = matrix_local[matrix_local > 0.0] # Keep only non zero values num_links = np.count_nonzero (matrix_local) # Address on source grid : all links points to the same LMDZ point. src_address_local = np.ones(num_links, dtype=np.int32 )*index_src # Address on destination grid : each NEMO point with non zero link dst_address_local = np.where ( key_repartition[n_bas].ravel() > 0.0, orc_index, 0) dst_address_local = dst_address_local[dst_address_local > 0] # Append to global tabs remap_matrix = np.append ( remap_matrix, matrix_local ) src_address = np.append ( src_address , src_address_local ) dst_address = np.append ( dst_address , dst_address_local ) # #print ( 'Sum of remap_matrix : {:12.3e}'.format(np.sum(matrix_local)) ) print ( 'Point in atmospheric grid : {:4d} -- num_links: {:6d}'.format(index_src, num_links) ) print (' ') ## End of loop print (' ') # num_links = np.count_nonzero (remap_matrix) print ( 'Total num_links : {:10d}'.format(num_links) ) # Diag : interpolate uniform field src_field = np.zeros(shape=(src_grid_size)) for n_bas in np.arange ( nb_zone ) : index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1 src_field[index_src-1] = n_bas dst_field = np.zeros(shape=(dst_grid_size,)) for link in np.arange (num_links) : dst_field [dst_address [link]-1] += remap_matrix[link] * src_field[src_address[link]-1] ### ===== Writing the weights file, for OASIS MCT ========================================== # Output file calving = myargs.output f_calving = netCDF4.Dataset ( calving, 'w', format=FmtNetcdf ) print ('Output file: ' + calving ) f_calving.Conventions = "CF-1.6" f_calving.source = "IPSL Earth system model" f_calving.group = "ICMC IPSL Climate Modelling Center" f_calving.Institution = "IPSL https.//www.ipsl.fr" f_calving.Ocean = dst_Name + " https://www.nemo-ocean.eu" f_calving.Atmosphere = src_Name + " http://lmdz.lmd.jussieu.fr" if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.repartition_file f_calving.associatedFiles = grids + " " + areas + " " + masks f_calving.directory = os.getcwd () f_calving.description = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX" f_calving.title = calving f_calving.Program = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) f_calving.repartitionType = myargs.type if myargs.type in [ 'iceberg', 'iceshelf' ] : f_calving.repartitionFile = myargs.repartition_file f_calving.repartitionVar = repartitionVar f_calving.gridsFile = grids f_calving.areasFile = areas f_calving.masksFile = masks f_calving.timeStamp = time.asctime() f_calving.HOSTNAME = platform.node() #f_calving.LOGNAME = os.getlogin() f_calving.Python = "Python version " + platform.python_version() f_calving.OS = platform.system() f_calving.release = platform.release() f_calving.hardware = platform.machine() f_calving.conventions = "SCRIP" if src_name == 'lmd' : f_calving.source_grid = "curvilinear" if src_name == 'ico' : f_calving.source_grid = "unstructured" f_calving.dest_grid = "curvilinear" f_calving.Model = "IPSL CM6" f_calving.SVN_Author = "$Author$" f_calving.SVN_Date = "$Date$" f_calving.SVN_Revision = "$Revision$" f_calving.SVN_Id = "$Id$" f_calving.SVN_HeadURL = "$HeadURL$" d_nb_zone = f_calving.createDimension ('nb_zone' , nb_zone ) d_num_links = f_calving.createDimension ('num_links' , num_links ) d_num_wgts = f_calving.createDimension ('num_wgts' , 1 ) d_src_grid_size = f_calving.createDimension ('src_grid_size' , src_grid_size ) d_src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0] ) d_src_grid_rank = f_calving.createDimension ('src_grid_rank' , 2 ) d_dst_grid_size = f_calving.createDimension ('dst_grid_size' , dst_grid_size ) d_dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] ) d_dst_grid_rank = f_calving.createDimension ('dst_grid_rank' , 2 ) v_remap_matrix = f_calving.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') ) v_src_address = f_calving.createVariable ( 'src_address' , 'i4', ('num_links',) ) v_src_address.convention = "Fortran style addressing, starting at 1" v_dst_address = f_calving.createVariable ( 'dst_address' , 'i4', ('num_links',) ) v_dst_address.convention = "Fortran style addressing, starting at 1" v_remap_matrix[:] = remap_matrix v_src_address [:] = src_address v_dst_address [:] = dst_address v_src_grid_dims = f_calving.createVariable ( 'src_grid_dims' , 'i4', ('src_grid_rank',) ) v_src_grid_center_lon = f_calving.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) ) v_src_grid_center_lat = f_calving.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) ) v_src_grid_center_lon.units='degrees_east' ; v_src_grid_center_lon.long_name='Longitude' v_src_grid_center_lon.long_name='longitude' ; v_src_grid_center_lon.bounds="src_grid_corner_lon" v_src_grid_center_lat.units='degrees_north' ; v_src_grid_center_lat.long_name='Latitude' v_src_grid_center_lat.long_name='latitude ' ; v_src_grid_center_lat.bounds="src_grid_corner_lat" v_src_grid_corner_lon = f_calving.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners') ) v_src_grid_corner_lat = f_calving.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners') ) v_src_grid_corner_lon.units="degrees_east" v_src_grid_corner_lat.units="degrees_north" v_src_grid_area = f_calving.createVariable ( 'src_grid_area' , 'f8', ('src_grid_size',) ) v_src_grid_area.long_name="Grid area" ; v_src_grid_area.standard_name="cell_area" ; v_src_grid_area.units="m2" v_src_grid_imask = f_calving.createVariable ( 'src_grid_imask' , 'f8', ('src_grid_size',) ) v_src_grid_imask.long_name="Land-sea mask" ; v_src_grid_imask.units="Land:1, Ocean:0" v_src_grid_dims [:] = ( src_jpi, src_jpi ) v_src_grid_center_lon[:] = src_lon[:].ravel() v_src_grid_center_lat[:] = src_lat[:].ravel() v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) v_src_grid_area [:] = src_srf[:].ravel() v_src_grid_imask [:] = src_msk[:].ravel() # -- v_dst_grid_dims = f_calving.createVariable ( 'dst_grid_dims' , 'i4', ('dst_grid_rank',) ) v_dst_grid_center_lon = f_calving.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) ) v_dst_grid_center_lat = f_calving.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) ) v_dst_grid_center_lon.units='degrees_east' ; v_dst_grid_center_lon.long_name='Longitude' v_dst_grid_center_lon.long_name='longitude' ; v_dst_grid_center_lon.bounds="dst_grid_corner_lon" v_dst_grid_center_lat.units='degrees_north' ; v_dst_grid_center_lat.long_name='Latitude' v_dst_grid_center_lat.long_name='latitude' ; v_dst_grid_center_lat.bounds="dst_grid_corner_lat" v_dst_grid_corner_lon = f_calving.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners') ) v_dst_grid_corner_lat = f_calving.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners') ) v_dst_grid_corner_lon.units="degrees_east" v_dst_grid_corner_lat.units="degrees_north" v_dst_grid_area = f_calving.createVariable ( 'dst_grid_area' , 'f8', ('dst_grid_size',) ) v_dst_grid_area.long_name="Grid area" ; v_dst_grid_area.standard_name="cell_area" ; v_dst_grid_area.units="m2" v_dst_grid_imask = f_calving.createVariable ( 'dst_grid_imask' , 'f8', ('dst_grid_size',) ) v_dst_grid_imask.long_name="Land-sea mask" ; v_dst_grid_imask.units="Land:1, Ocean:0" v_dst_bassin = f_calving.createVariable ( 'dst_bassin' , 'f8', ('nb_zone', 'dst_grid_size',) ) v_dst_repartition = f_calving.createVariable ( 'dst_repartition' , 'f8', ('nb_zone', 'dst_grid_size',) ) v_dst_southLimit = f_calving.createVariable ('dst_southLimit', 'f4', ('nb_zone',) ) v_dst_northLimit = f_calving.createVariable ('dst_northLimit', 'f4', ('nb_zone',) ) v_dst_southLimit[:] = np.min(limit_lat, axis=(1,) ) v_dst_northLimit[:] = np.max(limit_lat, axis=(1,) ) v_dst_grid_dims [:] = ( dst_jpi, dst_jpi ) v_dst_grid_center_lon[:] = dst_lon[:].ravel() v_dst_grid_center_lat[:] = dst_lat[:].ravel() v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) v_dst_grid_area [:] = dst_srf[:].ravel() v_dst_grid_imask [:] = dst_msk[:].ravel() v_dst_bassin[:] = basin_msk.reshape ( (nb_zone,dst_grid_size) ) v_dst_repartition[:] = key_repartition.reshape( (nb_zone,dst_grid_size) ) # Additionnal fields for debugging # ================================ v_dst_grid_imaskutil = f_calving.createVariable ( 'dst_grid_imaskutil' , 'f8', ('dst_grid_size',) ) v_dst_grid_imaskutil.long_name="Land-sea mask" ; v_dst_grid_imaskutil.units="Land:1, Ocean:0" v_dst_grid_imaskclose = f_calving.createVariable ( 'dst_grid_imaskclose' , 'f8', ('dst_grid_size',) ) v_dst_grid_imaskclose.long_name="Land-sea mask" ; v_dst_grid_imaskclose.units="Land:1, Ocean:0" v_dst_grid_imaskutil [:] = dst_mskutil[:].ravel() v_dst_grid_imaskclose[:] = dst_closed[:].ravel() v_dst_lon_addressed = f_calving.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) ) v_dst_lat_addressed = f_calving.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) ) v_dst_lon_addressed.long_name="Longitude" ; v_dst_lon_addressed.standard_name="longitude" ; v_dst_lon_addressed.units="degrees_east" v_dst_lat_addressed.long_name="Latitude" ; v_dst_lat_addressed.standard_name="latitude" ; v_dst_lat_addressed.units="degrees_north" v_dst_area_addressed = f_calving.createVariable ( 'dst_area_addressed', 'f8', ('num_links',) ) v_dst_area_addressed.long_name="Grid area" ; v_dst_area_addressed.standard_name="cell_area" ; v_dst_grid_area.units="m2" v_dst_imask_addressed = f_calving.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) ) v_dst_imask_addressed.long_name="Land-sea mask" ; v_dst_imask_addressed.units="Land:1, Ocean:0" v_dst_imaskutil_addressed = f_calving.createVariable ( 'dst_imaskutil_addressed', 'i4', ('num_links',) ) v_dst_imaskutil_addressed.long_name="Land-sea mask" ; v_dst_imaskutil_addressed.units="Land:1, Ocean:0" v_dst_lon_addressed [:] = dst_lon.ravel()[dst_address-1].ravel() v_dst_lat_addressed [:] = dst_lat.ravel()[dst_address-1].ravel() v_dst_area_addressed [:] = dst_srf[:].ravel()[dst_address-1].ravel() v_dst_imask_addressed[:] = dst_msk[:].ravel()[dst_address-1].ravel() v_dst_imaskutil_addressed[:] = dst_mskutil.ravel()[dst_address-1].ravel() v_src_field = f_calving.createVariable ( 'src_field', 'f8', ('src_grid_size',) ) v_dst_field = f_calving.createVariable ( 'dst_field', 'f8', ('dst_grid_size',) ) v_src_field[:] = src_field v_dst_field[:] = dst_field f_calving.close () ### ===== That's all Folk's !! ==============================================