Changeset 6066


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

O.M. : changes in MOSAIX

Location:
TOOLS/MOSAIX
Files:
3 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## ====================================================================================== 

  • TOOLS/MOSAIX/CreateWeightsMask.bash

    r6064 r6066  
    6969# ================================= 
    7070#CplModel=ORCA2.3xLMD9695 
    71 CplModel=ORCA2.3xICO30 
    72 #CplModel=ORCA2.3xICO40 
     71#CplModel=ORCA2.3xICO30 
     72CplModel=ORCA2.3xICO40 
    7373#CplModel=eORCA1.2xLMD144142 
    7474#CplModel=eORCA1.2xLMD256256 
    … …  
    800800              --atmQuantity=${runOff_atmQuantity} --oceQuantity=${runOff_oceQuantity} --ocePerio=${OcePerio} 
    801801fi 
    802 exit 
     802 
    803803## 
    804804echo ${Titre}"Calving weights"${Norm} 
    … …  
    809809        cp /ccc/work/cont003/gencmip6/deshayej/eORCA_R025_runoff_v1.2.nc . 
    810810        ${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 
    812814        ${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 
    814818        ${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 
    816822        ;; 
    817823         
    818         ( eORCA1.2 ) 
     824        ( eORCA1.2 | eORCA1.4 ) 
    819825        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 
    826838        ;; 
    827839         
    828840        ( * ) 
    829841        ${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 
    831845        ;; 
    832846 

  • TOOLS/MOSAIX/RunoffWeights.py

    r6064 r6066  
    2626 
    2727## Modules 
    28 import netCDF4 
    29 import numpy as np 
     28import numpy as np, xarray as xr 
    3029import nemo 
    3130from scipy import ndimage 
    … …  
    5554 
    5655# Adding arguments 
    57 parser.add_argument ('--oce'          , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025', 'eORCA025.1'] ) 
     56parser.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'] ) 
    5858parser.add_argument ('--atm'          , help='atm model name', type=str, default='LMD9695'    ) 
    5959parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', type=int, default=1 ) 
    6060parser.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'] ) 
     61parser.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'] ) 
     63parser.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'] ) 
    6365parser.add_argument ('--searchRadius' , help='max distance to connect a land point to an ocean point (in km)', type=float, default=600.0 ) 
    6466parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) 
    … …  
    6769parser.add_argument ('--o2a'   , help='o2a file name'  , default='o2a.nc'   ) 
    6870parser.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'] ) 
     71parser.add_argument ('--fmt'   , help='NetCDF file format, using nco syntax', default='netcdf4', 
     72                         choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) 
    7073parser.add_argument ('--ocePerio'   , help='periodicity of ocean grid', type=int, default=0 ) 
    7174 
    … …  
    104107 
    105108# Ocean grid periodicity 
    106 oce_perio=myargs.ocePerio 
     109oce_perio = myargs.ocePerio 
    107110 
    108111### Read coordinates of all models 
    109112### 
    110113 
    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() 
     114diaFile    = xr.open_dataset ( o2a   ) 
     115gridFile   = xr.open_dataset ( grids ) 
     116areaFile   = xr.open_dataset ( areas ) 
     117maskFile   = xr.open_dataset ( masks ) 
     118 
     119o2aFrac             = diaFile ['OceFrac'].squeeze() 
    117120o2aFrac = np.where ( np.abs(o2aFrac) < 1E10, o2aFrac, 0.0) 
    118121 
    … …  
    121124atm_grid_rank = len(gridFile['t'+atm_n+'.lat'][:].shape) 
    122125 
    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() 
     126atm_grid_center_lat = gridFile['t'+atm_n+'.lat'].squeeze() 
     127atm_grid_center_lon = gridFile['t'+atm_n+'.lon'].squeeze() 
     128atm_grid_corner_lat = gridFile['t'+atm_n+'.cla'].squeeze() 
     129atm_grid_corner_lon = gridFile['t'+atm_n+'.clo'].squeeze() 
     130atm_grid_area       = areaFile['t'+atm_n+'.srf'].squeeze() 
     131atm_grid_imask      = 1-maskFile['t'+atm_n+'.msk'][:].squeeze() 
    129132atm_grid_dims       = gridFile['t'+atm_n+'.lat'][:].shape 
     133 
     134if 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     
     142if 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']) 
    130149 
    131150atm_perio = 0 
    … …  
    137156oce_grid_rank = len(gridFile['torc.lat'][:].shape) 
    138157 
    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() 
     158oce_grid_center_lat = gridFile['torc.lat'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 
     159oce_grid_center_lon = gridFile['torc.lon'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 
     160oce_grid_corner_lat = gridFile['torc.cla'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 
     161oce_grid_corner_lon = gridFile['torc.clo'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 
     162oce_grid_area       = areaFile['torc.srf'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 
     163oce_grid_imask      = 1-maskFile['torc.msk'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) 
    145164oce_grid_dims       = gridFile['torc.lat'][:].shape 
     165 
    146166if oce_perio == 0 : 
    147167    if oce_jpi ==  182 : oce_perio = 4 # ORCA 2 
    148168    if oce_jpi ==  362 : oce_perio = 6 # ORCA 1 
    149169    if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 
     170         
    150171print ("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() 
     172oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel() 
    152173oce_address = np.arange(oce_jpj*oce_jpi) 
    153174 
    154175print ("Fill closed sea with image processing library") 
    155176oce_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 ) 
     177oce_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 ) 
    157179oce_grid_imask = oce_grid_imask2D.ravel() 
    158180##  
    159181print ("Computes an ocean coastal band") 
    160182 
    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) ) 
     183oceLand2D  = np.reshape ( np.where (oce_grid_pmask < 0.5, True, False), (oce_jpj, oce_jpi) ) 
     184oceOcean2D = np.reshape ( np.where (oce_grid_pmask > 0.5, True, False), (oce_jpj, oce_jpi) ) 
    163185 
    164186NNocean = 1+2*oceCoastWidth 
    … …  
    168190 
    169191oceOceanFiltered = oceOceanFiltered2D.ravel() 
    170 oceLand  = oceLand2D.ravel() 
     192oceLand  = oceLand2D.ravel () 
    171193oceOcean = oceOcean2D.ravel() 
    172194oceCoast = oceCoast2D.ravel() 
    … …  
    232254        num_links = int(z_mask.sum()) 
    233255        if num_links == 0 : continue 
    234         z_area = oceCoast_grid_area[z_mask].sum() 
     256        z_area = oceCoast_grid_area[z_mask].sum().values 
    235257        poids = np.ones ((num_links),dtype=np.float64) / z_area 
    236258        if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja] 
    237259        if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] 
    238260        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() ) ) 
    240263        # 
    241264        matrix_local = poids 
    … …  
    254277print ("Write output file") 
    255278runoff = myargs.output 
    256 f_runoff = netCDF4.Dataset ( runoff, 'w', format=FmtNetcdf ) 
    257279print ('Output file: ' + runoff ) 
    258280 
    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 
     282remap_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 
     285src_address  = xr.DataArray ( atm_address.astype(np.int32)+1, dims = ['num_links'], 
     286                                 attrs={"convention": "Fortran style addressing, starting at 1"})  
     287dst_address  = xr.DataArray ( oce_address.astype(np.int32)+1, dims = ['num_links'], 
     288                                 attrs={"convention": "Fortran style addressing, starting at 1"}) 
     289 
     290src_grid_dims       = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32), dims = ['src_grid_rank',] ) 
     291src_grid_center_lon = xr.DataArray (atm_grid_center_lon.values , dims = ['src_grid_size',] ) 
     292src_grid_center_lat = xr.DataArray (atm_grid_center_lat.values , dims = ['src_grid_size',] ) 
     293 
     294src_grid_center_lon.attrs['units']='degrees_east'  ; src_grid_center_lon.attrs['long_name']='Longitude' 
     295src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" 
     296src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' 
     297src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" 
     298  
     299src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) 
     300src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) 
     301src_grid_corner_lon.attrs['units']="degrees_east" 
     302src_grid_corner_lat.attrs['units']="degrees_north" 
     303 
     304src_grid_area       =  xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] )  
     305src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" 
     306 
     307src_grid_imask      =  xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] )  
     308src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" 
     309 
     310src_grid_pmask      =  xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] )  
     311src_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; src_grid_pmask.attrs['units']="Land:1, Ocean:0" 
    342312 
    343313# -- 
    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] 
     314dst_grid_dims       = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32), dims = ['dst_grid_rank',] ) 
     315dst_grid_center_lon = xr.DataArray (oce_grid_center_lon.values, dims = ['dst_grid_size',] ) 
     316dst_grid_center_lat = xr.DataArray (oce_grid_center_lat.values, dims = ['dst_grid_size',] ) 
     317 
     318dst_grid_center_lon.attrs['units']='degrees_east'  ; dst_grid_center_lon.attrs['long_name']='Longitude' 
     319dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" 
     320dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' 
     321dst_grid_center_lat.attrs['long_name']='latitude ' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" 
     322 
     323dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 
     324dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) 
     325dst_grid_corner_lon.attrs['units']="degrees_east" 
     326dst_grid_corner_lat.attrs['units']="degrees_north" 
     327 
     328dst_grid_area       =  xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] )  
     329dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" 
     330 
     331dst_grid_imask      =  xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] )  
     332dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" 
     333 
     334dst_grid_pmask      =  xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] )  
     335dst_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; dst_grid_pmask.attrs['units']="Land:1, Ocean:0" 
     336 
     337src_lon_addressed   =  xr.DataArray (atm_grid_center_lon.values[atm_address]                  , dims = ['num_links'] ) 
     338src_lat_addressed   =  xr.DataArray (atm_grid_center_lat.values[atm_address]                  , dims = ['num_links'] ) 
     339src_area_addressed  =  xr.DataArray (atm_grid_area      .values[atm_address]                  , dims = ['num_links'] ) 
     340src_imask_addressed =  xr.DataArray (1-atm_grid_imask   .values[atm_address].astype(np.int32) , dims = ['num_links'] ) 
     341src_pmask_addressed =  xr.DataArray (1-atm_grid_pmask   .values[atm_address].astype(np.int32) , dims = ['num_links'] ) 
     342 
     343dst_lon_addressed   =  xr.DataArray (oce_grid_center_lon.values[atm_address], dims = ['num_links'] ) 
     344dst_lat_addressed   =  xr.DataArray (oce_grid_center_lat.values[oce_address], dims = ['num_links'] ) 
     345dst_area_addressed  =  xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32)      , dims = ['num_links'] ) 
     346dst_imask_addressed =  xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32)   , dims = ['num_links'] ) 
     347dst_pmask_addressed =  xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32)   , dims = ['num_links'] ) 
     348 
     349src_lon_addressed.attrs['long_name']="Longitude" ; src_lon_addressed.attrs['standard_name']="longitude" ; src_lon_addressed.attrs['units']="degrees_east" 
     350src_lat_addressed.attrs['long_name']="Latitude"  ; src_lat_addressed.attrs['standard_name']="latitude"  ; src_lat_addressed.attrs['units']="degrees_north" 
     351 
     352dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" 
     353dst_lat_addressed.attrs['long_name']="Latitude"  ; dst_lat_addressed.attrs['standard_name']="latitude"  ; dst_lat_addressed.attrs['units']="degrees_north" 
    397354 
    398355if 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 
     363atmCoast         = xr.DataArray (atmCoast.astype(np.int32)          , dims = ['src_grid_size',]) 
     364oceLand          = xr.DataArray (oceLand.astype(np.int32)           , dims = ['dst_grid_size',]) 
     365oceOcean         = xr.DataArray (oceOcean.astype(np.int32)          , dims = ['dst_grid_size',]) 
     366oceOceanFiltered = xr.DataArray (oceOceanFiltered.astype(np.float32), dims = ['dst_grid_size',]) 
     367oceCoast         = xr.DataArray (oceCoast.astype(np.int32)          , dims = ['dst_grid_size',]) 
     368 
     369 
     370f_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 
     404f_runoff.attrs['Conventions']     = "CF-1.6" 
     405f_runoff.attrs['source']          = "IPSL Earth system model" 
     406f_runoff.attrs['group']           = "ICMC IPSL Climate Modelling Center" 
     407f_runoff.attrs['Institution']     = "IPSL https.//www.ipsl.fr" 
     408f_runoff.attrs['Ocean']           = oce_Name + " https://www.nemo-ocean.eu" 
     409f_runoff.attrs['Atmosphere']      = atm_Name + " http://lmdz.lmd.jussieu.fr" 
     410f_runoff.attrs['associatedFiles'] = grids + " " + areas + " " + masks 
     411f_runoff.attrs['directory']       = os.getcwd () 
     412f_runoff.attrs['description']     = "Generated with RunoffWeights.py" 
     413f_runoff.attrs['title']           = runoff 
     414args = '' 
     415for i in np.arange(1,len(sys.argv[1:])) : args += sys.argv[i] + ' ' 
     416f_runoff.attrs['Program']         = "Generated by " + sys.argv[0] + " with flags " + args 
     417f_runoff.attrs['atmCoastWidth']   = "{:d} grid points".format(atmCoastWidth) 
     418f_runoff.attrs['oceCoastWidth']   = "{:d} grid points".format(oceCoastWidth) 
     419f_runoff.attrs['searchRadius']    = "{:.0f} km".format(searchRadius/1000.) 
     420f_runoff.attrs['atmQuantity']     = myargs.atmQuantity 
     421f_runoff.attrs['oceQuantity']     = myargs.oceQuantity 
     422f_runoff.attrs['gridsFile']       = grids 
     423f_runoff.attrs['areasFile']       = areas 
     424f_runoff.attrs['masksFile']       = masks 
     425f_runoff.attrs['o2aFile']         = o2a 
     426f_runoff.attrs['timeStamp']       = time.asctime() 
     427f_runoff.attrs['HOSTNAME']        = platform.node() 
     428f_runoff.attrs['LOGNAME']         = os.getlogin() 
     429f_runoff.attrs['Python']          = "Python version " +  platform.python_version() 
     430f_runoff.attrs['OS']              = platform.system() 
     431f_runoff.attrs['release']         = platform.release() 
     432f_runoff.attrs['hardware']        = platform.machine() 
     433f_runoff.attrs['conventions']     = "SCRIP" 
     434f_runoff.attrs['source_grid']     = "curvilinear" 
     435f_runoff.attrs['dest_grid']       = "curvilinear" 
     436f_runoff.attrs['Model']           = "IPSL CM6" 
     437f_runoff.attrs['SVN_Author']      = "$Author$" 
     438f_runoff.attrs['SVN_Date']        = "$Date$" 
     439f_runoff.attrs['SVN_Revision']    = "$Revision$" 
     440f_runoff.attrs['SVN_Id']          = "$Id$" 
     441f_runoff.attrs['SVN_HeadURL']     = "$HeadURL$" 
     442 
     443f_runoff.to_netcdf ( runoff, mode='w', format=FmtNetcdf ) 
    429444f_runoff.close () 
    430445 
     446## 
     447 
    431448print ('That''s all folks !') 
    432  
    433449## ====================================================================================== 
Note: See TracChangeset for help on using the changeset viewer.