Changeset 4159

Timestamp:

12/03/18 11:19:11 (5 years ago)

Author:

omamce

Message:

O.M. : correct eORC025 case for calving weights

Location:

TOOLS/MOSAIX

Files:

• CalvingWeights.py (modified) (13 diffs)
• CreateWeightsMask.bash (modified) (20 diffs)

TOOLS/MOSAIX/CalvingWeights.py

@@ -33 +33 @@
 parser = argparse.ArgumentParser (
     description = """Compute calving weights""",
-    epilog='-------- This is the end of the help message --------')
+    epilog='-------- End of the help message --------')
 
 # Adding arguments
-parser.add_argument ('--oce'   , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025'] )
-parser.add_argument ('--atm'   , help='atm model name (ICO* or LMD*)', type=str, default='ICO40'    )
-parser.add_argument ('--type'  , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full'  )
-parser.add_argument ('--dir'   , help='Directory of inout file', type=str, default='./' )
-parser.add_argument ('--isf_icb', help='Data files with iceberg and iceshelf melting', type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' )
+parser.add_argument ('--oce'     , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025'] )
+parser.add_argument ('--atm'     , help='atm model name (ICO* or LMD*)', type=str, default='ICO40'    )
+parser.add_argument ('--type'    , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full'  )
+parser.add_argument ('--dir'     , help='Directory of input file', type=str, default='./' )
+parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' )
+parser.add_argument ('--repartition_var' , help='Variable name for iceshelf'      , type=str, default=None)
 
 # Parse command line
@@ -46 +47 @@
 
 # Model Names
-src_Name = myargs.atm
-dst_Name = myargs.oce
+src_Name = myargs.atm ; dst_Name = myargs.oce
+    
+# Default vars
+if myargs.repartition_var == None :
+    # Runoff from Antarctica iceshelves (Depoorter, 2013)
+    if myargs.type == 'iceshelf' : repartitionVar = 'sornfisf'
+        
+    # Runoff from Antarctica iceberg (Depoorter, 2013)
+    if myargs.type == 'iceberg' :  repartitionVar = 'Icb_Flux'
+        
+else :                             repartitionVar = myargs.repartition_var
 
 # Latitude limits of each calving zone
 limit_lat = ( (40.0, 90.0), (-50.0, 40.0), ( -90.0, -50.0) )
+
 nb_zone = len(limit_lat)
-
 ### ==========================================================================================
 
@@ -79 +89 @@
 f_grids = netCDF4.Dataset ( grids, "r" )
 
-src_lon = f_grids.variables [ 't' + src_name + '.lon' ][:]
-src_lat = f_grids.variables [ 't' + src_name + '.lat' ][:]
-dst_lon = f_grids.variables [ 't' + dst_name + '.lon' ][:]
-dst_lat = f_grids.variables [ 't' + dst_name + '.lat' ][:]
-
-src_msk = f_masks.variables [ 't' + src_name + '.msk' ][:]
-dst_msk = f_masks.variables [ 't' + dst_name + '.msk' ][:]
+src_lon = f_grids.variables ['t' + src_name + '.lon'][:]
+src_lat = f_grids.variables ['t' + src_name + '.lat'][:]
+dst_lon = f_grids.variables ['t' + dst_name + '.lon'][:]
+dst_lat = f_grids.variables ['t' + dst_name + '.lat'][:]
+
+src_msk = f_masks.variables ['t' + src_name + '.msk'][:]
+dst_msk = f_masks.variables ['t' + dst_name + '.msk'][:]
 dst_mskutil = 1-dst_msk # Reversing the convention : 0 on continent, 1 on ocean
-
-# Periodicityt masking for NEMO
+print ('dst_msk     : ' + str(np.sum(dst_msk)))
+print ('dst_mskutil : ' + str(np.sum(dst_mskutil)))
+
+# Periodicity masking for NEMO
 if dst_Name == 'ORCA2.3'  : nperio_dst = 4
 if dst_Name == 'eORCA1.2' : nperio_dst = 6
+if dst_Name == 'ORCA025'  : nperio_dst = 6
 if dst_Name == 'eORCA025' : nperio_dst = 6
+print ('nperio_dst: ' + str(nperio_dst) )
 dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T' )
-
-# Fill Closed seas : preparation by closing some straits
+print ('dst_mskutil : ' + str(np.sum(dst_mskutil)))
+
+##
+## Fill Closed seas and other
+## ==================================================
+
+# Preparation by closing some straits
+# -----------------------------------
+if dst_Name == 'eORCA025' :
+    print ('Filling some seas for eORCA025')
+    # Set Gibraltar strait to 0 to fill Mediterranean sea
+    dst_mskutil[838:841,1125] = 0
+    # Set Bal-El-Mandeb to 0 to fill Red Sea
+    dst_mskutil[736,1321:1324] = 0
+    # Set Stagerak to 0 to fill Baltic Sea
+    dst_mskutil[961,1179:1182] = 0
+    # Set Ormuz Strait to 0 to fill Arabo-Persian Gulf
+    dst_mskutil[794:798,1374] = 0
+    # Set Hudson Strait to 0 to fill Hudson Bay
+    dst_mskutil[997:1012,907] = 0
+    
 if dst_Name == 'eORCA1.2' :
+    print ('Filling some seas for eORCA1.2')
     # Set Gibraltar strait to 0 to fill Mediterrean sea
     dst_mskutil[240, 283] = 0
@@ -108 +142 @@
     
 if dst_Name == 'ORCA2.3' :
+    print ('Filling some seas for ORCA2.3')
     # Set Gibraltar strait to 0 to fill Mediterrean sea
     dst_mskutil[101,139] = 0
@@ -120 +155 @@
     dst_mskutil[87:89,166] = 0
 
-# Fill closed sea with image processing library
+dst_closed = dst_mskutil
+
+# Fill closed seas with image processing library
+# ----------------------------------------------
 dst_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' )
 
 # Surfaces
-src_srf = f_areas.variables [ 't' + src_name + '.srf' ]
-dst_srf = f_areas.variables [ 't' + dst_name + '.srf' ]
+src_srf = f_areas.variables ['t' + src_name + '.srf']
+dst_srf = f_areas.variables ['t' + dst_name + '.srf']
 dst_srfutil = dst_srf * np.float64 (dst_mskutil)
 
 dst_srfutil_sum = np.sum( dst_srfutil)
 
-src_clo = f_grids.variables [ 't' + src_name + '.clo' ][:]
-src_cla = f_grids.variables [ 't' + src_name + '.cla' ][:]
-dst_clo = f_grids.variables [ 't' + dst_name + '.clo' ][:]
-dst_cla = f_grids.variables [ 't' + dst_name + '.cla' ][:]
+src_clo = f_grids.variables ['t' + src_name + '.clo'][:]
+src_cla = f_grids.variables ['t' + src_name + '.cla'][:]
+dst_clo = f_grids.variables ['t' + dst_name + '.clo'][:]
+dst_cla = f_grids.variables ['t' + dst_name + '.cla'][:]
 
 # Indices
 ( src_jpj, src_jpi) = src_lat.shape ; src_grid_size = src_jpj*src_jpi
 ( dst_jpj, dst_jpi) = dst_lat.shape ; dst_grid_size = dst_jpj*dst_jpi
-orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 ## Fortran indices (starting at one)
+orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 # Fortran indices (starting at one)
 
 ### ===== Reading needed data ==================================================
 if myargs.type in ['iceberg', 'iceshelf' ]: 
-    # Reading data file for calving geometry around Antarctica
-    print ( 'Opening ' + myargs.isf_icb )
-    f_icb = netCDF4.Dataset ( myargs.isf_icb, "r" )
-    if myargs.type == 'iceshelf' :
-        # Runoff from Antarctica iceshelves (Depoorter, 2013)
-        repartition = np.sum ( f_icb.variables ['sornfisf'][:], axis=0 )
-    if myargs.type == 'iceberg' :
-        # Freshwater flux from icebergs melting
-        repartition = np.sum ( f_icb.variables ['Icb_flux'][:], axis=0 )
+    # Reading data file for calving or iceberg geometry around Antarctica
+    print ( 'Opening ' + myargs.repartition_file)
+    f_repartition = netCDF4.Dataset ( myargs.repartition_file, "r" )
+    repartition = np.sum ( f_repartition.variables [repartitionVar][:], axis=0 )
 
 ## Before loop on basins
@@ -160 +193 @@
 
 ### ===== Starting loop on basins==============================================
-## Initialise some fields
+
+# Initialise some fields
 remap_matrix = np.empty ( shape=(0), dtype=np.float64 )
 src_address  = np.empty ( shape=(0), dtype=np.int32   )
 dst_address  = np.empty ( shape=(0), dtype=np.int32   )
+
+basin_msk       = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float32)
+key_repartition = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float64)
 
 ## Loop
@@ -175 +212 @@
     if myargs.type == 'iceberg'  and     south : ok_run = True  ; print ('Applying iceberg to south' )
     if myargs.type == 'iceshelf' and     south : ok_run = True  ; print ('Applying iceshelf to south' )
-    if myargs.type == 'iceberg'  and not south : ok_run = False ; print ('Skipping iceberg : not south ')
-    if myargs.type == 'iceshelf' and not south : ok_run = False ; print ('Skipping iceshelf : not south ')
-    if myargs.type == 'nosouth'  and     south : ok_run = False ; print ('Skipping south : nosouth case' )
-    if myargs.type == 'nosouth'  and not south : ok_run = True  ; print ('Running : not in south, uniform repartition')
+    if myargs.type == 'iceberg'  and not south : ok_run = False ; print ('Skipping iceberg: not south ')
+    if myargs.type == 'iceshelf' and not south : ok_run = False ; print ('Skipping iceshelf: not south ')
+    if myargs.type == 'nosouth'  and     south : ok_run = False ; print ('Skipping south: nosouth case' )
+    if myargs.type == 'nosouth'  and not south : ok_run = True  ; print ('Running: not in south, uniform repartition')
     if myargs.type == 'full'                   : ok_run = True  ; print ('Running general trivial case, uniform repartition' )
 
     if ok_run :
+        # Calving integral send to one point per latitude bands
         index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1
    
         # Basin mask
-        basin_msk = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 )
+        basin_msk[n_bas] = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 )
 
         # Repartition pattern
-        if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition = repartition * basin_msk
-        else :                                       key_repartition =               basin_msk
+        if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition[n_bas] = repartition * basin_msk[n_bas]
+        else :                                       key_repartition[n_bas] =               basin_msk[n_bas]
 
         # Integral and normalisation
-        sum_repartition = np.sum ( key_repartition * dst_srfutil )
+        sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil )
         key_repartition = key_repartition / sum_repartition
     
