Ignore:
Timestamp:
02/24/22 15:17:45 (2 years ago)
Author:
omamce
Message:

O.M. : changes in MOSAIX

File:
1 edited

Legend:

Unmodified
Added
Removed

  • TOOLS/MOSAIX/CalvingWeights.py

    r6064 r6066  
    1414##  personal. 
    1515## 
    16 import netCDF4, numpy as np 
     16import numpy as np, xarray as xr 
    1717import sys, os, platform, argparse, textwrap, time 
    1818from scipy import ndimage 
    … …  
    3636 
    3737# Adding arguments 
    38 parser.add_argument ('--oce'     , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025', 'eORCA025.1'] ) 
     38parser.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'] ) 
    3940parser.add_argument ('--atm'     , help='atm model name (ICO* or LMD*)', type=str, default='ICO40'    ) 
    4041parser.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='./' ) 
    4243parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' ) 
    4344parser.add_argument ('--repartition_var' , help='Variable name for iceshelf'      , type=str, default=None) 
     45parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) 
     46parser.add_argument ('--areas' , help='masks file name', default='areas.nc' ) 
     47parser.add_argument ('--masks' , help='areas file name', default='masks.nc' ) 
     48parser.add_argument ('--o2a'   , help='o2a file name'  , default='o2a.nc'   ) 
    4449parser.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 ) 
     50parser.add_argument ('--fmt'     , help='NetCDF file format, using nco syntax', default='netcdf4', 
     51                         choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 
     52parser.add_argument ('--ocePerio', help='periodicity of ocean grid', type=int, default=0 ) 
    4753 
    4854# Parse command line 
    … …  
    8995 
    9096# 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' 
     97grids = myargs.grids 
     98areas = myargs.areas 
     99masks = myargs.masks 
     100o2a   = myargs.o2a 
    94101 
    95102print ('Opening ' + areas) 
    96 f_areas = netCDF4.Dataset ( areas, "r" ) 
     103f_areas = xr.open_dataset ( areas ) 
    97104print ('Opening ' + masks) 
    98 f_masks = netCDF4.Dataset ( masks, "r" ) 
     105f_masks = xr.open_dataset ( masks ) 
    99106print ('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'][:] 
     107f_grids = xr.open_dataset ( grids ) 
     108 
     109src_lon = f_grids ['t' + src_name + '.lon'] 
     110src_lat = f_grids ['t' + src_name + '.lat'] 
     111dst_lon = f_grids ['t' + dst_name + '.lon'] 
     112dst_lat = f_grids ['t' + dst_name + '.lat'] 
     113 
     114src_msk = f_masks ['t' + src_name + '.msk'] 
     115dst_msk = f_masks ['t' + dst_name + '.msk'] 
    109116dst_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))) 
     117print ('dst_msk     : ', dst_msk.sum().values ) 
     118print ('dst_mskutil : ', dst_mskutil.sum().values ) 
    112119 
    113120# Ocean grid periodicity 
    … …  
    116123# Periodicity masking for NEMO 
    117124if 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         
    123129print ('nperio_dst: ' + str(nperio_dst) ) 
    124130dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T', sval=0 ) 
    … …  
    144150    dst_mskutil[997:1012,907] = 0 
    145151     
    146 if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4' ] : 
     152if dst_Name in [ 'eORCA1.2', 'eORCA1.3', 'eORCA1.4', 'eORCA1.4.2' ] : 
    147153    print ('Filling some seas for ' + dst_Name) 
    148154    # Set Gibraltar strait to 0 to fill Mediterrean sea 
    … …  
    184190dst_srfutil_sum = np.sum( dst_srfutil) 
    185191 
    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'][:] 
     192src_clo = f_grids ['t' + src_name + '.clo'] 
     193src_cla = f_grids ['t' + src_name + '.cla'] 
     194dst_clo = f_grids ['t' + dst_name + '.clo'] 
     195dst_cla = f_grids ['t' + dst_name + '.cla'] 
    190196 
    191197# Indices 
    … …  
    198204    # Reading data file for calving or iceberg geometry around Antarctica 
    199205    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 ) 
    202208 
    203209## Before loop on basins 
    … …  
    246252 
    247253        # 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 
    249255        key_repartition = key_repartition / sum_repartition 
    250256     
    251257        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) ) 
    253259     
    254260        # Basin surface (modulated by repartition key) 
    … …  
    257263        # Weights and links 
    258264        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. ) 
    260266        matrix_local   = matrix_local[matrix_local > 0.0] # Keep only non zero values 
    261267        num_links      = np.count_nonzero (matrix_local) 
    … …  
    297303# Output file 
    298304calving = myargs.output 
    299 f_calving = netCDF4.Dataset ( calving, 'w', format=FmtNetcdf ) 
    300305     
    301306print ('Output file: ' + calving ) 
    302307 
    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() 
     308remap_matrix = xr.DataArray (np.reshape(remap_matrix, (num_links, 1)), dims = ['num_links', 'num_wgts'] ) 
     309src_address  = xr.DataArray (src_address.astype(np.int32) , dims = ['num_links',] ) 
     310src_address.attrs['convention'] = "Fortran style addressing, starting at 1" 
     311dst_address  = xr.DataArray (dst_address.astype(np.int32) , dims =  ['num_links',] ) 
     312dst_address.attrs['convention'] = "Fortran style addressing, starting at 1" 
     313 
     314src_grid_dims       = xr.DataArray ( np.array (( src_jpi, src_jpi ), dtype=np.int32), dims =  ['src_grid_rank',] ) 
     315src_grid_center_lon = xr.DataArray ( src_lon.values.ravel(), dims = ['src_grid_size',] ) 
     316src_grid_center_lat = xr.DataArray ( src_lat.values.ravel(), dims = ['src_grid_size',] ) 
     317src_grid_center_lon.attrs['units']='degrees_east'  ; src_grid_center_lon.attrs['long_name']='Longitude' 
     318src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon"  
     319src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' 
     320src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat"  
     321src_grid_corner_lon = xr.DataArray ( src_clo.values.reshape( (src_jpi*src_jpj, src_clo.shape[0]) ), dims = ['src_grid_size', 'src_grid_corners'] ) 
     322src_grid_corner_lat = xr.DataArray ( src_cla.values.reshape( (src_jpi*src_jpj, src_clo.shape[0]) ), dims = ['src_grid_size', 'src_grid_corners' ] ) 
     323src_grid_corner_lon.attrs['units']="degrees_east" 
     324src_grid_corner_lat.attrs['units']="degrees_north" 
     325src_grid_area       = xr.DataArray ( src_srf.values.ravel(), dims = ['src_grid_size',] ) 
     326src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" 
     327src_grid_imask      = xr.DataArray ( src_msk.values.ravel().astype(np.int32) , dims = ['src_grid_size',] ) 
     328src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" 
    385329 
    386330# -- 
    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) ) 
     331dst_grid_dims       = xr.DataArray ( np.array(( dst_jpi, dst_jpi), dtype=np.int32) , dims = ['dst_grid_rank', ]) 
     332dst_grid_center_lon = xr.DataArray ( dst_lon.values.ravel(), dims = ['dst_grid_size', ]) 
     333dst_grid_center_lat = xr.DataArray ( dst_lat.values.ravel(), dims = ['dst_grid_size', ]) 
     334dst_grid_center_lon.attrs['units']='degrees_east'  ; dst_grid_center_lon.attrs['long_name']='Longitude' 
     335dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon"  
     336dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' 
     337dst_grid_center_lat.attrs['long_name']='latitude'  ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat"  
     338dst_grid_corner_lon = xr.DataArray ( dst_clo.values.reshape(dst_jpi*dst_jpj, dst_clo.shape[0] ), dims = ['dst_grid_size', 'dst_grid_corners'] ) 
     339dst_grid_corner_lat = xr.DataArray ( dst_cla.values.reshape(dst_jpi*dst_jpj, dst_clo.shape[0] ), dims = ['dst_grid_size', 'dst_grid_corners'] ) 
     340dst_grid_corner_lon.attrs['units']="degrees_east" 
     341dst_grid_corner_lat.attrs['units']="degrees_north" 
     342dst_grid_area       = xr.DataArray ( dst_srf.values.ravel(), dims = ['dst_grid_size', ] ) 
     343dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 
     344dst_grid_imask      = xr.DataArray ( dst_msk.values.ravel(), dims =  ['dst_grid_size',]  ) 
     345dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" 
     346 
     347dst_bassin      = xr.DataArray ( basin_msk.reshape ( (nb_zone,dst_grid_size) ) , dims = ['nb_zone', 'dst_grid_size',  ] ) 
     348dst_repartition = xr.DataArray ( key_repartition.reshape( (nb_zone,dst_grid_size) ) , dims = ['nb_zone', 'dst_grid_size',  ] ) 
     349 
     350dst_southLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , dims = ['nb_zone', ] ) 
     351dst_northLimit = xr.DataArray ( np.min(limit_lat, axis=(1,)).astype(np.float32) , dims = ['nb_zone', ] ) 
    422352 
    423353# Additionnal fields for debugging 
    424354# ================================ 
    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 !! ============================================== 
     355dst_grid_imaskutil  = xr.DataArray ( dst_mskutil.ravel().astype(np.float32), dims = ['dst_grid_size',] ) 
     356dst_grid_imaskutil.attrs['long_name']="Land-sea mask" ; dst_grid_imaskutil.attrs['units']="Land:1, Ocean:0" 
     357dst_grid_imaskclose = xr.DataArray ( dst_closed.values.ravel(), dims = ['dst_grid_size',] ) 
     358dst_grid_imaskclose.attrs['long_name']="Land-sea mask" ; dst_grid_imaskclose.attrs['units']="Land:1, Ocean:0" 
     359 
     360dst_lon_addressed = xr.DataArray ( dst_lon.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 
     361dst_lat_addressed = xr.DataArray ( dst_lat.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 
     362dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" 
     363dst_lat_addressed.attrs['long_name']="Latitude"  ; dst_lat_addressed.attrs['standard_name']="latitude"  ; dst_lat_addressed.attrs['units']="degrees_north" 
     364dst_area_addressed  = xr.DataArray ( dst_srf.values.ravel()[dst_address-1].ravel(), dims = ['num_links',] ) 
     365dst_area_addressed.attrs['long_name']="Grid area" ; dst_area_addressed.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 
     366dst_imask_addressed = xr.DataArray ( dst_msk.values.ravel()[dst_address-1].ravel(), dims = ['num_links', ] ) 
     367dst_imask_addressed.attrs['long_name']="Land-sea mask" ; dst_imask_addressed.attrs['units']="Land:1, Ocean:0" 
     368dst_imaskutil_addressed = xr.DataArray ( dst_mskutil.ravel()[dst_address-1].astype(np.float32), dims = ['num_links'] ) 
     369dst_imaskutil_addressed.attrs['long_name']="Land-sea mask" ; dst_imaskutil_addressed.attrs['units']="Land:1, Ocean:0" 
     370 
     371src_field = xr.DataArray ( src_field, dims = ['src_grid_size',] ) 
     372dst_field = xr.DataArray ( dst_field, dims = ['dst_grid_size',] ) 
     373 
     374f_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 
     407f_calving.attrs['Conventions']     = "CF-1.6" 
     408f_calving.attrs['source']          = "IPSL Earth system model" 
     409f_calving.attrs['group']           = "ICMC IPSL Climate Modelling Center" 
     410f_calving.attrs['Institution']     = "IPSL https.//www.ipsl.fr" 
     411f_calving.attrs['Ocean']           = dst_Name + " https://www.nemo-ocean.eu" 
     412f_calving.attrs['Atmosphere']      = src_Name + " http://lmdz.lmd.jussieu.fr" 
     413if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.attrs['originalFiles'] = myargs.repartition_file 
     414f_calving.attrs['associatedFiles'] = grids + " " + areas + " " + masks 
     415f_calving.attrs['directory']       = os.getcwd () 
     416f_calving.attrs['description']     = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX" 
     417f_calving.attrs['title']           = calving 
     418f_calving.attrs['Program']         = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) 
     419f_calving.attrs['repartitionType'] = myargs.type 
     420if myargs.type in [ 'iceberg', 'iceshelf' ] : 
     421    f_calving.attrs['repartitionFile'] = myargs.repartition_file 
     422    f_calving.attrs['repartitionVar']  = repartitionVar 
     423f_calving.attrs['gridsFile']       = grids 
     424f_calving.attrs['areasFile']       = areas 
     425f_calving.attrs['masksFile']       = masks 
     426f_calving.attrs['timeStamp']       = time.asctime() 
     427f_calving.attrs['HOSTNAME']        = platform.node() 
     428f_calving.attrs['LOGNAME']         = os.getlogin() 
     429f_calving.attrs['Python']          = "Python version " +  platform.python_version() 
     430f_calving.attrs['OS']              = platform.system() 
     431f_calving.attrs['release']         = platform.release() 
     432f_calving.attrs['hardware']        = platform.machine() 
     433f_calving.attrs['conventions']     = "SCRIP" 
     434if src_name == 'lmd' : f_calving.attrs['source_grid'] = "curvilinear" 
     435if src_name == 'ico' : f_calving.attrs['source_grid'] = "unstructured" 
     436f_calving.attrs['dest_grid']       = "curvilinear" 
     437f_calving.attrs['Model']           = "IPSL CM6" 
     438f_calving.attrs['SVN_Author']      = "$Author$" 
     439f_calving.attrs['SVN_Date']        = "$Date$" 
     440f_calving.attrs['SVN_Revision']    = "$Revision$" 
     441f_calving.attrs['SVN_Id']          = "$Id$" 
     442f_calving.attrs['SVN_HeadURL']     = "$HeadURL$" 
     443 
     444## 
     445f_calving.to_netcdf ( calving, mode='w', format=FmtNetcdf ) 
     446f_calving.close() 
     447 
     448print ('That''s all folks !') 
     449## ====================================================================================== 
Note: See TracChangeset for help on using the changeset viewer.