source: TOOLS/MOSAIX/CalvingWeights.py @ 4199

Last change on this file since 4199 was 4199, checked in by omamce, 3 years ago

O.M. : MOASAIX, suppress generation of NetCDF3 files

  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 22.1 KB
Line 
1### ===========================================================================
2###
3### Compute calving weights.
4###
5### ===========================================================================
6##
7##  Warning, to install, configure, run, use any of Olivier Marti's
8##  software or to read the associated documentation you'll need at least
9##  one (1) brain in a reasonably working order. Lack of this implement
10##  will void any warranties (either express or implied).
11##  O. Marti assumes no responsability for errors, omissions,
12##  data loss, or any other consequences caused directly or indirectly by
13##  the usage of his software by incorrectly or partially configured
14##  personal.
15##
16import netCDF4, numpy as np
17import sys, os, platform, argparse, textwrap, time
18from scipy import ndimage
19import nemo
20
21## SVN information
22__Author__   = "$Author$"
23__Date__     = "$Date$"
24__Revision__ = "$Revision$"
25__Id__       = "$Id$"
26__HeadURL__  = "$HeadURL$"
27__SVN_Date__ = "$SVN_Date: $"
28
29###
30
31### ===== Handling command line parameters ==================================================
32# Creating a parser
33parser = argparse.ArgumentParser (
34    description = """Compute calving weights""",
35    epilog='-------- End of the help message --------')
36
37# Adding arguments
38parser.add_argument ('--oce'     , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025'] )
39parser.add_argument ('--atm'     , help='atm model name (ICO* or LMD*)', type=str, default='ICO40'    )
40parser.add_argument ('--type'    , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full'  )
41parser.add_argument ('--dir'     , help='Directory of input file', type=str, default='./' )
42parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' )
43parser.add_argument ('--repartition_var' , help='Variable name for iceshelf'      , type=str, default=None)
44parser.add_argument ('--output'  , help='output rmp file name', default='rmp_tlmd_to_torc_calving.nc' )
45parser.add_argument ('--fmt'     , help='NetCDF file format, using nco syntax', default='netcdf4', choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] )
46
47# Parse command line
48myargs = parser.parse_args()
49
50# Model Names
51src_Name = myargs.atm ; dst_Name = myargs.oce
52   
53# Default vars
54if myargs.repartition_var == None :
55    # Runoff from Antarctica iceshelves (Depoorter, 2013)
56    if myargs.type == 'iceshelf' : repartitionVar = 'sornfisf'
57       
58    # Runoff from Antarctica iceberg (Depoorter, 2013)
59    if myargs.type == 'iceberg' :  repartitionVar = 'Icb_Flux'
60       
61else :                             repartitionVar = myargs.repartition_var
62
63# Latitude limits of each calving zone
64limit_lat = ( (40.0, 90.0), (-50.0, 40.0), ( -90.0, -50.0) )
65
66nb_zone = len(limit_lat)
67
68if myargs.fmt == 'classic'         : FmtNetcdf = 'CLASSIC'
69if myargs.fmt == 'netcdf3'         : FmtNetcdf = 'CLASSIC'
70if myargs.fmt == '64bit'           : FmtNetcdf = 'NETCDF3_64BIT_OFFSET'
71if myargs.fmt == '64bit_data'      : FmtNetcdf = 'NETCDF3_64BIT_DATA'
72if myargs.fmt == '64bit_offset'    : FmtNetcdf = 'NETCDF3_64BIT_OFFSET'
73if myargs.fmt == 'netcdf4'         : FmtNetcdf = 'NETCDF4'
74if myargs.fmt == 'netcdf4_classic' : FmtNetcdf = 'NETCDF4_CLASSIC'
75### ==========================================================================================
76
77# Model short names
78src_name = None ; dst_name = None
79if src_Name.count('ICO')  != 0 : src_name = 'ico'
80if src_Name.count('LMD')  != 0 : src_name = 'lmd'
81if dst_Name.count('ORCA') != 0 : dst_name = 'orc'
82
83CplModel = dst_Name + 'x' + src_Name
84
85print ('Atm names : ' + src_name + ' ' + src_Name )
86print ('Oce names : ' + dst_name + ' ' + dst_Name )
87print (' ')
88
89# Reading domains characteristics
90areas = myargs.dir + '/areas_' + dst_Name + 'x' + src_Name + '.nc' 
91masks = myargs.dir + '/masks_' + dst_Name + 'x' + src_Name + '.nc'
92grids = myargs.dir + '/grids_' + dst_Name + 'x' + src_Name + '.nc'
93
94print ('Opening ' + areas)
95f_areas = netCDF4.Dataset ( areas, "r" )
96print ('Opening ' + masks)
97f_masks = netCDF4.Dataset ( masks, "r" )
98print ('Opening ' + grids)
99f_grids = netCDF4.Dataset ( grids, "r" )
100
101src_lon = f_grids.variables ['t' + src_name + '.lon'][:]
102src_lat = f_grids.variables ['t' + src_name + '.lat'][:]
103dst_lon = f_grids.variables ['t' + dst_name + '.lon'][:]
104dst_lat = f_grids.variables ['t' + dst_name + '.lat'][:]
105
106src_msk = f_masks.variables ['t' + src_name + '.msk'][:]
107dst_msk = f_masks.variables ['t' + dst_name + '.msk'][:]
108dst_mskutil = 1-dst_msk # Reversing the convention : 0 on continent, 1 on ocean
109print ('dst_msk     : ' + str(np.sum(dst_msk)))
110print ('dst_mskutil : ' + str(np.sum(dst_mskutil)))
111
112# Periodicity masking for NEMO
113if dst_Name == 'ORCA2.3'  : nperio_dst = 4
114if dst_Name == 'eORCA1.2' : nperio_dst = 6
115if dst_Name == 'ORCA025'  : nperio_dst = 6
116if dst_Name == 'eORCA025' : nperio_dst = 6
117print ('nperio_dst: ' + str(nperio_dst) )
118dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T' )
119print ('dst_mskutil : ' + str(np.sum(dst_mskutil)))
120
121##
122## Fill Closed seas and other
123## ==================================================
124
125# Preparation by closing some straits
126# -----------------------------------
127if dst_Name == 'eORCA025' :
128    print ('Filling some seas for eORCA025')
129    # Set Gibraltar strait to 0 to fill Mediterranean sea
130    dst_mskutil[838:841,1125] = 0
131    # Set Bal-El-Mandeb to 0 to fill Red Sea
132    dst_mskutil[736,1321:1324] = 0
133    # Set Stagerak to 0 to fill Baltic Sea
134    dst_mskutil[961,1179:1182] = 0
135    # Set Ormuz Strait to 0 to fill Arabo-Persian Gulf
136    dst_mskutil[794:798,1374] = 0
137    # Set Hudson Strait to 0 to fill Hudson Bay
138    dst_mskutil[997:1012,907] = 0
139   
140if dst_Name == 'eORCA1.2' :
141    print ('Filling some seas for eORCA1.2')
142    # Set Gibraltar strait to 0 to fill Mediterrean sea
143    dst_mskutil[240, 283] = 0
144    # Set Bal-El-Mandeb to 0 to fill Red Sea
145    dst_mskutil[211:214, 332] = 0
146    # Set Stagerak to 0 to fill Baltic Sea
147    dst_mskutil[272:276, 293] = 0
148    # Set Ormuz Strait to 0 to fill Arabo-Persian Gulf
149    dst_mskutil[227:230, 345] = 0
150    # Set Hudson Strait to 0 to fill Hudson Bay
151    dst_mskutil[284,222:225] = 0
152   
153if dst_Name == 'ORCA2.3' :
154    print ('Filling some seas for ORCA2.3')
155    # Set Gibraltar strait to 0 to fill Mediterrean sea
156    dst_mskutil[101,139] = 0
157    # Set Black Sea to zero. At the edge of the domain : binary_fill_holes fails
158    dst_mskutil[ 99:111,   0:  5] = 0
159    dst_mskutil[106:111, 173:182] = 0
160    # Set Stagerak to 0 to fill Baltic Sea
161    dst_mskutil[115,144] = 0
162    # Set Hudson Strait to 0 to fill Hudson Bay
163    dst_mskutil[120:123,110] = 0
164    # Set Bal-El-Mandeb to 0 to fill Red Sea
165    dst_mskutil[87:89,166] = 0
166
167dst_closed = dst_mskutil
168
169# Fill closed seas with image processing library
170# ----------------------------------------------
171dst_mskutil = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(dst_mskutil, nperio=nperio_dst, cd_type='T')), nperio=nperio_dst, cd_type='T' )
172
173# Surfaces
174src_srf = f_areas.variables ['t' + src_name + '.srf']
175dst_srf = f_areas.variables ['t' + dst_name + '.srf']
176dst_srfutil = dst_srf * np.float64 (dst_mskutil)
177
178dst_srfutil_sum = np.sum( dst_srfutil)
179
180src_clo = f_grids.variables ['t' + src_name + '.clo'][:]
181src_cla = f_grids.variables ['t' + src_name + '.cla'][:]
182dst_clo = f_grids.variables ['t' + dst_name + '.clo'][:]
183dst_cla = f_grids.variables ['t' + dst_name + '.cla'][:]
184
185# Indices
186( src_jpj, src_jpi) = src_lat.shape ; src_grid_size = src_jpj*src_jpi
187( dst_jpj, dst_jpi) = dst_lat.shape ; dst_grid_size = dst_jpj*dst_jpi
188orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 # Fortran indices (starting at one)
189
190### ===== Reading needed data ==================================================
191if myargs.type in ['iceberg', 'iceshelf' ]: 
192    # Reading data file for calving or iceberg geometry around Antarctica
193    print ( 'Opening ' + myargs.repartition_file)
194    f_repartition = netCDF4.Dataset ( myargs.repartition_file, "r" )
195    repartition = np.sum ( f_repartition.variables [repartitionVar][:], axis=0 )
196
197## Before loop on basins
198remap_matrix = np.empty ( shape=(0), dtype=np.float64 )
199src_address  = np.empty ( shape=(0), dtype=np.int32   )
200dst_address  = np.empty ( shape=(0), dtype=np.int32   )
201
202print (' ')
203
204### ===== Starting loop on basins==============================================
205
206# Initialise some fields
207remap_matrix = np.empty ( shape=(0), dtype=np.float64 )
208src_address  = np.empty ( shape=(0), dtype=np.int32   )
209dst_address  = np.empty ( shape=(0), dtype=np.int32   )
210
211basin_msk       = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float32)
212key_repartition = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float64)
213
214## Loop
215for n_bas in np.arange ( nb_zone ) :
216    south = False ; ok_run = False
217    lat_south = np.min(limit_lat[n_bas]) ; lat_north = np.max(limit_lat[n_bas])
218    if lat_south <= -60.0 : south = True
219   
220    print ( 'basin: {:2d} -- Latitudes: {:+.1f} {:+.1f} --'.format(n_bas, lat_south, lat_north) )
221    ##
222    if myargs.type == 'iceberg'  and     south : ok_run = True  ; print ('Applying iceberg to south' )
223    if myargs.type == 'iceshelf' and     south : ok_run = True  ; print ('Applying iceshelf to south' )
224    if myargs.type == 'iceberg'  and not south : ok_run = False ; print ('Skipping iceberg: not south ')
225    if myargs.type == 'iceshelf' and not south : ok_run = False ; print ('Skipping iceshelf: not south ')
226    if myargs.type == 'nosouth'  and     south : ok_run = False ; print ('Skipping south: nosouth case' )
227    if myargs.type == 'nosouth'  and not south : ok_run = True  ; print ('Running: not in south, uniform repartition')
228    if myargs.type == 'full'                   : ok_run = True  ; print ('Running general trivial case, uniform repartition' )
229
230    if ok_run :
231        # Calving integral send to one point per latitude bands
232        index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1
233   
234        # Basin mask
235        basin_msk[n_bas] = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 )
236
237        # Repartition pattern
238        if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition[n_bas] = repartition * basin_msk[n_bas]
239        else :                                       key_repartition[n_bas] =               basin_msk[n_bas]
240
241        # Integral and normalisation
242        sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil )
243        key_repartition = key_repartition / sum_repartition
244   
245        print ( 'Sum of repartition key                      : {:12.3e}'.format (np.sum (key_repartition[n_bas]               )) )
246        print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] * dst_srfutil )) )
247   
248        # Basin surface (modulated by repartition key)
249        basin_srfutil = np.sum ( key_repartition[n_bas] * dst_srfutil )
250           
251        # Weights and links
252        poids = 1.0 / basin_srfutil
253        matrix_local   = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids, 0. )
254        matrix_local   = matrix_local[matrix_local > 0.0] # Keep only non zero values
255        num_links      = np.count_nonzero (matrix_local)
256        # Address on source grid : all links points to the same LMDZ point.
257        src_address_local = np.ones(num_links, dtype=np.int32 )*index_src
258        # Address on destination grid : each NEMO point with non zero link
259        dst_address_local = np.where ( key_repartition[n_bas].ravel() > 0.0, orc_index, 0)
260        dst_address_local = dst_address_local[dst_address_local > 0]
261        # Append to global tabs
262        remap_matrix = np.append ( remap_matrix, matrix_local      )
263        src_address  = np.append ( src_address , src_address_local )
264        dst_address  = np.append ( dst_address , dst_address_local )
265        #
266        #print ( 'Sum of remap_matrix : {:12.3e}'.format(np.sum(matrix_local)) )
267        print ( 'Point in atmospheric grid : {:4d} -- num_links: {:6d}'.format(index_src, num_links) ) 
268        print ('  ')
269
270
271
272## End of loop
273print (' ')
274
275#
276num_links = np.count_nonzero (remap_matrix)
277print ( 'Total num_links : {:10d}'.format(num_links) )
278
279# Diag : interpolate uniform field
280src_field = np.zeros(shape=(src_grid_size))
281for n_bas in np.arange ( nb_zone ) :
282    index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1
283    src_field[index_src-1] = n_bas
284
285dst_field = np.zeros(shape=(dst_grid_size,))
286for link in np.arange (num_links) :
287    dst_field [dst_address [link]-1] +=  remap_matrix[link] * src_field[src_address[link]-1]
288
289### ===== Writing the weights file, for OASIS MCT ==========================================
290
291# Output file
292calving = myargs.output
293f_calving = netCDF4.Dataset ( calving, 'w', format=FmtNetcdf )
294   
295print ('Output file: ' + calving )
296
297f_calving.Conventions     = "CF-1.6"
298f_calving.source          = "IPSL Earth system model"
299f_calving.group           = "ICMC IPSL Climate Modelling Center"
300f_calving.Institution     = "IPSL https.//www.ipsl.fr"
301f_calving.Ocean           = dst_Name + " https://www.nemo-ocean.eu"
302f_calving.Atmosphere      = src_Name + " http://lmdz.lmd.jussieu.fr"
303if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.repartition_file
304f_calving.associatedFiles = grids + " " + areas + " " + masks
305f_calving.directory       = os.getcwd ()
306f_calving.description     = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX"
307f_calving.title           = calving
308f_calving.Program         = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:])
309f_calving.repartitionType = myargs.type
310if myargs.type in [ 'iceberg', 'iceshelf' ] :
311    f_calving.repartitionFile = myargs.repartition_file
312    f_calving.repartitionVar  = repartitionVar
313f_calving.gridsFile       = grids
314f_calving.areasFile       = areas
315f_calving.masksFile       = masks
316f_calving.timeStamp       = time.asctime()
317f_calving.uuid            = f_areas.uuid
318f_calving.HOSTNAME        = platform.node()
319#f_calving.LOGNAME         = os.getlogin()
320f_calving.Python          = "Python version " +  platform.python_version()
321f_calving.OS              = platform.system()
322f_calving.release         = platform.release()
323f_calving.hardware        = platform.machine()
324f_calving.conventions     = "SCRIP"
325if src_name == 'lmd' : f_calving.source_grid = "curvilinear"
326if src_name == 'ico' : f_calving.source_grid = "unstructured"
327f_calving.dest_grid       = "curvilinear"
328f_calving.Model           = "IPSL CM6"
329f_calving.SVN_Author      = "$Author$"
330f_calving.SVN_Date        = "$Date$"
331f_calving.SVN_Revision    = "$Revision$"
332f_calving.SVN_Id          = "$Id$"
333f_calving.SVN_HeadURL     = "$HeadURL$"
334
335d_nb_zone   = f_calving.createDimension ('nb_zone'         ,   nb_zone )
336d_num_links = f_calving.createDimension ('num_links'       , num_links )
337d_num_wgts  = f_calving.createDimension ('num_wgts'        ,         1 )
338
339d_src_grid_size    = f_calving.createDimension ('src_grid_size'   , src_grid_size )
340d_src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0]  )
341d_src_grid_rank    = f_calving.createDimension ('src_grid_rank'   ,        2  )
342
343d_dst_grid_size    = f_calving.createDimension ('dst_grid_size'   , dst_grid_size )
344d_dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] )
345d_dst_grid_rank    = f_calving.createDimension ('dst_grid_rank'   ,        2  )
346
347v_remap_matrix = f_calving.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') )
348v_src_address  = f_calving.createVariable ( 'src_address' , 'i4', ('num_links',) )
349v_src_address  = f_calving.createVariable ( 'dst_address' , 'i4', ('num_links',) )
350
351v_remap_matrix[:] = remap_matrix
352v_src_address [:] = src_address
353v_src_address [:] = src_address
354
355v_src_grid_dims       = f_calving.createVariable ( 'src_grid_dims'      , 'i4', ('src_grid_rank',) )
356v_src_grid_center_lon = f_calving.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) )
357v_src_grid_center_lat = f_calving.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) )
358v_src_grid_center_lon.units='degrees_east'  ; v_src_grid_center_lon.long_name='Longitude'
359v_src_grid_center_lon.long_name='longitude' ; v_src_grid_center_lon.bounds="src_grid_corner_lon" 
360v_src_grid_center_lat.units='degrees_north' ; v_src_grid_center_lat.long_name='Latitude'
361v_src_grid_center_lat.long_name='latitude ' ; v_src_grid_center_lat.bounds="src_grid_corner_lat" 
362v_src_grid_corner_lon = f_calving.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners')  )
363v_src_grid_corner_lat = f_calving.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners')  )
364v_src_grid_corner_lon.units="degrees_east"
365v_src_grid_corner_lat.units="degrees_north"
366v_src_grid_area       = f_calving.createVariable ( 'src_grid_area'      , 'f8', ('src_grid_size',)  )
367v_src_grid_area.long_name="Grid area" ; v_src_grid_area.standard_name="cell_area" ; v_src_grid_area.units="m2"
368v_src_grid_imask      = f_calving.createVariable ( 'src_grid_imask'     , 'f8', ('src_grid_size',)  )
369v_src_grid_imask.long_name="Land-sea mask" ; v_src_grid_imask.units="Land:1, Ocean:0"
370
371v_src_grid_dims      [:] = ( src_jpi, src_jpi ) 
372v_src_grid_center_lon[:] = src_lon[:].ravel()
373v_src_grid_center_lat[:] = src_lat[:].ravel()
374v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) )
375v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) )
376v_src_grid_area      [:] = src_srf[:].ravel()
377v_src_grid_imask     [:] = src_msk[:].ravel()
378
379# --
380
381v_dst_grid_dims       = f_calving.createVariable ( 'dst_grid_dims'      , 'i4', ('dst_grid_rank',) )
382v_dst_grid_center_lon = f_calving.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) )
383v_dst_grid_center_lat = f_calving.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) )
384v_dst_grid_center_lon.units='degrees_east'  ; v_dst_grid_center_lon.long_name='Longitude'
385v_dst_grid_center_lon.long_name='longitude' ; v_dst_grid_center_lon.bounds="dst_grid_corner_lon" 
386v_dst_grid_center_lat.units='degrees_north' ; v_dst_grid_center_lat.long_name='Latitude'
387v_dst_grid_center_lat.long_name='latitude'  ; v_dst_grid_center_lat.bounds="dst_grid_corner_lat" 
388v_dst_grid_corner_lon = f_calving.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
389v_dst_grid_corner_lat = f_calving.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
390v_dst_grid_corner_lon.units="degrees_east"
391v_dst_grid_corner_lat.units="degrees_north"
392v_dst_grid_area       = f_calving.createVariable ( 'dst_grid_area'      , 'f8', ('dst_grid_size',)  )
393v_dst_grid_area.long_name="Grid area" ; v_dst_grid_area.standard_name="cell_area" ; v_dst_grid_area.units="m2"
394v_dst_grid_imask      = f_calving.createVariable ( 'dst_grid_imask'     , 'f8', ('dst_grid_size',)  )
395v_dst_grid_imask.long_name="Land-sea mask" ; v_dst_grid_imask.units="Land:1, Ocean:0"
396
397v_dst_bassin      = f_calving.createVariable ( 'dst_bassin'      , 'f8', ('nb_zone', 'dst_grid_size',)  )
398v_dst_repartition = f_calving.createVariable ( 'dst_repartition' , 'f8', ('nb_zone', 'dst_grid_size',)  )
399
400v_dst_southLimit = f_calving.createVariable ('dst_southLimit', 'f4', ('nb_zone',) )
401v_dst_northLimit = f_calving.createVariable ('dst_northLimit', 'f4', ('nb_zone',) )
402v_dst_southLimit[:] = np.min(limit_lat, axis=(1,) ) 
403v_dst_northLimit[:] = np.max(limit_lat, axis=(1,) ) 
404
405v_dst_grid_dims      [:] = ( dst_jpi, dst_jpi ) 
406v_dst_grid_center_lon[:] = dst_lon[:].ravel()
407v_dst_grid_center_lat[:] = dst_lat[:].ravel()
408v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) )
409v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) )
410v_dst_grid_area      [:] = dst_srf[:].ravel()
411v_dst_grid_imask     [:] = dst_msk[:].ravel()
412
413v_dst_bassin[:]      = basin_msk.reshape      ( (nb_zone,dst_grid_size) )
414v_dst_repartition[:] = key_repartition.reshape( (nb_zone,dst_grid_size) )
415
416# Additionnal fields for debugging
417# ================================
418v_dst_grid_imaskutil      = f_calving.createVariable ( 'dst_grid_imaskutil'     , 'f8', ('dst_grid_size',)  )
419v_dst_grid_imaskutil.long_name="Land-sea mask" ; v_dst_grid_imaskutil.units="Land:1, Ocean:0"
420v_dst_grid_imaskclose      = f_calving.createVariable ( 'dst_grid_imaskclose'     , 'f8', ('dst_grid_size',)  )
421v_dst_grid_imaskclose.long_name="Land-sea mask" ; v_dst_grid_imaskclose.units="Land:1, Ocean:0"
422v_dst_grid_imaskutil [:] = dst_mskutil[:].ravel()
423v_dst_grid_imaskclose[:] = dst_closed[:].ravel()
424
425v_dst_lon_addressed = f_calving.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) )
426v_dst_lat_addressed = f_calving.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) )
427v_dst_lon_addressed.long_name="Longitude" ; v_dst_lon_addressed.standard_name="longitude" ; v_dst_lon_addressed.units="degrees_east"
428v_dst_lat_addressed.long_name="Latitude"  ; v_dst_lat_addressed.standard_name="latitude"  ; v_dst_lat_addressed.units="degrees_north"
429v_dst_area_addressed  = f_calving.createVariable ( 'dst_area_addressed', 'f8', ('num_links',) )
430v_dst_area_addressed.long_name="Grid area" ; v_dst_area_addressed.standard_name="cell_area" ; v_dst_grid_area.units="m2"
431v_dst_imask_addressed = f_calving.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) )
432v_dst_imask_addressed.long_name="Land-sea mask" ; v_dst_imask_addressed.units="Land:1, Ocean:0"
433v_dst_imaskutil_addressed = f_calving.createVariable ( 'dst_imaskutil_addressed', 'i4', ('num_links',) )
434v_dst_imaskutil_addressed.long_name="Land-sea mask" ; v_dst_imaskutil_addressed.units="Land:1, Ocean:0"
435
436v_dst_lon_addressed  [:] = dst_lon.ravel()[dst_address-1].ravel()
437v_dst_lat_addressed  [:] = dst_lat.ravel()[dst_address-1].ravel()
438v_dst_area_addressed [:] = dst_srf[:].ravel()[dst_address-1].ravel()
439v_dst_imask_addressed[:] = dst_msk[:].ravel()[dst_address-1].ravel()
440v_dst_imaskutil_addressed[:] = dst_mskutil.ravel()[dst_address-1].ravel()
441
442v_src_field = f_calving.createVariable ( 'src_field', 'f8', ('src_grid_size',) )
443v_dst_field = f_calving.createVariable ( 'dst_field', 'f8', ('dst_grid_size',) )
444v_src_field[:] = src_field
445v_dst_field[:] = dst_field
446
447f_calving.close ()
448
449
450### ===== That's all Folk's !! ==============================================
Note: See TracBrowser for help on using the repository browser.