source: TOOLS/MOSAIX/RunoffWeights.py @ 6823

Last change on this file since 6823 was 6666, checked in by omamce, 8 months ago

O.M. : MOSAIX

Improved code with pylint analysis
  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 25.7 KB
Line 
1# -*- Mode: python -*-
2#!/usr/bin/env python3
3### ===========================================================================
4###
5### Compute runoff weights.
6### For LMDZ only. Not suitable for DYNAMICO
7###
8### ===========================================================================
9##
10##  MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt"
11##  file for an english version of the licence and
12##  "Licence_CeCILL_V2-fr.txt" for a french version.
13##
14##  Permission is hereby granted, free of charge, to any person or
15##  organization obtaining a copy of the software and accompanying
16##  documentation covered by this license (the "Software") to use,
17##  reproduce, display, distribute, execute, and transmit the
18##  Software, and to prepare derivative works of the Software, and to
19##  permit third-parties to whom the Software is furnished to do so,
20##  all subject to the following:
21##
22##  Warning, to install, configure, run, use any of MOSAIX software or
23##  to read the associated documentation you'll need at least one (1)
24##  brain in a reasonably working order. Lack of this implement will
25##  void any warranties (either express or implied).  Authors assumes
26##  no responsability for errors, omissions, data loss, or any other
27##  consequences caused directly or indirectly by the usage of his
28##  software by incorrectly or partially configured
29##
30##
31'''Python code to generates weights for run-off. Launched by `CreateWeights.bash`.
32
33More information with python `RunOffWeights.py -h`.
34
35- Defines a coastal band on land in the atmosphere model. For LMDZ
36  rectilinear grids, the width of this band is parametrized, in grid
37  points. For ico grids, the band contains only point with a
38  fraction of ocean in ]0,1[.
39
40- Defines a coastal band on ocean in the ocean model. The width of
41  this band is parametrized.
42
43- An atmosphere point of the coastal band is linked to an ocean
44  point in the coastal band if the distance between the two is less
45  than searchRadius.
46
47## SVN information
48Author   = "$Author$"
49Date     = "$Date$"
50Revision = "$Revision$"
51Id       = "$Id$"
52HeadURL  = "$HeadURL$"
53'''
54
55## Modules
56import sys
57import os
58import platform
59import getpass
60import argparse
61import time
62import numpy as np
63import xarray as xr
64from scipy import ndimage
65import nemo
66
67# SVN information
68__SVN__ = ({
69    'Author'   : '$Author$',
70    'Date'     : '$Date$',
71    'Revision' : '$Revision$',
72    'Id'       : '$Id$',
73    'HeadURL'  : '$HeadURL$',
74    'SVN_Date' : '$SVN_Date: $',
75    })
76##
77
78## Useful constants
79zero    = np.dtype('float64').type(0.0)
80zone    = np.dtype('float64').type(1.0)
81epsfrac = np.dtype('float64').type(1.0E-10)
82pi      = np.pi
83rad     = pi/np.dtype('float64').type(180.0)  # Conversion from degrees to radian
84ra      = np.dtype('float64').type(6371229.0) # Earth radius
85
86## Functions
87def geodist (plon1, plat1, plon2, plat2) :
88    '''Distance between two points (on sphere)'''
89    zs = ( np.sin (rad*plat1) * np.sin (rad*plat2) +
90           np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) )
91    zs = np.maximum (-zone, np.minimum (zone, zs))
92    zgeodist =  np.arccos (zs)
93    return zgeodist
94
95### ===== Reading command line parameters ==================================================
96# Creating a parser
97parser = argparse.ArgumentParser (
98    description = '''Compute calving weights''',
99    epilog='-------- End of the help message --------')
100
101# Adding arguments
102parser.add_argument ('--oce'          , help='oce model name', type=str, default='eORCA1.2',
103                     choices=['ORCA2.3',  'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2',
104                              'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] )
105parser.add_argument ('--atm'          , help='atm model name', type=str, default='LMD9695'    )
106parser.add_argument ('--atmCoastWidth', type=int, default=1,
107                         help='width of the coastal band in atmosphere (in grid points)' )
108parser.add_argument ('--oceCoastWidth',
109                         help='width of the coastal band in ocean (in grid points)'     ,
110                         type=int, default=2 )
111parser.add_argument ('--atmQuantity'  , type=str, default='Quantity',
112                         choices=['Quantity', 'Surfacic'],
113                         help='Quantity if atm provides quantities (m/s), Surfacic if atm provided flux (m/s/m2)' )
114parser.add_argument ('--oceQuantity'  , type=str, default='Surfacic',
115                         choices=['Quantity', 'Surfacic'],
116                         help='Quantity if oce requires quantities (ks/s), Surfacic if oce requires flux (m/s/m2)' )
117parser.add_argument ('--searchRadius' , type=float, default=600.0,
118                         help='max distance to connect a land point to an ocean point (in km)' )
119parser.add_argument ('--grids' , help='grids file name', default='grids.nc' )
120parser.add_argument ('--areas' , help='masks file name', default='areas.nc' )
121parser.add_argument ('--masks' , help='areas file name', default='masks.nc' )
122parser.add_argument ('--o2a'   , help='o2a file name'  , default='o2a.nc'   )
123parser.add_argument ('--output', help='output rmp file name',
124                         default='rmp_tlmd_to_torc_runoff.nc' )
125parser.add_argument ('--fmt'   , default='netcdf4',
126                         help='NetCDF file format, using nco syntax',
127                         choices=['classic', 'netcdf3', '64bit', '64bit_data',
128                                  '64bit_data', 'netcdf4', 'netcdf4_classsic'] )
129parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=float,
130                         default=0, choices=nemo.NPERIO_VALID_RANGE)
131parser.add_argument ('--modelName' , help='Name of model', type=str, default=None)
132
133# Parse command line
134myargs = parser.parse_args()
135
136#
137grids = myargs.grids
138areas = myargs.areas
139masks = myargs.masks
140o2a   = myargs.o2a
141
142# Model Names
143atm_Name   = myargs.atm
144oce_Name   = myargs.oce
145model_name = myargs.modelName
146# Width of the coastal band (land points) in the atmopshere
147atmCoastWidth = myargs.atmCoastWidth
148# Width of the coastal band (ocean points) in the ocean
149oceCoastWidth = myargs.oceCoastWidth
150searchRadius  = myargs.searchRadius * 1000.0 # From km to meters
151# Netcdf format
152if myargs.fmt == 'classic'         : FmtNetcdf = 'CLASSIC'
153if myargs.fmt == 'netcdf3'         : FmtNetcdf = 'CLASSIC'
154if myargs.fmt == '64bit'           : FmtNetcdf = 'NETCDF3_64BIT_OFFSET'
155if myargs.fmt == '64bit_data'      : FmtNetcdf = 'NETCDF3_64BIT_DATA'
156if myargs.fmt == '64bit_offset'    : FmtNetcdf = 'NETCDF3_64BIT_OFFSET'
157if myargs.fmt == 'netcdf4'         : FmtNetcdf = 'NETCDF4'
158if myargs.fmt == 'netcdf4_classic' : FmtNetcdf = 'NETCDF4_CLASSIC'
159
160#
161if atm_Name.find('LMD') >= 0 : atm_n = 'lmd' ; atmDomainType = 'rectilinear'
162if atm_Name.find('ICO') >= 0 : atm_n = 'ico' ; atmDomainType = 'unstructured'
163
164print ('atmQuantity : ' + str (myargs.atmQuantity) )
165print ('oceQuantity : ' + str (myargs.oceQuantity) )
166
167# Ocean grid periodicity
168oce_perio = myargs.ocePerio
169
170### Read coordinates of all models
171###
172
173diaFile    = xr.open_dataset ( o2a   )
174gridFile   = xr.open_dataset ( grids )
175areaFile   = xr.open_dataset ( areas )
176maskFile   = xr.open_dataset ( masks )
177
178o2aFrac             = diaFile ['OceFrac'].squeeze()
179o2aFrac = np.where ( np.abs(o2aFrac) < 1E10, o2aFrac, 0.0)
180
181(atm_nvertex, atm_jpj, atm_jpi) = gridFile['t'+atm_n+'.clo'][:].shape
182atm_grid_size = atm_jpj*atm_jpi
183atm_grid_rank = len(gridFile['t'+atm_n+'.lat'][:].shape)
184
185atm_grid_center_lat = gridFile['t'+atm_n+'.lat'].squeeze()
186atm_grid_center_lon = gridFile['t'+atm_n+'.lon'].squeeze()
187atm_grid_corner_lat = gridFile['t'+atm_n+'.cla'].squeeze()
188atm_grid_corner_lon = gridFile['t'+atm_n+'.clo'].squeeze()
189atm_grid_area       = areaFile['t'+atm_n+'.srf'].squeeze()
190atm_grid_imask      = 1-maskFile['t'+atm_n+'.msk'][:].squeeze()
191atm_grid_dims       = gridFile['t'+atm_n+'.lat'][:].shape
192
193if atmDomainType == 'unstructured' :
194    atm_grid_center_lat = atm_grid_center_lat.rename ({'ycell':'cell'})
195    atm_grid_center_lon = atm_grid_center_lon.rename ({'ycell':'cell'})
196    atm_grid_corner_lat = atm_grid_corner_lat.rename ({'ycell':'cell'})
197    atm_grid_corner_lon = atm_grid_corner_lon.rename ({'ycell':'cell'})
198    atm_grid_area       = atm_grid_area.rename  ({'ycell':'cell'})
199    atm_grid_imask      = atm_grid_imask.rename ({'ycell':'cell'})
200
201if atmDomainType == 'rectilinear' :
202    atm_grid_center_lat = atm_grid_center_lat.stack (cell=['y', 'x'])
203    atm_grid_center_lon = atm_grid_center_lon.stack (cell=['y', 'x'])
204    atm_grid_corner_lat = atm_grid_corner_lat.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'})
205    atm_grid_corner_lon = atm_grid_corner_lon.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'})
206    atm_grid_area       = atm_grid_area.stack       (cell=['y', 'x'])
207    atm_grid_imask      = atm_grid_imask.stack      (cell=['y', 'x'])
208
209atm_perio = 0
210atm_grid_pmask = atm_grid_imask
211
212atm_address = np.arange(atm_jpj*atm_jpi)
213
214(oce_nvertex, oce_jpj, oce_jpi) = gridFile['torc.cla'][:].shape ; jpon=oce_jpj*oce_jpj
215oce_grid_size = oce_jpj*oce_jpi
216oce_grid_rank = len(gridFile['torc.lat'][:].shape)
217
218oce_grid_center_lat = gridFile['torc.lat'].stack(oce_grid_size=['y_grid_T', 'x_grid_T'])
219oce_grid_center_lon = gridFile['torc.lon'].stack(oce_grid_size=['y_grid_T', 'x_grid_T'])
220oce_grid_corner_lat = gridFile['torc.cla'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T'])
221oce_grid_corner_lon = gridFile['torc.clo'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T'])
222oce_grid_area       = areaFile['torc.srf'].stack(oce_grid_size=['y_grid_T', 'x_grid_T'])
223oce_grid_imask      = 1-maskFile['torc.msk'].stack(oce_grid_size=['y_grid_T', 'x_grid_T'])
224oce_grid_dims       = gridFile['torc.lat'][:].shape
225
226# imask : includes periodicity point
227# pmask : exclude periodicity points
228
229#if oce_perio == 0 :
230#    if oce_jpi ==  182 : oce_perio = 4 # ORCA 2
231#    if oce_jpi ==  362 : oce_perio = 6 # ORCA 1
232#    if oce_jpi == 1442 : oce_perio = 6 # ORCA 025
233oce_preio = nemo.__guess_nperio__ (oce_jpj, oce_jpi, nperio=oce_perio)
234
235print ( f'Oce NPERIO parameter : {oce_perio}' )
236oce_grid_imask = nemo.lbc      (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)),
237                                nperio=oce_perio, cd_type='T'        ).ravel()
238oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask       , (oce_jpj,oce_jpi)),
239                                nperio=oce_perio, cd_type='T', sval=0).ravel()
240oce_address = np.arange(oce_jpj*oce_jpi)
241
242print ('Fill closed seas with image processing library')
243oce_grid_imask2D = np.reshape (oce_grid_imask, (oce_jpj,oce_jpi))
244oce_grid_imask2D = 1-nemo.fill_closed_seas (1-oce_grid_imask2D, nperio=oce_perio, cd_type='T')
245#oce_grid_imask2D = nemo.lbc ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')),
246#                                       nperio=oce_perio, cd_type='T' )
247oce_grid_imask = oce_grid_imask2D.ravel ()
248
249##
250print ('Computes an ocean coastal band with image processing library')
251oceLand2D  = np.reshape ( np.where (oce_grid_imask < 0.5, True, False), (oce_jpj, oce_jpi) )
252oceOcean2D = np.reshape ( np.where (oce_grid_imask > 0.5, True, False), (oce_jpj, oce_jpi) )
253
254NNocean = 1+2*oceCoastWidth
255oceOceanFiltered2D = ndimage.uniform_filter (oceOcean2D.astype(float), size=NNocean)
256oceCoast2D = np.where (oceOceanFiltered2D<(1.0-0.5/(NNocean**2)),True,False) & oceOcean2D
257oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D, (oce_jpj,oce_jpi)),
258                            nperio=oce_perio, cd_type='T', sval=0.).ravel()
259
260oceOceanFiltered = oceOceanFiltered2D.ravel()
261oceLand  = oceLand2D.ravel ()
262oceOcean = oceOcean2D.ravel()
263oceCoast = oceCoast2D.ravel()
264
265print ( f'Number of points in oceLand  : {oceLand.sum():8d}' )
266print ( f'Number of points in oceOcean : {oceOcean.sum():8d}' )
267print ( f'Number of points in oceCoast : {oceCoast.sum():8d}' )
268
269# Arrays with coastal points only
270oceCoast_grid_center_lon = oce_grid_center_lon[oceCoast]
271oceCoast_grid_center_lat = oce_grid_center_lat[oceCoast]
272oceCoast_grid_area       = oce_grid_area      [oceCoast]
273oceCoast_grid_imask      = oce_grid_imask     [oceCoast]
274oceCoast_grid_pmask      = oce_grid_pmask     [oceCoast]
275oceCoast_address         = oce_address        [oceCoast]
276
277print ('Computes an atmosphere coastal band' )
278atmLand      = np.where (o2aFrac[:] < epsfrac       , True, False)
279atmLandFrac  = np.where (o2aFrac[:] < zone-epsfrac  , True, False)
280atmFrac      = np.where (o2aFrac[:] > epsfrac       , True, False) & \
281                         np.where (o2aFrac[:] < (zone-epsfrac), True, False)
282atmOcean     = np.where (o2aFrac[:] < (zone-epsfrac), True, False)
283atmOceanFrac = np.where (o2aFrac[:] > epsfrac       , True, False)
284
285## For LMDZ only !!
286if atmDomainType == 'rectilinear' :
287    print ('Extend coastal band ' )
288    NNatm = 1+2*atmCoastWidth
289    atmLand2D = np.reshape ( atmLand, ( atm_jpj, atm_jpi) )
290
291    atmLandFiltered2D = ndimage.uniform_filter(atmLand2D.astype(float), size=NNatm)
292    atmCoast2D = np.where (atmLandFiltered2D<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac
293
294    atmLandFiltered = atmLandFiltered2D.ravel()
295    atmCoast = atmCoast2D.ravel()
296
297    print ( f'Number of points in atmLand  : {atmLand.sum():8d}'  )
298    print ( f'Number of points in atmOcean : {atmOcean.sum():8d}' )
299    print ( f'Number of points in atmCoast : {atmCoast.sum():8d}' )
300
301else :
302    atmCoast = atmFrac
303
304
305# Arrays with coastal points only
306atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast]
307atmCoast_grid_center_lat = atm_grid_center_lat[atmCoast]
308atmCoast_grid_area       = atm_grid_area      [atmCoast]
309atmCoast_grid_imask      = atm_grid_imask     [atmCoast]
310atmCoast_grid_pmask      = atm_grid_pmask     [atmCoast]
311atmCoast_address         = atm_address        [atmCoast]
312
313# Initialisations before the loop
314remap_matrix = np.empty ( shape=(0), dtype=np.float64 )
315atm_address  = np.empty ( shape=(0), dtype=np.int32   )
316oce_address  = np.empty ( shape=(0), dtype=np.int32   )
317
318## Loop on atmosphere coastal points
319if searchRadius > 0. :
320    print ('Loop on atmosphere coastal points')
321    for ja in np.arange (len(atmCoast_grid_pmask)) :
322        z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja],
323                               oceCoast_grid_center_lon, oceCoast_grid_center_lat)
324        z_mask = np.where (z_dist*ra < searchRadius, True, False)
325        num_links = int (z_mask.sum())
326        if num_links == 0 : continue
327        z_area = oceCoast_grid_area[z_mask].sum().values
328        poids = np.ones ((num_links),dtype=np.float64) / z_area
329        if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja]
330        if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask]
331        if  ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print
332            print ( f'ja:{ja:8d} num_links:{num_links:8d} z_area:{z_area:8.4e} '+
333                    f'atm area:{atm_grid_area[ja].values:8.4e} weights sum:{poids.sum():8.4e}' )
334        #
335        matrix_local = poids
336        atm_address_local = np.ones (num_links, dtype=np.int32 ) * atmCoast_address[ja]
337        # Address on destination grid
338        oce_address_local = oceCoast_address [z_mask]
339        # Append to global arrays
340        remap_matrix = np.append ( remap_matrix, matrix_local      )
341        atm_address  = np.append ( atm_address , atm_address_local )
342        oce_address  = np.append ( oce_address , oce_address_local )
343
344    print ('End of loop')
345
346num_links = remap_matrix.shape[0]
347
348print ('Write output file')
349runoff = myargs.output
350print ('Output file: ' + runoff )
351
352remap_matrix = xr.DataArray ( np.reshape(remap_matrix, (num_links, 1)),
353                                  dims = ['num_links', 'num_wgts'] )
354
355# OASIS uses Fortran style indexing, starting at one
356src_address  = xr.DataArray ( atm_address.astype(np.int32)+1, dims = ['num_links'],
357                                 attrs={'convention': 'Fortran style addressing, starting at 1'})
358dst_address  = xr.DataArray ( oce_address.astype(np.int32)+1, dims = ['num_links'],
359                                 attrs={'convention': 'Fortran style addressing, starting at 1'})
360
361src_grid_dims       = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32),
362                                        dims = ['src_grid_rank',] )
363src_grid_center_lon = xr.DataArray (atm_grid_center_lon.values , dims = ['src_grid_size',] )
364src_grid_center_lat = xr.DataArray (atm_grid_center_lat.values , dims = ['src_grid_size',] )
365
366src_grid_center_lon.attrs.update ( {'units':'degrees_east' , 'long_name':'Longitude',
367                                        'standard_name':'longitude',
368                                        'bounds':'src_grid_corner_lon'} )
369src_grid_center_lat.attrs.update ( {'units':'degrees_north', 'long_name':'Latitude' ,
370                                        'standard_name':'latitude' ,
371                                        'bounds':'src_grid_corner_lat'} )
372
373src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(),
374                                        dims = [ 'src_grid_size', 'src_grid_corners'] )
375src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(),
376                                        dims = [ 'src_grid_size', 'src_grid_corners'] )
377src_grid_corner_lon.attrs['units']='degrees_east'
378src_grid_corner_lat.attrs['units']='degrees_north'
379
380src_grid_area       =  xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] )
381src_grid_area.attrs.update ( {'long_name':'Grid area', 'standard_name':'cell_area',
382                                  'units':'m2' } )
383src_grid_imask      =  xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] )
384src_grid_imask.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0'} )
385
386src_grid_pmask      =  xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] )
387src_grid_pmask.attrs.update ( {'long_name':'Land-sea mask (periodicity removed)',
388                                   'units':'Land:1, Ocean:0' } )
389
390# --
391dst_grid_dims       = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32),
392                                        dims = ['dst_grid_rank',] )
393dst_grid_center_lon = xr.DataArray (oce_grid_center_lon.values, dims = ['dst_grid_size',] )
394dst_grid_center_lat = xr.DataArray (oce_grid_center_lat.values, dims = ['dst_grid_size',] )
395
396dst_grid_center_lon.attrs.update ( {'units':'degrees_east' , 'long_name':'Longitude',
397                                        'standard_name':'longitude',
398                                        'bounds':'dst_grid_corner_lon' } )
399dst_grid_center_lat.attrs.update ( {'units':'degrees_north', 'long_name':'Latitude' ,
400                                        'standard_name':'latitude' ,
401                                        'bounds':'dst_grid_corner_lat' } )
402
403dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values),
404                                        dims = [ 'dst_grid_size', 'dst_grid_corners'] )
405dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values),
406                                        dims = [ 'dst_grid_size', 'dst_grid_corners'] )
407dst_grid_corner_lon.attrs['units']='degrees_east'
408dst_grid_corner_lat.attrs['units']='degrees_north'
409
410dst_grid_area       =  xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] )
411dst_grid_area.attrs.update ( {'long_name':'Grid area', 'standard_name':'cell_area', 'units':'m2' })
412
413dst_grid_imask      =  xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] )
414dst_grid_imask.attrs.update ( {'long_name':'Land-sea mask', 'units':'Land:1, Ocean:0'} )
415
416dst_grid_pmask      =  xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] )
417dst_grid_pmask.attrs.update ( {'long_name':'Land-sea mask (periodicity removed)',
418                                   'units':'Land:1, Ocean:0'} )
419
420src_lon_addressed   =  xr.DataArray (atm_grid_center_lon.values[atm_address] ,
421                                         dims = ['num_links'] )
422src_lat_addressed   =  xr.DataArray (atm_grid_center_lat.values[atm_address] ,
423                                         dims = ['num_links'] )
424src_area_addressed  =  xr.DataArray (atm_grid_area      .values[atm_address] ,
425                                         dims = ['num_links'] )
426src_imask_addressed =  xr.DataArray (1-atm_grid_imask   .values[atm_address].astype(np.int32) ,
427                                         dims = ['num_links'] )
428src_pmask_addressed =  xr.DataArray (1-atm_grid_pmask   .values[atm_address].astype(np.int32) ,
429                                         dims = ['num_links'] )
430
431dst_lon_addressed   =  xr.DataArray (oce_grid_center_lon.values[atm_address],
432                                         dims = ['num_links'] )
433dst_lat_addressed   =  xr.DataArray (oce_grid_center_lat.values[oce_address],
434                                         dims = ['num_links'] )
435dst_area_addressed  =  xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32) ,
436                                         dims = ['num_links'] )
437dst_imask_addressed =  xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32)   ,
438                                         dims = ['num_links'] )
439dst_pmask_addressed =  xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32)   ,
440                                         dims = ['num_links'] )
441
442src_lon_addressed.attrs.update ({'units':'degrees_east' , 'long_name':'Longitude',
443                                     'standard_name':'longitude' })
444src_lat_addressed.attrs.update ({'units':'degrees_north', 'long_name':'Latitude' ,
445                                     'standard_name':'latitude'  })
446
447dst_lon_addressed.attrs.update ({'units':'degrees_east' , 'long_name':'Longitude',
448                                     'standard_name':'longitude' })
449dst_lat_addressed.attrs.update ({'units':'degrees_north', 'long_name':'Latitude' ,
450                                     'standard_name':'latitude'  })
451
452if atmDomainType == 'rectilinear' :
453    atmLand         = xr.DataArray ( atmLand.ravel()     , dims = ['src_grid_size',] )
454    atmLandFiltered = xr.DataArray ( atmLandFrac.ravel() , dims = ['src_grid_size',] )
455    atmLandFrac     = xr.DataArray ( atmFrac.ravel()     , dims = ['src_grid_size',] )
456    atmFrac         = xr.DataArray ( atmFrac.ravel()     , dims = ['src_grid_size',] )
457    atmOcean        = xr.DataArray ( atmOcean.ravel()    , dims = ['src_grid_size',] )
458    atmOceanFrac    = xr.DataArray ( atmOceanFrac.ravel(), dims = ['src_grid_size',] )
459
460atmCoast         = xr.DataArray (atmCoast.astype(np.int32)          , dims = ['src_grid_size',])
461oceLand          = xr.DataArray (oceLand.astype(np.int32)           , dims = ['dst_grid_size',])
462oceOcean         = xr.DataArray (oceOcean.astype(np.int32)          , dims = ['dst_grid_size',])
463oceOceanFiltered = xr.DataArray (oceOceanFiltered.astype(np.float32), dims = ['dst_grid_size',])
464oceCoast         = xr.DataArray (oceCoast.astype(np.int32)          , dims = ['dst_grid_size',])
465
466
467f_runoff = xr.Dataset ( {
468    'remap_matrix'        : remap_matrix,
469        'src_address'         : src_address,
470        'dst_address'         : dst_address,
471        'src_grid_dims'       : src_grid_dims,
472        'src_grid_center_lon' : src_grid_center_lon,
473        'src_grid_center_lat' : src_grid_center_lat,
474    'src_grid_corner_lon' : src_grid_corner_lon,
475    'src_grid_corner_lat' : src_grid_corner_lat,
476        'src_grid_area'       : src_grid_area,
477        'src_grid_imask'      : src_grid_imask,
478        'src_grid_pmask'      : src_grid_pmask,
479        'dst_grid_dims'       : dst_grid_dims,
480        'dst_grid_center_lon' : dst_grid_center_lon,
481        'dst_grid_center_lat' : dst_grid_center_lat,
482        'dst_grid_corner_lon' : dst_grid_corner_lon,
483        'dst_grid_corner_lat' : dst_grid_corner_lat,
484        'dst_grid_area'       : dst_grid_area,
485        'dst_grid_imask'      : dst_grid_imask,
486        'dst_grid_pmask'      : dst_grid_pmask,
487        'src_lon_addressed'   : src_lon_addressed,
488        'src_lat_addressed'   : src_lat_addressed,
489        'src_area_addressed'  : src_area_addressed,
490        'dst_lon_addressed'   : dst_lon_addressed,
491        'dst_lat_addressed'   : dst_lat_addressed,
492        'dst_area_addressed'  : dst_area_addressed,
493        'dst_imask_addressed' : dst_imask_addressed,
494        'dst_pmask_addressed' : dst_pmask_addressed,
495        'atmCoast'            : atmCoast,
496        'oceLand'             : oceLand,
497        'oceOcean'            : oceOcean,
498        'oceOceanFiltered'    : oceOceanFiltered,
499        'oceCoast'            : oceCoast
500    } )
501
502f_runoff.attrs.update ( {
503    'Conventions'     : 'CF-1.6',
504    'source'          : 'IPSL Earth system model',
505    'group'           : 'ICMC IPSL Climate Modelling Center',
506    'Institution'     : 'IPSL https.//www.ipsl.fr',
507    'Ocean'           : f'{oce_Name} https://www.nemo-ocean.eu',
508    'Atmosphere'      : f'{atm_Name} http://lmdz.lmd.jussieu.fr',
509    'associatedFiles' : f'{grids} {areas} {masks}',
510    'description'     : 'Generated with RunoffWeights.py',
511    'title'           : runoff,
512    'Program'         : f'Generated by {sys.argv[0]} with flags ' + ' '.join (sys.argv[1:]),
513    'atmCoastWidth'   : f'{atmCoastWidth:d} grid points',
514    'oceCoastWidth'   : f'{oceCoastWidth:d} grid points',
515    'searchRadius'    : f'{searchRadius/1000.:.0f} km',
516    'atmQuantity'     : myargs.atmQuantity,
517    'oceQuantity'     : myargs.oceQuantity,
518    'gridsFile'       : grids,
519    'areasFile'       : areas,
520    'masksFile'       : masks,
521    'o2aFile'         : o2a,
522    'timeStamp'       : time.asctime (),
523    'directory'       : os.getcwd (),
524    'HOSTNAME'        : platform.node (),
525    'LOGNAME'         : getpass.getuser(),
526    'Python'          : f'Python version {platform.python_version ()}',
527    'OS'              : platform.system (),
528    'release'         : platform.release (),
529    'hardware'        : platform.machine (),
530    'conventions'     : 'SCRIP',
531    'source_grid'     : 'curvilinear',
532    'dest_grid'       : 'curvilinear',
533    'Model'           : model_name,
534    'SVN_Author'      : '$Author$',
535    'SVN_Date'        : '$Date$',
536    'SVN_Revision'    : '$Revision$',
537    'SVN_Id'          : '$Id$',
538    'SVN_HeadURL'     : '$HeadURL$',
539    } )
540
541f_runoff.to_netcdf ( runoff, mode='w', format=FmtNetcdf )
542f_runoff.close ()
543
544##
545
546print ("That''s all folks !")
547## ======================================================================================
Note: See TracBrowser for help on using the repository browser.