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 | # SVN information |
---|
32 | __Author__ = "$Author$" |
---|
33 | __Date__ = "$Date$" |
---|
34 | __Revision__ = "$Revision$" |
---|
35 | __Id__ = "$Id$" |
---|
36 | __HeadURL__ = "$HeadURL$" |
---|
37 | __SVN_Date__ = "$SVN_Date: $" |
---|
38 | ## |
---|
39 | |
---|
40 | ## Modules |
---|
41 | import numpy as np, xarray as xr |
---|
42 | import nemo |
---|
43 | from scipy import ndimage |
---|
44 | import sys, os, platform, argparse, textwrap, time |
---|
45 | |
---|
46 | ## Useful constants |
---|
47 | zero = np.dtype('float64').type(0.0) |
---|
48 | zone = np.dtype('float64').type(1.0) |
---|
49 | epsfrac = np.dtype('float64').type(1.0E-10) |
---|
50 | pi = np.pi |
---|
51 | rad = pi/np.dtype('float64').type(180.0) # Conversion from degrees to radian |
---|
52 | ra = np.dtype('float64').type(6371229.0) # Earth radius |
---|
53 | |
---|
54 | ## Functions |
---|
55 | def geodist (plon1, plat1, plon2, plat2) : |
---|
56 | """Distance between two points (on sphere)""" |
---|
57 | zs = np.sin (rad*plat1) * np.sin (rad*plat2) + np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) |
---|
58 | zs = np.maximum (-zone, np.minimum (zone, zs)) |
---|
59 | geodist = np.arccos (zs) |
---|
60 | return geodist |
---|
61 | |
---|
62 | ### ===== Reading command line parameters ================================================== |
---|
63 | # Creating a parser |
---|
64 | parser = argparse.ArgumentParser ( |
---|
65 | description = """Compute calving weights""", |
---|
66 | epilog='-------- End of the help message --------') |
---|
67 | |
---|
68 | # Adding arguments |
---|
69 | parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', |
---|
70 | choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) |
---|
71 | parser.add_argument ('--atm' , help='atm model name', type=str, default='LMD9695' ) |
---|
72 | parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', type=int, default=1 ) |
---|
73 | parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)' , type=int, default=2 ) |
---|
74 | parser.add_argument ('--atmQuantity' , help='Quantity if atm provides quantities (m/s), Surfacic if atm provided flux (m/s/m2)' , type=str, |
---|
75 | default='Quantity', choices=['Quantity', 'Surfacic'] ) |
---|
76 | parser.add_argument ('--oceQuantity' , help='Quantity if oce requires quantities (ks/s), Surfacic if oce requires flux (m/s/m2)', type=str, |
---|
77 | default='Surfacic', choices=['Quantity', 'Surfacic'] ) |
---|
78 | parser.add_argument ('--searchRadius' , help='max distance to connect a land point to an ocean point (in km)', type=float, default=600.0 ) |
---|
79 | parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) |
---|
80 | parser.add_argument ('--areas' , help='masks file name', default='areas.nc' ) |
---|
81 | parser.add_argument ('--masks' , help='areas file name', default='masks.nc' ) |
---|
82 | parser.add_argument ('--o2a' , help='o2a file name' , default='o2a.nc' ) |
---|
83 | parser.add_argument ('--output', help='output rmp file name', default='rmp_tlmd_to_torc_runoff.nc' ) |
---|
84 | parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', |
---|
85 | choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) |
---|
86 | parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=float, default=0, choices=nemo.nperio_valid_range) |
---|
87 | |
---|
88 | # Parse command line |
---|
89 | myargs = parser.parse_args() |
---|
90 | |
---|
91 | # |
---|
92 | grids = myargs.grids |
---|
93 | areas = myargs.areas |
---|
94 | masks = myargs.masks |
---|
95 | o2a = myargs.o2a |
---|
96 | |
---|
97 | # Model Names |
---|
98 | atm_Name = myargs.atm |
---|
99 | oce_Name = myargs.oce |
---|
100 | # Width of the coastal band (land points) in the atmopshere |
---|
101 | atmCoastWidth = myargs.atmCoastWidth |
---|
102 | # Width of the coastal band (ocean points) in the ocean |
---|
103 | oceCoastWidth = myargs.oceCoastWidth |
---|
104 | searchRadius = myargs.searchRadius * 1000.0 # From km to meters |
---|
105 | # Netcdf format |
---|
106 | if myargs.fmt == 'classic' : FmtNetcdf = 'CLASSIC' |
---|
107 | if myargs.fmt == 'netcdf3' : FmtNetcdf = 'CLASSIC' |
---|
108 | if myargs.fmt == '64bit' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' |
---|
109 | if myargs.fmt == '64bit_data' : FmtNetcdf = 'NETCDF3_64BIT_DATA' |
---|
110 | if myargs.fmt == '64bit_offset' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' |
---|
111 | if myargs.fmt == 'netcdf4' : FmtNetcdf = 'NETCDF4' |
---|
112 | if myargs.fmt == 'netcdf4_classic' : FmtNetcdf = 'NETCDF4_CLASSIC' |
---|
113 | |
---|
114 | # |
---|
115 | if atm_Name.find('LMD') >= 0 : atm_n = 'lmd' ; atmDomainType = 'rectilinear' |
---|
116 | if atm_Name.find('ICO') >= 0 : atm_n = 'ico' ; atmDomainType = 'unstructured' |
---|
117 | |
---|
118 | print ('atmQuantity : ' + str (myargs.atmQuantity) ) |
---|
119 | print ('oceQuantity : ' + str (myargs.oceQuantity) ) |
---|
120 | |
---|
121 | # Ocean grid periodicity |
---|
122 | oce_perio = myargs.ocePerio |
---|
123 | |
---|
124 | ### Read coordinates of all models |
---|
125 | ### |
---|
126 | |
---|
127 | diaFile = xr.open_dataset ( o2a ) |
---|
128 | gridFile = xr.open_dataset ( grids ) |
---|
129 | areaFile = xr.open_dataset ( areas ) |
---|
130 | maskFile = xr.open_dataset ( masks ) |
---|
131 | |
---|
132 | o2aFrac = diaFile ['OceFrac'].squeeze() |
---|
133 | o2aFrac = np.where ( np.abs(o2aFrac) < 1E10, o2aFrac, 0.0) |
---|
134 | |
---|
135 | (atm_nvertex, atm_jpj, atm_jpi) = gridFile['t'+atm_n+'.clo'][:].shape |
---|
136 | atm_grid_size = atm_jpj*atm_jpi |
---|
137 | atm_grid_rank = len(gridFile['t'+atm_n+'.lat'][:].shape) |
---|
138 | |
---|
139 | atm_grid_center_lat = gridFile['t'+atm_n+'.lat'].squeeze() |
---|
140 | atm_grid_center_lon = gridFile['t'+atm_n+'.lon'].squeeze() |
---|
141 | atm_grid_corner_lat = gridFile['t'+atm_n+'.cla'].squeeze() |
---|
142 | atm_grid_corner_lon = gridFile['t'+atm_n+'.clo'].squeeze() |
---|
143 | atm_grid_area = areaFile['t'+atm_n+'.srf'].squeeze() |
---|
144 | atm_grid_imask = 1-maskFile['t'+atm_n+'.msk'][:].squeeze() |
---|
145 | atm_grid_dims = gridFile['t'+atm_n+'.lat'][:].shape |
---|
146 | |
---|
147 | if atmDomainType == 'unstructured' : |
---|
148 | atm_grid_center_lat = atm_grid_center_lat.rename ({'ycell':'cell'}) |
---|
149 | atm_grid_center_lon = atm_grid_center_lon.rename ({'ycell':'cell'}) |
---|
150 | atm_grid_corner_lat = atm_grid_corner_lat.rename ({'ycell':'cell'}) |
---|
151 | atm_grid_corner_lon = atm_grid_corner_lon.rename ({'ycell':'cell'}) |
---|
152 | atm_grid_area = atm_grid_area.rename ({'ycell':'cell'}) |
---|
153 | atm_grid_imask = atm_grid_imask.rename ({'ycell':'cell'}) |
---|
154 | |
---|
155 | if atmDomainType == 'rectilinear' : |
---|
156 | atm_grid_center_lat = atm_grid_center_lat.stack (cell=['y', 'x']) |
---|
157 | atm_grid_center_lon = atm_grid_center_lon.stack (cell=['y', 'x']) |
---|
158 | atm_grid_corner_lat = atm_grid_corner_lat.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'}) |
---|
159 | atm_grid_corner_lon = atm_grid_corner_lon.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'}) |
---|
160 | atm_grid_area = atm_grid_area.stack (cell=['y', 'x']) |
---|
161 | atm_grid_imask = atm_grid_imask.stack (cell=['y', 'x']) |
---|
162 | |
---|
163 | atm_perio = 0 |
---|
164 | atm_grid_pmask = atm_grid_imask |
---|
165 | atm_address = np.arange(atm_jpj*atm_jpi) |
---|
166 | |
---|
167 | (oce_nvertex, oce_jpj, oce_jpi) = gridFile['torc.cla'][:].shape ; jpon=oce_jpj*oce_jpj |
---|
168 | oce_grid_size = oce_jpj*oce_jpi |
---|
169 | oce_grid_rank = len(gridFile['torc.lat'][:].shape) |
---|
170 | |
---|
171 | oce_grid_center_lat = gridFile['torc.lat'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
172 | oce_grid_center_lon = gridFile['torc.lon'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
173 | oce_grid_corner_lat = gridFile['torc.cla'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
174 | oce_grid_corner_lon = gridFile['torc.clo'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
175 | oce_grid_area = areaFile['torc.srf'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
176 | oce_grid_imask = 1-maskFile['torc.msk'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
177 | oce_grid_dims = gridFile['torc.lat'][:].shape |
---|
178 | |
---|
179 | if oce_perio == 0 : |
---|
180 | if oce_jpi == 182 : oce_perio = 4 # ORCA 2 |
---|
181 | if oce_jpi == 362 : oce_perio = 6 # ORCA 1 |
---|
182 | if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 |
---|
183 | |
---|
184 | print ("Oce NPERIO parameter : {:}".format(oce_perio)) |
---|
185 | oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel() |
---|
186 | oce_address = np.arange(oce_jpj*oce_jpi) |
---|
187 | |
---|
188 | print ("Fill closed sea with image processing library") |
---|
189 | oce_grid_imask2D = np.reshape(oce_grid_pmask,(oce_jpj,oce_jpi)) |
---|
190 | oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), |
---|
191 | nperio=oce_perio, cd_type='T', sval=0 ) |
---|
192 | oce_grid_imask = oce_grid_imask2D.ravel() |
---|
193 | ## |
---|
194 | print ("Computes an ocean coastal band") |
---|
195 | |
---|
196 | oceLand2D = np.reshape ( np.where (oce_grid_pmask < 0.5, True, False), (oce_jpj, oce_jpi) ) |
---|
197 | oceOcean2D = np.reshape ( np.where (oce_grid_pmask > 0.5, True, False), (oce_jpj, oce_jpi) ) |
---|
198 | |
---|
199 | NNocean = 1+2*oceCoastWidth |
---|
200 | oceOceanFiltered2D = ndimage.uniform_filter(oceOcean2D.astype(float), size=NNocean) |
---|
201 | oceCoast2D = np.where (oceOceanFiltered2D<(1.0-0.5/(NNocean**2)),True,False) & oceOcean2D |
---|
202 | oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D,(oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T').ravel() |
---|
203 | |
---|
204 | oceOceanFiltered = oceOceanFiltered2D.ravel() |
---|
205 | oceLand = oceLand2D.ravel () |
---|
206 | oceOcean = oceOcean2D.ravel() |
---|
207 | oceCoast = oceCoast2D.ravel() |
---|
208 | |
---|
209 | print ('Number of points in oceLand : {:8d}'.format (oceLand.sum()) ) |
---|
210 | print ('Number of points in oceOcean : {:8d}'.format (oceOcean.sum()) ) |
---|
211 | print ('Number of points in oceCoast : {:8d}'.format (oceCoast.sum()) ) |
---|
212 | |
---|
213 | # Arrays with coastal points only |
---|
214 | oceCoast_grid_center_lon = oce_grid_center_lon[oceCoast] |
---|
215 | oceCoast_grid_center_lat = oce_grid_center_lat[oceCoast] |
---|
216 | oceCoast_grid_area = oce_grid_area [oceCoast] |
---|
217 | oceCoast_grid_imask = oce_grid_imask [oceCoast] |
---|
218 | oceCoast_grid_pmask = oce_grid_pmask [oceCoast] |
---|
219 | oceCoast_address = oce_address [oceCoast] |
---|
220 | |
---|
221 | print ("Computes an atmosphere coastal band " ) |
---|
222 | atmLand = np.where (o2aFrac[:] < epsfrac , True, False) |
---|
223 | atmLandFrac = np.where (o2aFrac[:] < zone-epsfrac , True, False) |
---|
224 | atmFrac = np.where (o2aFrac[:] > epsfrac , True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) |
---|
225 | atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False) |
---|
226 | atmOceanFrac = np.where (o2aFrac[:] > epsfrac , True, False) |
---|
227 | |
---|
228 | ## For LMDZ only !! |
---|
229 | if atmDomainType == 'rectilinear' : |
---|
230 | print ("Extend coastal band " ) |
---|
231 | NNatm = 1+2*atmCoastWidth |
---|
232 | atmLand2D = np.reshape ( atmLand, ( atm_jpj, atm_jpi) ) |
---|
233 | |
---|
234 | atmLandFiltered2D = ndimage.uniform_filter(atmLand2D.astype(float), size=NNatm) |
---|
235 | atmCoast2D = np.where (atmLandFiltered2D<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac |
---|
236 | |
---|
237 | atmLandFiltered = atmLandFiltered2D.ravel() |
---|
238 | atmCoast = atmCoast2D.ravel() |
---|
239 | |
---|
240 | print ('Number of points in atmLand : {:8d}'.format (atmLand.sum()) ) |
---|
241 | print ('Number of points in atmOcean : {:8d}'.format (atmOcean.sum()) ) |
---|
242 | print ('Number of points in atmCoast : {:8d}'.format (atmCoast.sum()) ) |
---|
243 | |
---|
244 | else : |
---|
245 | atmCoast = atmFrac |
---|
246 | |
---|
247 | |
---|
248 | # Arrays with coastal points only |
---|
249 | atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast] |
---|
250 | atmCoast_grid_center_lat = atm_grid_center_lat[atmCoast] |
---|
251 | atmCoast_grid_area = atm_grid_area [atmCoast] |
---|
252 | atmCoast_grid_imask = atm_grid_imask [atmCoast] |
---|
253 | atmCoast_grid_pmask = atm_grid_pmask [atmCoast] |
---|
254 | atmCoast_address = atm_address [atmCoast] |
---|
255 | |
---|
256 | # Initialisations before the loop |
---|
257 | remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) |
---|
258 | atm_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
259 | oce_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
260 | |
---|
261 | ## Loop on atmosphere coastal points |
---|
262 | if searchRadius > 0. : |
---|
263 | print ("Loop on atmosphere coastal points") |
---|
264 | for ja in np.arange(len(atmCoast_grid_pmask)) : |
---|
265 | z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) |
---|
266 | z_mask = np.where (z_dist*ra < searchRadius, True, False) |
---|
267 | num_links = int(z_mask.sum()) |
---|
268 | if num_links == 0 : continue |
---|
269 | z_area = oceCoast_grid_area[z_mask].sum().values |
---|
270 | poids = np.ones ((num_links),dtype=np.float64) / z_area |
---|
271 | if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja] |
---|
272 | if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] |
---|
273 | if ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print |
---|
274 | print ( 'ja:{:8d}, num_links:{:8d}, z_area:{:8.4e}, atm area:{:8.4e}, weights sum:{:8.4e} ' |
---|
275 | .format(ja, num_links, z_area, atm_grid_area[ja].values, poids.sum() ) ) |
---|
276 | # |
---|
277 | matrix_local = poids |
---|
278 | atm_address_local = np.ones(num_links, dtype=np.int32 ) * atmCoast_address[ja] |
---|
279 | # Address on destination grid |
---|
280 | oce_address_local = oceCoast_address[z_mask] |
---|
281 | # Append to global arrays |
---|
282 | remap_matrix = np.append ( remap_matrix, matrix_local ) |
---|
283 | atm_address = np.append ( atm_address , atm_address_local ) |
---|
284 | oce_address = np.append ( oce_address , oce_address_local ) |
---|
285 | |
---|
286 | print ('End of loop') |
---|
287 | |
---|
288 | num_links = remap_matrix.shape[0] |
---|
289 | |
---|
290 | print ("Write output file") |
---|
291 | runoff = myargs.output |
---|
292 | print ('Output file: ' + runoff ) |
---|
293 | |
---|
294 | |
---|
295 | remap_matrix = xr.DataArray ( np.reshape(remap_matrix, (num_links, 1)), dims = ['num_links', 'num_wgts'] ) |
---|
296 | |
---|
297 | # OASIS uses Fortran style indexing, starting at one |
---|
298 | src_address = xr.DataArray ( atm_address.astype(np.int32)+1, dims = ['num_links'], |
---|
299 | attrs={"convention": "Fortran style addressing, starting at 1"}) |
---|
300 | dst_address = xr.DataArray ( oce_address.astype(np.int32)+1, dims = ['num_links'], |
---|
301 | attrs={"convention": "Fortran style addressing, starting at 1"}) |
---|
302 | |
---|
303 | src_grid_dims = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32), dims = ['src_grid_rank',] ) |
---|
304 | src_grid_center_lon = xr.DataArray (atm_grid_center_lon.values , dims = ['src_grid_size',] ) |
---|
305 | src_grid_center_lat = xr.DataArray (atm_grid_center_lat.values , dims = ['src_grid_size',] ) |
---|
306 | |
---|
307 | src_grid_center_lon.attrs['units']='degrees_east' ; src_grid_center_lon.attrs['long_name']='Longitude' |
---|
308 | src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" |
---|
309 | src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' |
---|
310 | src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" |
---|
311 | |
---|
312 | src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) |
---|
313 | src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) |
---|
314 | src_grid_corner_lon.attrs['units']="degrees_east" |
---|
315 | src_grid_corner_lat.attrs['units']="degrees_north" |
---|
316 | |
---|
317 | src_grid_area = xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] ) |
---|
318 | src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" |
---|
319 | |
---|
320 | src_grid_imask = xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] ) |
---|
321 | src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" |
---|
322 | |
---|
323 | src_grid_pmask = xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] ) |
---|
324 | src_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; src_grid_pmask.attrs['units']="Land:1, Ocean:0" |
---|
325 | |
---|
326 | # -- |
---|
327 | dst_grid_dims = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32), dims = ['dst_grid_rank',] ) |
---|
328 | dst_grid_center_lon = xr.DataArray (oce_grid_center_lon.values, dims = ['dst_grid_size',] ) |
---|
329 | dst_grid_center_lat = xr.DataArray (oce_grid_center_lat.values, dims = ['dst_grid_size',] ) |
---|
330 | |
---|
331 | dst_grid_center_lon.attrs['units']='degrees_east' ; dst_grid_center_lon.attrs['long_name']='Longitude' |
---|
332 | dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" |
---|
333 | dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' |
---|
334 | dst_grid_center_lat.attrs['long_name']='latitude ' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" |
---|
335 | |
---|
336 | dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) |
---|
337 | dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) |
---|
338 | dst_grid_corner_lon.attrs['units']="degrees_east" |
---|
339 | dst_grid_corner_lat.attrs['units']="degrees_north" |
---|
340 | |
---|
341 | dst_grid_area = xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] ) |
---|
342 | dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" |
---|
343 | |
---|
344 | dst_grid_imask = xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] ) |
---|
345 | dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" |
---|
346 | |
---|
347 | dst_grid_pmask = xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] ) |
---|
348 | dst_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; dst_grid_pmask.attrs['units']="Land:1, Ocean:0" |
---|
349 | |
---|
350 | src_lon_addressed = xr.DataArray (atm_grid_center_lon.values[atm_address] , dims = ['num_links'] ) |
---|
351 | src_lat_addressed = xr.DataArray (atm_grid_center_lat.values[atm_address] , dims = ['num_links'] ) |
---|
352 | src_area_addressed = xr.DataArray (atm_grid_area .values[atm_address] , dims = ['num_links'] ) |
---|
353 | src_imask_addressed = xr.DataArray (1-atm_grid_imask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) |
---|
354 | src_pmask_addressed = xr.DataArray (1-atm_grid_pmask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) |
---|
355 | |
---|
356 | dst_lon_addressed = xr.DataArray (oce_grid_center_lon.values[atm_address], dims = ['num_links'] ) |
---|
357 | dst_lat_addressed = xr.DataArray (oce_grid_center_lat.values[oce_address], dims = ['num_links'] ) |
---|
358 | dst_area_addressed = xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
359 | dst_imask_addressed = xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
360 | dst_pmask_addressed = xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
361 | |
---|
362 | src_lon_addressed.attrs['long_name']="Longitude" ; src_lon_addressed.attrs['standard_name']="longitude" ; src_lon_addressed.attrs['units']="degrees_east" |
---|
363 | src_lat_addressed.attrs['long_name']="Latitude" ; src_lat_addressed.attrs['standard_name']="latitude" ; src_lat_addressed.attrs['units']="degrees_north" |
---|
364 | |
---|
365 | dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" |
---|
366 | dst_lat_addressed.attrs['long_name']="Latitude" ; dst_lat_addressed.attrs['standard_name']="latitude" ; dst_lat_addressed.attrs['units']="degrees_north" |
---|
367 | |
---|
368 | if atmDomainType == 'rectilinear' : |
---|
369 | atmLand = xr.DataArray ( atmLand.ravel() , dims = ['src_grid_size',] ) |
---|
370 | atmLandFiltered = xr.DataArray ( atmLandFrac.ravel() , dims = ['src_grid_size',] ) |
---|
371 | atmLandFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) |
---|
372 | atmFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) |
---|
373 | atmOcean = xr.DataArray ( atmOcean.ravel() , dims = ['src_grid_size',] ) |
---|
374 | atmOceanFrac = xr.DataArray ( atmOceanFrac.ravel(), dims = ['src_grid_size',] ) |
---|
375 | |
---|
376 | atmCoast = xr.DataArray (atmCoast.astype(np.int32) , dims = ['src_grid_size',]) |
---|
377 | oceLand = xr.DataArray (oceLand.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
378 | oceOcean = xr.DataArray (oceOcean.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
379 | oceOceanFiltered = xr.DataArray (oceOceanFiltered.astype(np.float32), dims = ['dst_grid_size',]) |
---|
380 | oceCoast = xr.DataArray (oceCoast.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
381 | |
---|
382 | |
---|
383 | f_runoff = xr.Dataset ( { |
---|
384 | 'remap_matrix' : remap_matrix, |
---|
385 | 'src_address' : src_address, |
---|
386 | 'dst_address' : dst_address, |
---|
387 | 'src_grid_dims' : src_grid_dims, |
---|
388 | 'src_grid_center_lon' : src_grid_center_lon, |
---|
389 | 'src_grid_center_lat' : src_grid_center_lat, |
---|
390 | 'src_grid_corner_lon' : src_grid_corner_lon, |
---|
391 | 'src_grid_corner_lat' : src_grid_corner_lat, |
---|
392 | 'src_grid_area' : src_grid_area, |
---|
393 | 'src_grid_area' : src_grid_area, |
---|
394 | 'src_grid_pmask' : src_grid_pmask, |
---|
395 | 'dst_grid_dims' : dst_grid_dims, |
---|
396 | 'dst_grid_center_lon' : dst_grid_center_lon, |
---|
397 | 'st_grid_center_lat' : dst_grid_center_lat, |
---|
398 | 'dst_grid_corner_lon' : dst_grid_corner_lon, |
---|
399 | 'dst_grid_corner_lat' : dst_grid_corner_lat, |
---|
400 | 'dst_grid_area' : dst_grid_area, |
---|
401 | 'dst_grid_imask' : dst_grid_imask, |
---|
402 | 'dst_grid_pmask' : dst_grid_pmask, |
---|
403 | 'src_lon_addressed' : src_lon_addressed, |
---|
404 | 'src_lat_addressed' : src_lat_addressed, |
---|
405 | 'src_area_addressed' : src_area_addressed, |
---|
406 | 'dst_lon_addressed' : dst_lon_addressed, |
---|
407 | 'dst_lat_addressed' : dst_lat_addressed, |
---|
408 | 'dst_area_addressed' : dst_area_addressed, |
---|
409 | 'dst_imask_addressed' : dst_imask_addressed, |
---|
410 | 'dst_pmask_addressed' : dst_pmask_addressed, |
---|
411 | 'atmCoast' : atmCoast, |
---|
412 | 'oceLand' : oceLand, |
---|
413 | 'oceOcean' : oceOcean, |
---|
414 | 'oceOceanFiltered' : oceOceanFiltered, |
---|
415 | 'oceCoast' : oceCoast |
---|
416 | } ) |
---|
417 | |
---|
418 | f_runoff.attrs['Conventions'] = "CF-1.6" |
---|
419 | f_runoff.attrs['source'] = "IPSL Earth system model" |
---|
420 | f_runoff.attrs['group'] = "ICMC IPSL Climate Modelling Center" |
---|
421 | f_runoff.attrs['Institution'] = "IPSL https.//www.ipsl.fr" |
---|
422 | f_runoff.attrs['Ocean'] = oce_Name + " https://www.nemo-ocean.eu" |
---|
423 | f_runoff.attrs['Atmosphere'] = atm_Name + " http://lmdz.lmd.jussieu.fr" |
---|
424 | f_runoff.attrs['associatedFiles'] = grids + " " + areas + " " + masks |
---|
425 | f_runoff.attrs['description'] = "Generated with RunoffWeights.py" |
---|
426 | f_runoff.attrs['title'] = runoff |
---|
427 | f_runoff.attrs['Program'] = "Generated by " + sys.argv[0] + " with flags " + ' '.join (sys.argv[1:]) |
---|
428 | f_runoff.attrs['atmCoastWidth'] = "{:d} grid points".format(atmCoastWidth) |
---|
429 | f_runoff.attrs['oceCoastWidth'] = "{:d} grid points".format(oceCoastWidth) |
---|
430 | f_runoff.attrs['searchRadius'] = "{:.0f} km".format(searchRadius/1000.) |
---|
431 | f_runoff.attrs['atmQuantity'] = myargs.atmQuantity |
---|
432 | f_runoff.attrs['oceQuantity'] = myargs.oceQuantity |
---|
433 | f_runoff.attrs['gridsFile'] = grids |
---|
434 | f_runoff.attrs['areasFile'] = areas |
---|
435 | f_runoff.attrs['masksFile'] = masks |
---|
436 | f_runoff.attrs['o2aFile'] = o2a |
---|
437 | f_runoff.attrs['timeStamp'] = time.asctime () |
---|
438 | try : f_calving.attrs['directory'] = os.getcwd () |
---|
439 | except : pass |
---|
440 | try : f_runoff.attrs['HOSTNAME'] = platform.node () |
---|
441 | except : pass |
---|
442 | try : f_runoff.attrs['LOGNAME'] = os.getlogin () |
---|
443 | except : pass |
---|
444 | try : f_runoff.attrs['Python'] = "Python version " + platform.python_version () |
---|
445 | except : pass |
---|
446 | try : f_runoff.attrs['OS'] = platform.system () |
---|
447 | except : pass |
---|
448 | try : f_runoff.attrs['release'] = platform.release () |
---|
449 | except : pass |
---|
450 | try : f_runoff.attrs['hardware'] = platform.machine () |
---|
451 | except : pass |
---|
452 | f_runoff.attrs['conventions'] = "SCRIP" |
---|
453 | f_runoff.attrs['source_grid'] = "curvilinear" |
---|
454 | f_runoff.attrs['dest_grid'] = "curvilinear" |
---|
455 | f_runoff.attrs['Model'] = "IPSL CM6" |
---|
456 | f_runoff.attrs['SVN_Author'] = "$Author$" |
---|
457 | f_runoff.attrs['SVN_Date'] = "$Date$" |
---|
458 | f_runoff.attrs['SVN_Revision'] = "$Revision$" |
---|
459 | f_runoff.attrs['SVN_Id'] = "$Id$" |
---|
460 | f_runoff.attrs['SVN_HeadURL'] = "$HeadURL$" |
---|
461 | |
---|
462 | f_runoff.to_netcdf ( runoff, mode='w', format=FmtNetcdf ) |
---|
463 | f_runoff.close () |
---|
464 | |
---|
465 | ## |
---|
466 | |
---|
467 | print ('That''s all folks !') |
---|
468 | ## ====================================================================================== |
---|