source: TOOLS/MOSAIX/CalvingWeights.py @ 4148

Last change on this file since 4148 was 4146, checked in by omamce, 6 years ago

O.M. : small corrections on diagnostics only

  • Property svn:keywords set to Date Revision HeadURL Author Id
File size: 17.0 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='-------- This is the 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 inout file', type=str, default='./' )
42parser.add_argument ('--isf_icb', help='Data files with iceberg and iceshelf melting', type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' )
43
44# Parse command line
45myargs = parser.parse_args()
46
47# Model Names
48src_Name = myargs.atm
49dst_Name = myargs.oce
50
51# Latitude limits of each calving zone
52limit_lat = ( (40.0, 90.0), (-50.0, 40.0), ( -90.0, -50.0) )
53nb_zone = len(limit_lat)
54
55### ==========================================================================================
56
57# Model short names
58src_name = None ; dst_name = None
59if src_Name.count('ICO')  != 0 : src_name = 'ico'
60if src_Name.count('LMD')  != 0 : src_name = 'lmd'
61if dst_Name.count('ORCA') != 0 : dst_name = 'orc'
62
63CplModel = dst_Name + 'x' + src_Name
64
65print ('Atm names : ' + src_name + ' ' + src_Name )
66print ('Oce names : ' + dst_name + ' ' + dst_Name )
67print (' ')
68
69# Reading domains characteristics
70areas = myargs.dir + '/areas_' + dst_Name + 'x' + src_Name + '.nc' 
71masks = myargs.dir + '/masks_' + dst_Name + 'x' + src_Name + '.nc'
72grids = myargs.dir + '/grids_' + dst_Name + 'x' + src_Name + '.nc'
73
74print ('Opening ' + areas)
75f_areas = netCDF4.Dataset ( areas, "r" )
76print ('Opening ' + masks)
77f_masks = netCDF4.Dataset ( masks, "r" )
78print ('Opening ' + grids)
79f_grids = netCDF4.Dataset ( grids, "r" )
80
81src_lon = f_grids.variables [ 't' + src_name + '.lon' ][:]
82src_lat = f_grids.variables [ 't' + src_name + '.lat' ][:]
83dst_lon = f_grids.variables [ 't' + dst_name + '.lon' ][:]
84dst_lat = f_grids.variables [ 't' + dst_name + '.lat' ][:]
85
86src_msk = f_masks.variables [ 't' + src_name + '.msk' ][:]
87dst_msk = f_masks.variables [ 't' + dst_name + '.msk' ][:]
88dst_mskutil = 1-dst_msk # Reversing the convention : 0 on continent, 1 on ocean
89
90# Periodicityt masking for NEMO
91if dst_Name == 'ORCA2.3'  : nperio_dst = 4
92if dst_Name == 'eORCA1.2' : nperio_dst = 6
93if dst_Name == 'eORCA025' : nperio_dst = 6
94dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T' )
95
96# Fill Closed seas : preparation by closing some straits
97if dst_Name == 'eORCA1.2' :
98    # Set Gibraltar strait to 0 to fill Mediterrean sea
99    dst_mskutil[240, 283] = 0
100    # Set Bal-El-Mandeb to 0 to fill Red Sea
101    dst_mskutil[211:214, 332] = 0
102    # Set Stagerak to 0 to fill Baltic Sea
103    dst_mskutil[272:276, 293] = 0
104    # Set Ormuz Strait to 0 to fill Arabo-Persian Gulf
105    dst_mskutil[227:230, 345] = 0
106    # Set Hudson Strait to 0 to fill Hudson Bay
107    dst_mskutil[284,222:225] = 0
108   
109if dst_Name == 'ORCA2.3' :
110    # Set Gibraltar strait to 0 to fill Mediterrean sea
111    dst_mskutil[101,139] = 0
112    # Set Black Sea to zero. At the edge of the domain : binary_fill_holes fails
113    dst_mskutil[ 99:111,   0:  5] = 0
114    dst_mskutil[106:111, 173:182] = 0
115    # Set Stagerak to 0 to fill Baltic Sea
116    dst_mskutil[115,144] = 0
117    # Set Hudson Strait to 0 to fill Hudson Bay
118    dst_mskutil[120:123,110] = 0
119    # Set Bal-El-Mandeb to 0 to fill Red Sea
120    dst_mskutil[87:89,166] = 0
121
122# Fill closed sea with image processing library
123dst_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' )
124
125# Surfaces
126src_srf = f_areas.variables [ 't' + src_name + '.srf' ]
127dst_srf = f_areas.variables [ 't' + dst_name + '.srf' ]
128dst_srfutil = dst_srf * np.float64 (dst_mskutil)
129
130dst_srfutil_sum = np.sum( dst_srfutil)
131
132src_clo = f_grids.variables [ 't' + src_name + '.clo' ][:]
133src_cla = f_grids.variables [ 't' + src_name + '.cla' ][:]
134dst_clo = f_grids.variables [ 't' + dst_name + '.clo' ][:]
135dst_cla = f_grids.variables [ 't' + dst_name + '.cla' ][:]
136
137# Indices
138( src_jpj, src_jpi) = src_lat.shape ; src_grid_size = src_jpj*src_jpi
139( dst_jpj, dst_jpi) = dst_lat.shape ; dst_grid_size = dst_jpj*dst_jpi
140orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 ## Fortran indices (starting at one)
141
142### ===== Reading needed data ==================================================
143if myargs.type in ['iceberg', 'iceshelf' ]: 
144    # Reading data file for calving geometry around Antarctica
145    print ( 'Opening ' + myargs.isf_icb )
146    f_icb = netCDF4.Dataset ( myargs.isf_icb, "r" )
147    if myargs.type == 'iceshelf' :
148        # Runoff from Antarctica iceshelves (Depoorter, 2013)
149        repartition = np.sum ( f_icb.variables ['sornfisf'][:], axis=0 )
150    if myargs.type == 'iceberg' :
151        # Freshwater flux from icebergs melting
152        repartition = np.sum ( f_icb.variables ['Icb_flux'][:], axis=0 )
153
154## Before loop on basins
155remap_matrix = np.empty ( shape=(0), dtype=np.float64 )
156src_address  = np.empty ( shape=(0), dtype=np.int32   )
157dst_address  = np.empty ( shape=(0), dtype=np.int32   )
158
159print (' ')
160
161### ===== Starting loop on basins==============================================
162## Initialise some fields
163remap_matrix = np.empty ( shape=(0), dtype=np.float64 )
164src_address  = np.empty ( shape=(0), dtype=np.int32   )
165dst_address  = np.empty ( shape=(0), dtype=np.int32   )
166
167## Loop
168for n_bas in np.arange ( nb_zone ) :
169    south = False ; ok_run = False
170    lat_south = np.min(limit_lat[n_bas]) ; lat_north = np.max(limit_lat[n_bas])
171    if lat_south <= -60.0 : south = True
172   
173    print ( 'basin: {:2d} -- Latitudes: {:+.1f} {:+.1f} --'.format(n_bas, lat_south, lat_north) )
174    ##
175    if myargs.type == 'iceberg'  and     south : ok_run = True  ; print ('Applying iceberg to south' )
176    if myargs.type == 'iceshelf' and     south : ok_run = True  ; print ('Applying iceshelf to south' )
177    if myargs.type == 'iceberg'  and not south : ok_run = False ; print ('Skipping iceberg : not south ')
178    if myargs.type == 'iceshelf' and not south : ok_run = False ; print ('Skipping iceshelf : not south ')
179    if myargs.type == 'nosouth'  and     south : ok_run = False ; print ('Skipping south : nosouth case' )
180    if myargs.type == 'nosouth'  and not south : ok_run = True  ; print ('Running : not in south, uniform repartition')
181    if myargs.type == 'full'                   : ok_run = True  ; print ('Running general trivial case, uniform repartition' )
182
183    if ok_run :
184        index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1
185   
186        # Basin mask
187        basin_msk = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 )
188
189        # Repartition pattern
190        if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition = repartition * basin_msk
191        else :                                       key_repartition =               basin_msk
192
193        # Integral and normalisation
194        sum_repartition = np.sum ( key_repartition * dst_srfutil )
195        key_repartition = key_repartition / sum_repartition
196   
197        print ( 'Sum of repartition key                      : {:12.3e}'.format (np.sum (key_repartition               )) )
198        print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition * dst_srfutil )) )
199   
200        # Basin surface (modulated by repartition key)
201        basin_srfutil = np.sum ( key_repartition * dst_srfutil )
202           
203        # Weights and links
204        poids = 1.0 / basin_srfutil
205        matrix_local   = np.where ( basin_msk.ravel() > 0.5, key_repartition.ravel()*poids, 0. )
206        matrix_local   = matrix_local[matrix_local > 0.0] # Keep only non zero values
207        num_links      = np.count_nonzero (matrix_local)
208        # address on source grid : all links points to the same LMDZ point.
209        src_address_local = np.ones(num_links, dtype=np.int32 )*index_src
210        # address on destination grid : each NEMO point with non zero link
211        dst_address_local = np.where ( key_repartition.ravel() > 0.0, orc_index, 0)
212        dst_address_local = dst_address_local[dst_address_local > 0]
213        # Append to global tabs
214        remap_matrix = np.append ( remap_matrix, matrix_local      )
215        src_address  = np.append ( src_address , src_address_local )
216        dst_address  = np.append ( dst_address , dst_address_local )
217        #
218        #print ( 'Sum of remap_matrix : {:12.3e}'.format(np.sum(matrix_local)) )
219        print ( 'Point in atmospheric grid : {:4d} -- num_links: {:6d}'.format(index_src, num_links) ) 
220        print ('  ')
221
222## End of loop
223print (' ')
224
225#
226num_links = np.count_nonzero (remap_matrix)
227print ( 'Total num_links : {:10d}'.format(num_links) )
228
229### ===== Writing the weights file, for OASIS MCT ==========================================
230
231# Output file
232calving = 'rmp_t' + src_Name + '_to_' + 't' + dst_Name + '_calving_' + myargs.type + '.nc'
233f_calving = netCDF4.Dataset ( calving, 'w', format='NETCDF3_64BIT' )
234
235print ('Output file: ' + calving )
236
237f_calving.Conventions     = "CF-1.6"
238f_calving.source          = "IPSL Earth system model"
239f_calving.group           = "ICMC IPSL Climate Modelling Center"
240f_calving.Institution     = "IPSL https.//www.ipsl.fr"
241f_calving.Ocean           = dst_Name + " https://www.nemo-ocean.eu"
242f_calving.Atmosphere      = src_Name + " http://lmdz.lmd.jussieu.fr"
243if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.isf_icb
244f_calving.associatedFiles = grids + " " + areas + " " + masks
245f_calving.directory       = os.getcwd ()
246f_calving.description     = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX"
247f_calving.title           = calving
248f_calving.Program         = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:])
249f_calving.timeStamp       = time.asctime()
250f_calving.uuid            = f_areas.uuid
251f_calving.HOSTNAME        = platform.node()
252#f_calving.LOGNAME         = os.getlogin()
253f_calving.Python          = "Python version " +  platform.python_version()
254f_calving.OS              = platform.system()
255f_calving.release         = platform.release()
256f_calving.hardware        = platform.machine()
257f_calving.Comment         = "Preliminary attempt - Do not trust !"
258f_calving.conventions     = "SCRIP"
259if src_name == 'lmd' : f_calving.source_grid = "curvilinear"
260if src_name == 'ico' : f_calving.source_grid = "unstructured"
261f_calving.dest_grid       = "curvilinear"
262f_calving.Model           = "IPSL CM6"
263f_calving.SVN_Author      = "$Author$"
264f_calving.SVN_Date        = "$Date$"
265f_calving.SVN_Revision    = "$Revision$"
266f_calving.SVN_Id          = "$Id$"
267f_calving.SVN_HeadURL     = "$HeadURL$"
268
269num_links        = f_calving.createDimension ('num_links'       , num_links )
270num_wgts         = f_calving.createDimension ('num_wgts'        ,         1 )
271
272src_grid_size    = f_calving.createDimension ('src_grid_size'   , src_grid_size )
273src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0]  )
274src_grid_rank    = f_calving.createDimension ('src_grid_rank'   ,        2  )
275
276dst_grid_size    = f_calving.createDimension ('dst_grid_size'   , dst_grid_size )
277dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] )
278dst_grid_rank    = f_calving.createDimension ('dst_grid_rank'   ,        2  )
279
280v_remap_matrix = f_calving.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') )
281v_src_address  = f_calving.createVariable ( 'src_address' , 'i4', ('num_links',) )
282v_src_address  = f_calving.createVariable ( 'dst_address' , 'i4', ('num_links',) )
283
284v_remap_matrix[:] = remap_matrix
285v_src_address [:] = src_address
286v_src_address [:] = src_address
287
288v_src_grid_dims       = f_calving.createVariable ( 'src_grid_dims'      , 'i4', ('src_grid_rank',) )
289v_src_grid_center_lon = f_calving.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) )
290v_src_grid_center_lat = f_calving.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) )
291v_src_grid_center_lon.units='degrees_east'  ; v_src_grid_center_lon.long_name='Longitude' ; v_src_grid_center_lon.long_name='longitude' ; v_src_grid_center_lon.bounds="src_grid_corner_lon" 
292v_src_grid_center_lat.units='degrees_north' ; v_src_grid_center_lat.long_name='Latitude'  ; v_src_grid_center_lat.long_name='latitude ' ; v_src_grid_center_lat.bounds="src_grid_corner_lat" 
293v_src_grid_corner_lon = f_calving.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners')  )
294v_src_grid_corner_lat = f_calving.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners')  )
295v_src_grid_corner_lon.units="degrees_east"
296v_src_grid_corner_lat.units="degrees_north"
297v_src_grid_area       = f_calving.createVariable ( 'src_grid_area'      , 'f8', ('src_grid_size',)  )
298v_src_grid_area.long_name="Grid area" ; v_src_grid_area.standard_name="cell_area" ; v_src_grid_area.units="m2"
299v_src_grid_imask      = f_calving.createVariable ( 'src_grid_imask'     , 'f8', ('src_grid_size',)  )
300v_src_grid_imask.long_name="Land-sea mask" ; v_src_grid_imask.units="Land:1, Ocean:0"
301
302v_src_grid_dims      [:] = ( src_jpi, src_jpi ) 
303v_src_grid_center_lon[:] = src_lon[:].ravel()
304v_src_grid_center_lat[:] = src_lat[:].ravel()
305v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, src_grid_corners.__len__()) )
306v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, src_grid_corners.__len__()) )
307v_src_grid_area      [:] = src_srf[:].ravel()
308v_src_grid_imask     [:] = src_msk[:].ravel()
309
310# --
311
312v_dst_grid_dims       = f_calving.createVariable ( 'dst_grid_dims'      , 'i4', ('dst_grid_rank',) )
313v_dst_grid_center_lon = f_calving.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) )
314v_dst_grid_center_lat = f_calving.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) )
315v_dst_grid_center_lon.units='degrees_east'  ; v_dst_grid_center_lon.long_name='Longitude' ; v_dst_grid_center_lon.long_name='longitude' ; v_dst_grid_center_lon.bounds="dst_grid_corner_lon" 
316v_dst_grid_center_lat.units='degrees_north' ; v_dst_grid_center_lat.long_name='Latitude'  ; v_dst_grid_center_lat.long_name='latitude'  ; v_dst_grid_center_lat.bounds="dst_grid_corner_lat" 
317v_dst_grid_corner_lon = f_calving.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
318v_dst_grid_corner_lat = f_calving.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
319v_dst_grid_corner_lon.units="degrees_east"
320v_dst_grid_corner_lat.units="degrees_north"
321v_dst_grid_area       = f_calving.createVariable ( 'dst_grid_area'      , 'f8', ('dst_grid_size',)  )
322v_dst_grid_area.long_name="Grid area" ; v_dst_grid_area.standard_name="cell_area" ; v_dst_grid_area.units="m2"
323v_dst_grid_imask      = f_calving.createVariable ( 'dst_grid_imask'     , 'f8', ('dst_grid_size',)  )
324v_dst_grid_imask.long_name="Land-sea mask" ; v_dst_grid_imask.units="Land:1, Ocean:0"
325
326v_dst_grid_dims      [:] = ( dst_jpi, dst_jpi ) 
327v_dst_grid_center_lon[:] = dst_lon[:].ravel()
328v_dst_grid_center_lat[:] = dst_lat[:].ravel()
329v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, dst_grid_corners.__len__()) )
330v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, dst_grid_corners.__len__()) )
331v_dst_grid_area      [:] = dst_srf[:].ravel()
332v_dst_grid_imask     [:] = dst_msk[:].ravel()
333
334# For diags
335v_dst_lon_addressed = f_calving.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) )
336v_dst_lat_addressed = f_calving.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) )
337v_dst_lon_addressed.long_name="Longitude" ; v_dst_lon_addressed.standard_name="longitude" ; v_dst_lon_addressed.units="degrees_east"
338v_dst_lat_addressed.long_name="Latitude"  ; v_dst_lat_addressed.standard_name="latitude"  ; v_dst_lat_addressed.units="degrees_north"
339v_dst_lon_addressed [:] = dst_lon.ravel()[dst_address-1].ravel()
340v_dst_lat_addressed [:] = dst_lat.ravel()[dst_address-1].ravel()
341
342f_calving.close ()
343
344
345### ===== That's all Folk's !! ==============================================
Note: See TracBrowser for help on using the repository browser.