-        print ( 'Sum of repartition key                      : {:12.3e}'.format (np.sum (key_repartition               )) )
-        print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition * dst_srfutil )) )
+        print ( 'Sum of repartition key                      : {:12.3e}'.format (np.sum (key_repartition[n_bas]               )) )
+        print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] * dst_srfutil )) )
     
         # Basin surface (modulated by repartition key)
-        basin_srfutil = np.sum ( key_repartition * dst_srfutil )
+        basin_srfutil = np.sum ( key_repartition[n_bas] * dst_srfutil )
             
         # Weights and links
         poids = 1.0 / basin_srfutil
-        matrix_local   = np.where ( basin_msk.ravel() > 0.5, key_repartition.ravel()*poids, 0. )
+        matrix_local   = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids, 0. )
         matrix_local   = matrix_local[matrix_local > 0.0] # Keep only non zero values
         num_links      = np.count_nonzero (matrix_local)
-        # address on source grid : all links points to the same LMDZ point.
+        # Address on source grid : all links points to the same LMDZ point.
         src_address_local = np.ones(num_links, dtype=np.int32 )*index_src
-        # address on destination grid : each NEMO point with non zero link
-        dst_address_local = np.where ( key_repartition.ravel() > 0.0, orc_index, 0)
+        # Address on destination grid : each NEMO point with non zero link
+        dst_address_local = np.where ( key_repartition[n_bas].ravel() > 0.0, orc_index, 0)
         dst_address_local = dst_address_local[dst_address_local > 0]
         # Append to global tabs
