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/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.