@@ -241 +279 @@
 f_calving.Ocean           = dst_Name + " https://www.nemo-ocean.eu"
 f_calving.Atmosphere      = src_Name + " http://lmdz.lmd.jussieu.fr"
-if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.isf_icb
+if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.repartition_file
 f_calving.associatedFiles = grids + " " + areas + " " + masks
 f_calving.directory       = os.getcwd ()
@@ -267 +305 @@
 f_calving.SVN_HeadURL     = "$HeadURL$"
 
-num_links        = f_calving.createDimension ('num_links'       , num_links )
-num_wgts         = f_calving.createDimension ('num_wgts'        ,         1 )
-
-src_grid_size    = f_calving.createDimension ('src_grid_size'   , src_grid_size )
-src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0]  )
-src_grid_rank    = f_calving.createDimension ('src_grid_rank'   ,        2  )
-
-dst_grid_size    = f_calving.createDimension ('dst_grid_size'   , dst_grid_size )
-dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] )
-dst_grid_rank    = f_calving.createDimension ('dst_grid_rank'   ,        2  )
+d_nb_zone          = f_calving.createDimension ('nb_zone'         ,   nb_zone )
+d_num_links        = f_calving.createDimension ('num_links'       , num_links )
+d_num_wgts         = f_calving.createDimension ('num_wgts'        ,         1 )
+
+d_src_grid_size    = f_calving.createDimension ('src_grid_size'   , src_grid_size )
+d_src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0]  )
+D_src_grid_rank    = f_calving.createDimension ('src_grid_rank'   ,        2  )
+
+d_dst_grid_size    = f_calving.createDimension ('dst_grid_size'   , dst_grid_size )
+d_dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] )
+d_dst_grid_rank    = f_calving.createDimension ('dst_grid_rank'   ,        2  )
 
 v_remap_matrix = f_calving.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') )
@@ -289 +328 @@
 v_src_grid_center_lon = f_calving.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) )
 v_src_grid_center_lat = f_calving.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) )
-v_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" 
-v_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" 
+v_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" 
+v_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" 
 v_src_grid_corner_lon = f_calving.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners')  )
 v_src_grid_corner_lat = f_calving.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners')  )
@@ -303 +344 @@
 v_src_grid_center_lon[:] = src_lon[:].ravel()
 v_src_grid_center_lat[:] = src_lat[:].ravel()
-v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, src_grid_corners.__len__()) )
-v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, src_grid_corners.__len__()) )
+v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) )
+v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) )
 v_src_grid_area      [:] = src_srf[:].ravel()
 v_src_grid_imask     [:] = src_msk[:].ravel()
@@ -313 +354 @@
 v_dst_grid_center_lon = f_calving.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) )
 v_dst_grid_center_lat = f_calving.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) )
-v_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" 
-v_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" 
+v_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" 
+v_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" 
 v_dst_grid_corner_lon = f_calving.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
 v_dst_grid_corner_lat = f_calving.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners')  )
@@ -324 +367 @@
 v_dst_grid_imask.long_name="Land-sea mask" ; v_dst_grid_imask.units="Land:1, Ocean:0"
 
+v_dst_bassin      = f_calving.createVariable ( 'dst_bassin'      , 'f8', ('nb_zone', 'dst_grid_size',)  )
+v_dst_repartition = f_calving.createVariable ( 'dst_repartition' , 'f8', ('nb_zone', 'dst_grid_size',)  )
+
 v_dst_grid_dims      [:] = ( dst_jpi, dst_jpi ) 
 v_dst_grid_center_lon[:] = dst_lon[:].ravel()
 v_dst_grid_center_lat[:] = dst_lat[:].ravel()
-v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, dst_grid_corners.__len__()) )
-v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, dst_grid_corners.__len__()) )
+v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) )
+v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) )
 v_dst_grid_area      [:] = dst_srf[:].ravel()
 v_dst_grid_imask     [:] = dst_msk[:].ravel()
 
-# For diags
+v_dst_bassin[:]      = basin_msk.reshape      ( (nb_zone,dst_grid_size) )
+v_dst_repartition[:] = key_repartition.reshape( (nb_zone,dst_grid_size) )
+
+# Additoonnal fields for debugging
+# ================================
+v_dst_grid_imaskutil      = f_calving.createVariable ( 'dst_grid_imaskutil'     , 'f8', ('dst_grid_size',)  )
+v_dst_grid_imaskutil.long_name="Land-sea mask" ; v_dst_grid_imaskutil.units="Land:1, Ocean:0"
+v_dst_grid_imaskclose      = f_calving.createVariable ( 'dst_grid_imaskclose'     , 'f8', ('dst_grid_size',)  )
+v_dst_grid_imaskclose.long_name="Land-sea mask" ; v_dst_grid_imaskclose.units="Land:1, Ocean:0"
+v_dst_grid_imaskutil [:] = dst_mskutil[:].ravel()
+v_dst_grid_imaskclose[:] = dst_closed[:].ravel()
+
 v_dst_lon_addressed = f_calving.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) )
 v_dst_lat_addressed = f_calving.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) )
 v_dst_lon_addressed.long_name="Longitude" ; v_dst_lon_addressed.standard_name="longitude" ; v_dst_lon_addressed.units="degrees_east"
 v_dst_lat_addressed.long_name="Latitude"  ; v_dst_lat_addressed.standard_name="latitude"  ; v_dst_lat_addressed.units="degrees_north"
-v_dst_lon_addressed [:] = dst_lon.ravel()[dst_address-1].ravel()
-v_dst_lat_addressed [:] = dst_lat.ravel()[dst_address-1].ravel()
+v_dst_area_addressed  = f_calving.createVariable ( 'dst_area_addressed', 'f8', ('num_links',) )
+v_dst_area_addressed.long_name="Grid area" ; v_dst_area_addressed.standard_name="cell_area" ; v_dst_grid_area.units="m2"
+v_dst_imask_addressed = f_calving.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) )
+v_dst_imask_addressed.long_name="Land-sea mask" ; v_dst_imask_addressed.units="Land:1, Ocean:0"
+v_dst_imaskutil_addressed = f_calving.createVariable ( 'dst_imaskutil_addressed', 'i4', ('num_links',) )
+v_dst_imaskutil_addressed.long_name="Land-sea mask" ; v_dst_imaskutil_addressed.units="Land:1, Ocean:0"
+
+v_dst_lon_addressed  [:] = dst_lon.ravel()[dst_address-1].ravel()
+v_dst_lat_addressed  [:] = dst_lat.ravel()[dst_address-1].ravel()
+v_dst_area_addressed [:] = dst_srf[:].ravel()[dst_address-1].ravel()
+v_dst_imask_addressed[:] = dst_msk[:].ravel()[dst_address-1].ravel()
+v_dst_imaskutil_addressed[:] = dst_mskutil.ravel()[dst_address-1].ravel()
 
 f_calving.close ()

TOOLS/MOSAIX/CreateWeightsMask.bash

@@ -4 +4 @@
 #MSUB -e Out_WeightsMask    # Error output
 #MSUB -eo
-#MSUB -n 8                  # Number of processors
+#MSUB -n 16                 # Number of processors
 #MSUB -T 1800               # Time limit (seconds)
 #MSUB -q skylake
@@ -49 +49 @@
 ## Configuration
 ## ===========================================================================
+set -e
 
 #
@@ -54 +55 @@
 # ==============
 #OCE=ORCA2.3
-OCE=eORCA1.2
-#OCE=ORCA025
+#OCE=eORCA1.2
+OCE=eORCA025
 
 #ATM=ICO30
 #ATM=ICO40
 #ATM=ICO450
-#ATM=LMD9695
-ATM=LMD144142
+ATM=LMD9695
+#ATM=LMD144142
 
 # Default values, used to create ocean fraction on atmospheric grid 
-DefaultValues=( Direction=o2a,ocegrid=t,atmgrid=t,Order=1st,Quantity=false,Renormalize=false,useArea=false,maskSrc=true,maskDst=false )
+DefaultValues=( Direction=o2a,oceGrid=t,atmGrid=t,Order=1st,Quantity=false,Renormalize=false,useArea=false,maskSrc=true,maskDst=false )
 
 # List of weights to build
 CommandList=( \
-    Direction=o2a,Order=1st,Quantity=false,Renormalize=true,ocegrid=t,atmgrid=t,useArea=false,maskSrc=true,maskDst=true  \
-    Direction=a2o,Order=1st,Quantity=false,Renormalize=false,atmgrid=t,ocegrid=t,useArea=true,maskSrc=true,maskDst=true  \
-    Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmgrid=t,ocegrid=u,useArea=true,maskSrc=true,maskDst=true  \
-    Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmgrid=t,ocegrid=v,useArea=true,maskSrc=true,maskDst=true  \
-    Direction=a2o,Order=1st,Quantity=true,Renormalize=false,atmgrid=c,ocegrid=t,useArea=true,maskSrc=true,maskDst=true   \
-   )
+    Direction=o2a,Order=1st,Quantity=false,Renormalize=true,oceGrid=t,atmGrid=t,useArea=false,maskSrc=true,maskDst=true  \
+    Direction=a2o,Order=1st,Quantity=false,Renormalize=false,atmGrid=t,oceGrid=t,useArea=true,maskSrc=true,maskDst=true  \
+    Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmGrid=t,oceGrid=u,useArea=true,maskSrc=true,maskDst=true  \
+    Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmGrid=t,oceGrid=v,useArea=true,maskSrc=true,maskDst=true  \
+    Direction=a2o,Order=1st,Quantity=true,Renormalize=false,atmGrid=c,oceGrid=t,useArea=true,maskSrc=true,maskDst=true   \
+    )
+
+# Each item in CommandList describes the property of interpolation weights to generate. 
+# White spaces separate analysis. No spaces in any analysis.
+#
+# Keywords : 
+#   Direction   : o2a for ocean to atmosphere, a2o for atmosphere to ocean
+#   Order       : 1st or 2nd
+#   Quantity    : true if integrated quantity over a grid box, false for flux (quantity / m^2)
+#                 or intensive value (temperature, salinity, sea-ice fraction, ...)
+#   Renormalize : used when source grid is masked, to use values on non masked points only
+#   oceGrid     : t, u or v point for OPA grid
+#   atmGrid     : up to know, only t grid used in the atmosphere
+#   useArea     : if true area from the model metrics is used. If false, areas are computed by XIOS from grid corners
+#   maskSrc     : true to use the source grid mask, false to used all grid points
+#   maskDst     : true to use the destination grid mask, false to use all grid points
+#
 
 ## ===========================================================================
@@ -81 +98 @@
 ## ===========================================================================
 SUBMIT_DIR=$(pwd)
-# Function to handle command parameters
+
+# Functions to handle command parameters
+# ======================================
 function read_Command {
     # Decipher the command line to set bash variables 
@@ -92 +111 @@
     done
     local l_Command=${1}
-    for l_Element in $(echo ${l_Command} | tr "," "\n" )
-    do
+    for l_Element in $(echo ${l_Command} | tr "," "\n" ) ; do
         [[ "X${l_Debug}" = "Xtrue" ]] && echo ${l_Element}
         eval export ${l_Element}
@@ -101 +119 @@
 function setValues {
     #
-    ocegrid=${ocegrid,,}
-    atmgrid=${atmgrid,,}
-    OCEGRID=${ocegrid^^}
-    ATMGRID=${atmgrid^^}
+    oceGrid=${oceGrid,,} ; atmGrid=${atmGrid,,}
+    OCEGRID=${oceGrid^^} ; ATMGRID=${atmGrid^^}
 
     case ${Order} in
@@ -127 +143 @@
     case ${Direction} in
         ( o2a )
-        src=${oce} ; SRC=${OCE} ; srcGrid=${ocegrid} ; srcDomainType=${oceDomainType} ; SRCGRID=${OCEGRID} ; srcArea=area_grid_${OCEGRID}
-        dst=${atm} ; DST=${ATM} ; dstGrid=${atmgrid} ; dstDomainType=${atmDomainType} ; DSTGRID=${ATMGRID} ; dstArea=aire
+        src=${oce} ; SRC=${OCE} ; srcGrid=${oceGrid} ; srcDomainType=${oceDomainType} ; SRCGRID=${OCEGRID} ; srcArea=area_grid_${OCEGRID}
+        dst=${atm} ; DST=${ATM} ; dstGrid=${atmGrid} ; dstDomainType=${atmDomainType} ; DSTGRID=${ATMGRID} ; dstArea=aire
         ;; 
         ( a2o )
-        src=${atm} ; SRC=${ATM} ; srcGrid=${atmgrid} ; srcDomainType=${atmDomainType} ; SRCGRID=${ATMGRID} ; srcArea=aire
-        dst=${oce} ; DST=${OCE} ; dstGrid=${ocegrid} ; dstDomainType=${oceDomainType} ; DSTGRID=${OCEGRID} ; dstArea=area_grid_${OCEGRID}
+        src=${atm} ; SRC=${ATM} ; srcGrid=${atmGrid} ; srcDomainType=${atmDomainType} ; SRCGRID=${ATMGRID} ; srcArea=aire
+        dst=${oce} ; DST=${OCE} ; dstGrid=${oceGrid} ; dstDomainType=${oceDomainType} ; DSTGRID=${OCEGRID} ; dstArea=area_grid_${OCEGRID}
         ;;
     esac
@@ -166 +182 @@
     TMPDIR=${HOME}/Scratch/TMP
     SUBMIT_DIR=$(pwd)
-    MPIRUN=/opt/local/bin/mpirun-openmpi-gcc49 -n 2
+    MPIRUN=/opt/local/bin/mpirun -n 4
     ;; 
     ( * ) exit -1 ;;
@@ -184 +200 @@
     ( *ORC* )         oce=orc ; oceDomainType=curvilinear   ;;
 esac
+
 case ${ATM} in
     ( *dynamico*    ) atm=ico ; atmDomainType=unstructured  ;;
@@ -190 +207 @@
 esac
 
-case ${OCE} in
-    ( ORCA2.3*         ) OcePerio=4 ;;
-    ( ORCA1* | eORCA1* ) OcePerio=6 ;;
-    ( ORCA025*         ) OcePerio=6 ;;
+case ${OCE} in # Periodicity type of ORCA grid
+    ( ORCA2.3*             ) OcePerio=4 ;;
+    ( ORCA1* | eORCA1*     ) OcePerio=6 ;;
+    ( ORCA025* | eORCA025  ) OcePerio=6 ;;
 esac
 #
@@ -210 +227 @@
 
 ncks --overwrite --fl_fmt=${FMT_OASIS} --history ${OCE}_coordinates_mask.nc ${OCE}_coordinates_mask_${FMT_OASIS}.nc
-ncks --overwrite --fl_fmt=${FMT_OASIS} --history ${ATM}_grid.nc                    ${ATM}_grid_${FMT_OASIS}.nc
+ncks --overwrite --fl_fmt=${FMT_OASIS} --history ${ATM}_grid.nc             ${ATM}_grid_${FMT_OASIS}.nc
 #
-# Creates OCEAN C grid : redundant point removed to compute proper integrals # A passer dans CreateWeights
+# Creates OCEAN C grid : redundant points removed to compute proper integrals # A passer dans CreateWeights
 # --------------------------------------------------------------------------------------------------------
 cat <<EOF >add_c_grid.nco
@@ -296 +313 @@
 
 ##
-##  Creates mask on ATM grid
+##  Creates ocean fractions on ATM grid
 ## ===========================================================================
 cp dia_t${oce}_to_t${atm}_${Suffix}.nc  dia_t${oce}_to_t${atm}_${Suffix}_mask.nc
@@ -342 +359 @@
 ## Creates mask of coastal OCE points
 ## ===========================================================================
-python ComputeNemoCoast.py -n ${OcePerio} -i ${OCE}_coordinates_mask.nc # Creates OceCoastal
+ccc_mprun -n 1 python -u ComputeNemoCoast.py -n ${OcePerio} -i ${OCE}_coordinates_mask.nc # Creates OceCoastal
 ##
 ## Creates mask of coastal ATM points
@@ -357 +374 @@
 
 ##
-## Other weigths files
+## Other weights files
 ## ===========================================================================
 # Loop on commands
@@ -480 +497 @@
 NCO="$(ncks --version |& tail -1|sed 's/ncks //')"
 PYTHON_VER=$( python -c "import sys ; print (sys.version.split(' ')[0])" )
-for InFile in *${oce}_to_*${atm}_*.nc *${atm}_to_*${oce}_*.nc ${ATM}_grid_maskFrom_${OCE}.nc ${ATM}_grid_maskFrom_${OCE}_${FMT_OASIS}.nc; do
+for InFile in $(ls *${oce}_to_*${atm}_*.nc *${atm}_to_*${oce}_*.nc ${ATM}_grid_maskFrom_${OCE}.nc ${ATM}_grid_maskFrom_${OCE}_${FMT_OASIS}.nc 2>/dev/null) ; do
     ncatted --history \
             --attribute uuid,global,d,,                                           \
@@ -516 +533 @@
 done
 ##
-## Update and complete weights file to fit OASIS requested format
+## Update and complete weight files to fit OASIS requested format
 ## ===========================================================================
 cat <<EOF > add_dim.nco
@@ -567 +584 @@
 ## Add missing variables in rmp files
 ## ===========================================================================
-for rmpFile in rmp_?${atm}_to_[tuv]${oce}_*.nc rmp_[tuv]${oce}_to_t${atm}_*.nc* ; do
+for $(ls rmpFile in rmp_?${atm}_to_[tuv]${oce}_*.nc rmp_[tuv]${oce}_to_t${atm}_*.nc* 2>/dev/null) ; do
     echo ${rmpFile}
     a_to_o=false ; o_to_a=false
@@ -704 +721 @@
 
 ##
-## Creates and save auxiliary files for OASIS
-## ===========================================================================
+## Creates and save auxiliary files for OASIS : grids.nc, areas.nc and masks.nc
+## ============================================================================
 bash ${SUBMIT_DIR}/CreateOasisGrids.bash --oce ${OCE} --atm ${ATM}
 
@@ -713 +730 @@
 
 ##
-## Calving
+## Calving weights
 ## ===========================================================================
 case ${OCE} in
+    ( eORCA025 )
+    ccc_mprun -n 1 python -u CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=full     --dir=.
+    ;;
+    
     ( eORCA1.2 )
     cp ${R_IN}/OCE/NEMO/ORCA1_LIM3_PISCES/v3.6_stable/runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc .
-    for Type in nosouth iceberg iceshelf ; do
-        python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=${Type} --dir=. --isf_icb=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc
-    done
+    ccc_mprun -n 1 python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=nosouth  --dir=. 
+    ccc_mprun -n 1 python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=iceberg  --dir=. --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=Icb_flux
+    ccc_mprun -n 1 python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=iceshelf --dir=. --repartition_file=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc --repartition_var=sornfisf
     ;;
+    
     ( * )
-    python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=full --dir=.
+    ccc_mprun -n 1 python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=full     --dir=.
     ;;
 esac
@@ -780 +802 @@
 echo '$HeadURL$ ' >> README.txt
 
+# Compute checksums and add thmer to README
+# -----------------------------------------
 cat << EOF >> README.txt
 UUID common to all files : ${UUID}
 
-Files produced, with checksum produce by Unix command shasum (version $(shasum --version)) with default algorithm
+Files produced, with checksum produced by Unix command shasum (version $(shasum --version)) with default algorithm
 
 EOF