Changeset 4159


Ignore:
Timestamp:
12/03/18 11:19:11 (3 years ago)
Author:
omamce
Message:

O.M. : correct eORC025 case for calving weights

Location:
TOOLS/MOSAIX
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • TOOLS/MOSAIX/CalvingWeights.py

    r4146 r4159  
    3333parser = argparse.ArgumentParser ( 
    3434    description = """Compute calving weights""", 
    35     epilog='-------- This is the end of the help message --------') 
     35    epilog='-------- End of the help message --------') 
    3636 
    3737# Adding arguments 
    38 parser.add_argument ('--oce'   , help='oce model name', type=str, default='eORCA1.2', choices=['ORCA2.3', 'eORCA1.2', 'eORCA025'] ) 
    39 parser.add_argument ('--atm'   , help='atm model name (ICO* or LMD*)', type=str, default='ICO40'    ) 
    40 parser.add_argument ('--type'  , help='Type of distribution', type=str, choices=['iceshelf', 'iceberg', 'nosouth', 'full'], default='full'  ) 
    41 parser.add_argument ('--dir'   , help='Directory of inout file', type=str, default='./' ) 
    42 parser.add_argument ('--isf_icb', help='Data files with iceberg and iceshelf melting', type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' ) 
     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) 
    4344 
    4445# Parse command line 
     
    4647 
    4748# Model Names 
    48 src_Name = myargs.atm 
    49 dst_Name = myargs.oce 
     49src_Name = myargs.atm ; dst_Name = myargs.oce 
     50     
     51# Default vars 
     52if myargs.repartition_var == None : 
     53    # Runoff from Antarctica iceshelves (Depoorter, 2013) 
     54    if myargs.type == 'iceshelf' : repartitionVar = 'sornfisf' 
     55         
     56    # Runoff from Antarctica iceberg (Depoorter, 2013) 
     57    if myargs.type == 'iceberg' :  repartitionVar = 'Icb_Flux' 
     58         
     59else :                             repartitionVar = myargs.repartition_var 
    5060 
    5161# Latitude limits of each calving zone 
    5262limit_lat = ( (40.0, 90.0), (-50.0, 40.0), ( -90.0, -50.0) ) 
     63 
    5364nb_zone = len(limit_lat) 
    54  
    5565### ========================================================================================== 
    5666 
     
    7989f_grids = netCDF4.Dataset ( grids, "r" ) 
    8090 
    81 src_lon = f_grids.variables [ 't' + src_name + '.lon' ][:] 
    82 src_lat = f_grids.variables [ 't' + src_name + '.lat' ][:] 
    83 dst_lon = f_grids.variables [ 't' + dst_name + '.lon' ][:] 
    84 dst_lat = f_grids.variables [ 't' + dst_name + '.lat' ][:] 
    85  
    86 src_msk = f_masks.variables [ 't' + src_name + '.msk' ][:] 
    87 dst_msk = f_masks.variables [ 't' + dst_name + '.msk' ][:] 
     91src_lon = f_grids.variables ['t' + src_name + '.lon'][:] 
     92src_lat = f_grids.variables ['t' + src_name + '.lat'][:] 
     93dst_lon = f_grids.variables ['t' + dst_name + '.lon'][:] 
     94dst_lat = f_grids.variables ['t' + dst_name + '.lat'][:] 
     95 
     96src_msk = f_masks.variables ['t' + src_name + '.msk'][:] 
     97dst_msk = f_masks.variables ['t' + dst_name + '.msk'][:] 
    8898dst_mskutil = 1-dst_msk # Reversing the convention : 0 on continent, 1 on ocean 
    89  
    90 # Periodicityt masking for NEMO 
     99print ('dst_msk     : ' + str(np.sum(dst_msk))) 
     100print ('dst_mskutil : ' + str(np.sum(dst_mskutil))) 
     101 
     102# Periodicity masking for NEMO 
    91103if dst_Name == 'ORCA2.3'  : nperio_dst = 4 
    92104if dst_Name == 'eORCA1.2' : nperio_dst = 6 
     105if dst_Name == 'ORCA025'  : nperio_dst = 6 
    93106if dst_Name == 'eORCA025' : nperio_dst = 6 
     107print ('nperio_dst: ' + str(nperio_dst) ) 
    94108dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T' ) 
    95  
    96 # Fill Closed seas : preparation by closing some straits 
     109print ('dst_mskutil : ' + str(np.sum(dst_mskutil))) 
     110 
     111## 
     112## Fill Closed seas and other 
     113## ================================================== 
     114 
     115# Preparation by closing some straits 
     116# ----------------------------------- 
     117if dst_Name == 'eORCA025' : 
     118    print ('Filling some seas for eORCA025') 
     119    # Set Gibraltar strait to 0 to fill Mediterranean sea 
     120    dst_mskutil[838:841,1125] = 0 
     121    # Set Bal-El-Mandeb to 0 to fill Red Sea 
     122    dst_mskutil[736,1321:1324] = 0 
     123    # Set Stagerak to 0 to fill Baltic Sea 
     124    dst_mskutil[961,1179:1182] = 0 
     125    # Set Ormuz Strait to 0 to fill Arabo-Persian Gulf 
     126    dst_mskutil[794:798,1374] = 0 
     127    # Set Hudson Strait to 0 to fill Hudson Bay 
     128    dst_mskutil[997:1012,907] = 0 
     129     
    97130if dst_Name == 'eORCA1.2' : 
     131    print ('Filling some seas for eORCA1.2') 
    98132    # Set Gibraltar strait to 0 to fill Mediterrean sea 
    99133    dst_mskutil[240, 283] = 0 
     
    108142     
    109143if dst_Name == 'ORCA2.3' : 
     144    print ('Filling some seas for ORCA2.3') 
    110145    # Set Gibraltar strait to 0 to fill Mediterrean sea 
    111146    dst_mskutil[101,139] = 0 
     
    120155    dst_mskutil[87:89,166] = 0 
    121156 
    122 # Fill closed sea with image processing library 
     157dst_closed = dst_mskutil 
     158 
     159# Fill closed seas with image processing library 
     160# ---------------------------------------------- 
    123161dst_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' ) 
    124162 
    125163# Surfaces 
    126 src_srf = f_areas.variables [ 't' + src_name + '.srf' ] 
    127 dst_srf = f_areas.variables [ 't' + dst_name + '.srf' ] 
     164src_srf = f_areas.variables ['t' + src_name + '.srf'] 
     165dst_srf = f_areas.variables ['t' + dst_name + '.srf'] 
    128166dst_srfutil = dst_srf * np.float64 (dst_mskutil) 
    129167 
    130168dst_srfutil_sum = np.sum( dst_srfutil) 
    131169 
    132 src_clo = f_grids.variables [ 't' + src_name + '.clo' ][:] 
    133 src_cla = f_grids.variables [ 't' + src_name + '.cla' ][:] 
    134 dst_clo = f_grids.variables [ 't' + dst_name + '.clo' ][:] 
    135 dst_cla = f_grids.variables [ 't' + dst_name + '.cla' ][:] 
     170src_clo = f_grids.variables ['t' + src_name + '.clo'][:] 
     171src_cla = f_grids.variables ['t' + src_name + '.cla'][:] 
     172dst_clo = f_grids.variables ['t' + dst_name + '.clo'][:] 
     173dst_cla = f_grids.variables ['t' + dst_name + '.cla'][:] 
    136174 
    137175# Indices 
    138176( src_jpj, src_jpi) = src_lat.shape ; src_grid_size = src_jpj*src_jpi 
    139177( dst_jpj, dst_jpi) = dst_lat.shape ; dst_grid_size = dst_jpj*dst_jpi 
    140 orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 ## Fortran indices (starting at one) 
     178orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 # Fortran indices (starting at one) 
    141179 
    142180### ===== Reading needed data ================================================== 
    143181if 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 ) 
     182    # Reading data file for calving or iceberg geometry around Antarctica 
     183    print ( 'Opening ' + myargs.repartition_file) 
     184    f_repartition = netCDF4.Dataset ( myargs.repartition_file, "r" ) 
     185    repartition = np.sum ( f_repartition.variables [repartitionVar][:], axis=0 ) 
    153186 
    154187## Before loop on basins 
     
    160193 
    161194### ===== Starting loop on basins============================================== 
    162 ## Initialise some fields 
     195 
     196# Initialise some fields 
    163197remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) 
    164198src_address  = np.empty ( shape=(0), dtype=np.int32   ) 
    165199dst_address  = np.empty ( shape=(0), dtype=np.int32   ) 
     200 
     201basin_msk       = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float32) 
     202key_repartition = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float64) 
    166203 
    167204## Loop 
     
    175212    if myargs.type == 'iceberg'  and     south : ok_run = True  ; print ('Applying iceberg to south' ) 
    176213    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') 
     214    if myargs.type == 'iceberg'  and not south : ok_run = False ; print ('Skipping iceberg: not south ') 
     215    if myargs.type == 'iceshelf' and not south : ok_run = False ; print ('Skipping iceshelf: not south ') 
     216    if myargs.type == 'nosouth'  and     south : ok_run = False ; print ('Skipping south: nosouth case' ) 
     217    if myargs.type == 'nosouth'  and not south : ok_run = True  ; print ('Running: not in south, uniform repartition') 
    181218    if myargs.type == 'full'                   : ok_run = True  ; print ('Running general trivial case, uniform repartition' ) 
    182219 
    183220    if ok_run : 
     221        # Calving integral send to one point per latitude bands 
    184222        index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1 
    185223    
    186224        # Basin mask 
    187         basin_msk = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 ) 
     225        basin_msk[n_bas] = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 ) 
    188226 
    189227        # Repartition pattern 
    190         if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition = repartition * basin_msk 
    191         else :                                       key_repartition =               basin_msk 
     228        if myargs.type in ['iceberg', 'iceshelf' ] : key_repartition[n_bas] = repartition * basin_msk[n_bas] 
     229        else :                                       key_repartition[n_bas] =               basin_msk[n_bas] 
    192230 
    193231        # Integral and normalisation 
    194         sum_repartition = np.sum ( key_repartition * dst_srfutil ) 
     232        sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil ) 
    195233        key_repartition = key_repartition / sum_repartition 
    196234     
    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 )) ) 
     235        print ( 'Sum of repartition key                      : {:12.3e}'.format (np.sum (key_repartition[n_bas]               )) ) 
     236        print ( 'Integral (area weighted) of repartition key : {:12.3e}'.format (np.sum (key_repartition[n_bas] * dst_srfutil )) ) 
    199237     
    200238        # Basin surface (modulated by repartition key) 
    201         basin_srfutil = np.sum ( key_repartition * dst_srfutil ) 
     239        basin_srfutil = np.sum ( key_repartition[n_bas] * dst_srfutil ) 
    202240             
    203241        # Weights and links 
    204242        poids = 1.0 / basin_srfutil 
    205         matrix_local   = np.where ( basin_msk.ravel() > 0.5, key_repartition.ravel()*poids, 0. ) 
     243        matrix_local   = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids, 0. ) 
    206244        matrix_local   = matrix_local[matrix_local > 0.0] # Keep only non zero values 
    207245        num_links      = np.count_nonzero (matrix_local) 
    208         # address on source grid : all links points to the same LMDZ point. 
     246        # Address on source grid : all links points to the same LMDZ point. 
    209247        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) 
     248        # Address on destination grid : each NEMO point with non zero link 
     249        dst_address_local = np.where ( key_repartition[n_bas].ravel() > 0.0, orc_index, 0) 
    212250        dst_address_local = dst_address_local[dst_address_local > 0] 
    213251        # Append to global tabs 
     
    241279f_calving.Ocean           = dst_Name + " https://www.nemo-ocean.eu" 
    242280f_calving.Atmosphere      = src_Name + " http://lmdz.lmd.jussieu.fr" 
    243 if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.isf_icb 
     281if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.repartition_file 
    244282f_calving.associatedFiles = grids + " " + areas + " " + masks 
    245283f_calving.directory       = os.getcwd () 
     
    267305f_calving.SVN_HeadURL     = "$HeadURL$" 
    268306 
    269 num_links        = f_calving.createDimension ('num_links'       , num_links ) 
    270 num_wgts         = f_calving.createDimension ('num_wgts'        ,         1 ) 
    271  
    272 src_grid_size    = f_calving.createDimension ('src_grid_size'   , src_grid_size ) 
    273 src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0]  ) 
    274 src_grid_rank    = f_calving.createDimension ('src_grid_rank'   ,        2  ) 
    275  
    276 dst_grid_size    = f_calving.createDimension ('dst_grid_size'   , dst_grid_size ) 
    277 dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] ) 
    278 dst_grid_rank    = f_calving.createDimension ('dst_grid_rank'   ,        2  ) 
     307d_nb_zone          = f_calving.createDimension ('nb_zone'         ,   nb_zone ) 
     308d_num_links        = f_calving.createDimension ('num_links'       , num_links ) 
     309d_num_wgts         = f_calving.createDimension ('num_wgts'        ,         1 ) 
     310 
     311d_src_grid_size    = f_calving.createDimension ('src_grid_size'   , src_grid_size ) 
     312d_src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0]  ) 
     313D_src_grid_rank    = f_calving.createDimension ('src_grid_rank'   ,        2  ) 
     314 
     315d_dst_grid_size    = f_calving.createDimension ('dst_grid_size'   , dst_grid_size ) 
     316d_dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] ) 
     317d_dst_grid_rank    = f_calving.createDimension ('dst_grid_rank'   ,        2  ) 
    279318 
    280319v_remap_matrix = f_calving.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') ) 
     
    289328v_src_grid_center_lon = f_calving.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) ) 
    290329v_src_grid_center_lat = f_calving.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) ) 
    291 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"  
    292 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"  
     330v_src_grid_center_lon.units='degrees_east'  ; v_src_grid_center_lon.long_name='Longitude' 
     331v_src_grid_center_lon.long_name='longitude' ; v_src_grid_center_lon.bounds="src_grid_corner_lon"  
     332v_src_grid_center_lat.units='degrees_north' ; v_src_grid_center_lat.long_name='Latitude' 
     333v_src_grid_center_lat.long_name='latitude ' ; v_src_grid_center_lat.bounds="src_grid_corner_lat"  
    293334v_src_grid_corner_lon = f_calving.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners')  ) 
    294335v_src_grid_corner_lat = f_calving.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners')  ) 
     
    303344v_src_grid_center_lon[:] = src_lon[:].ravel() 
    304345v_src_grid_center_lat[:] = src_lat[:].ravel() 
    305 v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, src_grid_corners.__len__()) ) 
    306 v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, src_grid_corners.__len__()) ) 
     346v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) 
     347v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) 
    307348v_src_grid_area      [:] = src_srf[:].ravel() 
    308349v_src_grid_imask     [:] = src_msk[:].ravel() 
     
    313354v_dst_grid_center_lon = f_calving.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) ) 
    314355v_dst_grid_center_lat = f_calving.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) ) 
    315 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"  
    316 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"  
     356v_dst_grid_center_lon.units='degrees_east'  ; v_dst_grid_center_lon.long_name='Longitude' 
     357v_dst_grid_center_lon.long_name='longitude' ; v_dst_grid_center_lon.bounds="dst_grid_corner_lon"  
     358v_dst_grid_center_lat.units='degrees_north' ; v_dst_grid_center_lat.long_name='Latitude' 
     359v_dst_grid_center_lat.long_name='latitude'  ; v_dst_grid_center_lat.bounds="dst_grid_corner_lat"  
    317360v_dst_grid_corner_lon = f_calving.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners')  ) 
    318361v_dst_grid_corner_lat = f_calving.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners')  ) 
     
    324367v_dst_grid_imask.long_name="Land-sea mask" ; v_dst_grid_imask.units="Land:1, Ocean:0" 
    325368 
     369v_dst_bassin      = f_calving.createVariable ( 'dst_bassin'      , 'f8', ('nb_zone', 'dst_grid_size',)  ) 
     370v_dst_repartition = f_calving.createVariable ( 'dst_repartition' , 'f8', ('nb_zone', 'dst_grid_size',)  ) 
     371 
    326372v_dst_grid_dims      [:] = ( dst_jpi, dst_jpi )  
    327373v_dst_grid_center_lon[:] = dst_lon[:].ravel() 
    328374v_dst_grid_center_lat[:] = dst_lat[:].ravel() 
    329 v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, dst_grid_corners.__len__()) ) 
    330 v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, dst_grid_corners.__len__()) ) 
     375v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) 
     376v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) 
    331377v_dst_grid_area      [:] = dst_srf[:].ravel() 
    332378v_dst_grid_imask     [:] = dst_msk[:].ravel() 
    333379 
    334 # For diags 
     380v_dst_bassin[:]      = basin_msk.reshape      ( (nb_zone,dst_grid_size) ) 
     381v_dst_repartition[:] = key_repartition.reshape( (nb_zone,dst_grid_size) ) 
     382 
     383# Additoonnal fields for debugging 
     384# ================================ 
     385v_dst_grid_imaskutil      = f_calving.createVariable ( 'dst_grid_imaskutil'     , 'f8', ('dst_grid_size',)  ) 
     386v_dst_grid_imaskutil.long_name="Land-sea mask" ; v_dst_grid_imaskutil.units="Land:1, Ocean:0" 
     387v_dst_grid_imaskclose      = f_calving.createVariable ( 'dst_grid_imaskclose'     , 'f8', ('dst_grid_size',)  ) 
     388v_dst_grid_imaskclose.long_name="Land-sea mask" ; v_dst_grid_imaskclose.units="Land:1, Ocean:0" 
     389v_dst_grid_imaskutil [:] = dst_mskutil[:].ravel() 
     390v_dst_grid_imaskclose[:] = dst_closed[:].ravel() 
     391 
    335392v_dst_lon_addressed = f_calving.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) ) 
    336393v_dst_lat_addressed = f_calving.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) ) 
    337394v_dst_lon_addressed.long_name="Longitude" ; v_dst_lon_addressed.standard_name="longitude" ; v_dst_lon_addressed.units="degrees_east" 
    338395v_dst_lat_addressed.long_name="Latitude"  ; v_dst_lat_addressed.standard_name="latitude"  ; v_dst_lat_addressed.units="degrees_north" 
    339 v_dst_lon_addressed [:] = dst_lon.ravel()[dst_address-1].ravel() 
    340 v_dst_lat_addressed [:] = dst_lat.ravel()[dst_address-1].ravel() 
     396v_dst_area_addressed  = f_calving.createVariable ( 'dst_area_addressed', 'f8', ('num_links',) ) 
     397v_dst_area_addressed.long_name="Grid area" ; v_dst_area_addressed.standard_name="cell_area" ; v_dst_grid_area.units="m2" 
     398v_dst_imask_addressed = f_calving.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) ) 
     399v_dst_imask_addressed.long_name="Land-sea mask" ; v_dst_imask_addressed.units="Land:1, Ocean:0" 
     400v_dst_imaskutil_addressed = f_calving.createVariable ( 'dst_imaskutil_addressed', 'i4', ('num_links',) ) 
     401v_dst_imaskutil_addressed.long_name="Land-sea mask" ; v_dst_imaskutil_addressed.units="Land:1, Ocean:0" 
     402 
     403v_dst_lon_addressed  [:] = dst_lon.ravel()[dst_address-1].ravel() 
     404v_dst_lat_addressed  [:] = dst_lat.ravel()[dst_address-1].ravel() 
     405v_dst_area_addressed [:] = dst_srf[:].ravel()[dst_address-1].ravel() 
     406v_dst_imask_addressed[:] = dst_msk[:].ravel()[dst_address-1].ravel() 
     407v_dst_imaskutil_addressed[:] = dst_mskutil.ravel()[dst_address-1].ravel() 
    341408 
    342409f_calving.close () 
  • TOOLS/MOSAIX/CreateWeightsMask.bash

    r4153 r4159  
    44#MSUB -e Out_WeightsMask    # Error output 
    55#MSUB -eo 
    6 #MSUB -n 8                  # Number of processors 
     6#MSUB -n 16                 # Number of processors 
    77#MSUB -T 1800               # Time limit (seconds) 
    88#MSUB -q skylake 
     
    4949## Configuration 
    5050## =========================================================================== 
     51set -e 
    5152 
    5253# 
     
    5455# ============== 
    5556#OCE=ORCA2.3 
    56 OCE=eORCA1.2 
    57 #OCE=ORCA025 
     57#OCE=eORCA1.2 
     58OCE=eORCA025 
    5859 
    5960#ATM=ICO30 
    6061#ATM=ICO40 
    6162#ATM=ICO450 
    62 #ATM=LMD9695 
    63 ATM=LMD144142 
     63ATM=LMD9695 
     64#ATM=LMD144142 
    6465 
    6566# Default values, used to create ocean fraction on atmospheric grid  
    66 DefaultValues=( Direction=o2a,ocegrid=t,atmgrid=t,Order=1st,Quantity=false,Renormalize=false,useArea=false,maskSrc=true,maskDst=false ) 
     67DefaultValues=( Direction=o2a,oceGrid=t,atmGrid=t,Order=1st,Quantity=false,Renormalize=false,useArea=false,maskSrc=true,maskDst=false ) 
    6768 
    6869# List of weights to build 
    6970CommandList=( \ 
    70     Direction=o2a,Order=1st,Quantity=false,Renormalize=true,ocegrid=t,atmgrid=t,useArea=false,maskSrc=true,maskDst=true  \ 
    71     Direction=a2o,Order=1st,Quantity=false,Renormalize=false,atmgrid=t,ocegrid=t,useArea=true,maskSrc=true,maskDst=true  \ 
    72     Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmgrid=t,ocegrid=u,useArea=true,maskSrc=true,maskDst=true  \ 
    73     Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmgrid=t,ocegrid=v,useArea=true,maskSrc=true,maskDst=true  \ 
    74     Direction=a2o,Order=1st,Quantity=true,Renormalize=false,atmgrid=c,ocegrid=t,useArea=true,maskSrc=true,maskDst=true   \ 
    75    ) 
     71    Direction=o2a,Order=1st,Quantity=false,Renormalize=true,oceGrid=t,atmGrid=t,useArea=false,maskSrc=true,maskDst=true  \ 
     72    Direction=a2o,Order=1st,Quantity=false,Renormalize=false,atmGrid=t,oceGrid=t,useArea=true,maskSrc=true,maskDst=true  \ 
     73    Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmGrid=t,oceGrid=u,useArea=true,maskSrc=true,maskDst=true  \ 
     74    Direction=a2o,Order=2nd,Quantity=false,Renormalize=false,atmGrid=t,oceGrid=v,useArea=true,maskSrc=true,maskDst=true  \ 
     75    Direction=a2o,Order=1st,Quantity=true,Renormalize=false,atmGrid=c,oceGrid=t,useArea=true,maskSrc=true,maskDst=true   \ 
     76    ) 
     77 
     78# Each item in CommandList describes the property of interpolation weights to generate.  
     79# White spaces separate analysis. No spaces in any analysis. 
     80# 
     81# Keywords :  
     82#   Direction   : o2a for ocean to atmosphere, a2o for atmosphere to ocean 
     83#   Order       : 1st or 2nd 
     84#   Quantity    : true if integrated quantity over a grid box, false for flux (quantity / m^2) 
     85#                 or intensive value (temperature, salinity, sea-ice fraction, ...) 
     86#   Renormalize : used when source grid is masked, to use values on non masked points only 
     87#   oceGrid     : t, u or v point for OPA grid 
     88#   atmGrid     : up to know, only t grid used in the atmosphere 
     89#   useArea     : if true area from the model metrics is used. If false, areas are computed by XIOS from grid corners 
     90#   maskSrc     : true to use the source grid mask, false to used all grid points 
     91#   maskDst     : true to use the destination grid mask, false to use all grid points 
     92# 
    7693 
    7794## =========================================================================== 
     
    8198## =========================================================================== 
    8299SUBMIT_DIR=$(pwd) 
    83 # Function to handle command parameters 
     100 
     101# Functions to handle command parameters 
     102# ====================================== 
    84103function read_Command { 
    85104    # Decipher the command line to set bash variables  
     
    92111    done 
    93112    local l_Command=${1} 
    94     for l_Element in $(echo ${l_Command} | tr "," "\n" ) 
    95     do 
     113    for l_Element in $(echo ${l_Command} | tr "," "\n" ) ; do 
    96114        [[ "X${l_Debug}" = "Xtrue" ]] && echo ${l_Element} 
    97115        eval export ${l_Element} 
     
    101119function setValues { 
    102120    # 
    103     ocegrid=${ocegrid,,} 
    104     atmgrid=${atmgrid,,} 
    105     OCEGRID=${ocegrid^^} 
    106     ATMGRID=${atmgrid^^} 
     121    oceGrid=${oceGrid,,} ; atmGrid=${atmGrid,,} 
     122    OCEGRID=${oceGrid^^} ; ATMGRID=${atmGrid^^} 
    107123 
    108124    case ${Order} in 
     
    127143    case ${Direction} in 
    128144        ( o2a ) 
    129         src=${oce} ; SRC=${OCE} ; srcGrid=${ocegrid} ; srcDomainType=${oceDomainType} ; SRCGRID=${OCEGRID} ; srcArea=area_grid_${OCEGRID} 
    130         dst=${atm} ; DST=${ATM} ; dstGrid=${atmgrid} ; dstDomainType=${atmDomainType} ; DSTGRID=${ATMGRID} ; dstArea=aire 
     145        src=${oce} ; SRC=${OCE} ; srcGrid=${oceGrid} ; srcDomainType=${oceDomainType} ; SRCGRID=${OCEGRID} ; srcArea=area_grid_${OCEGRID} 
     146        dst=${atm} ; DST=${ATM} ; dstGrid=${atmGrid} ; dstDomainType=${atmDomainType} ; DSTGRID=${ATMGRID} ; dstArea=aire 
    131147        ;;  
    132148        ( a2o ) 
    133         src=${atm} ; SRC=${ATM} ; srcGrid=${atmgrid} ; srcDomainType=${atmDomainType} ; SRCGRID=${ATMGRID} ; srcArea=aire 
    134         dst=${oce} ; DST=${OCE} ; dstGrid=${ocegrid} ; dstDomainType=${oceDomainType} ; DSTGRID=${OCEGRID} ; dstArea=area_grid_${OCEGRID} 
     149        src=${atm} ; SRC=${ATM} ; srcGrid=${atmGrid} ; srcDomainType=${atmDomainType} ; SRCGRID=${ATMGRID} ; srcArea=aire 
     150        dst=${oce} ; DST=${OCE} ; dstGrid=${oceGrid} ; dstDomainType=${oceDomainType} ; DSTGRID=${OCEGRID} ; dstArea=area_grid_${OCEGRID} 
    135151        ;; 
    136152    esac 
     
    166182    TMPDIR=${HOME}/Scratch/TMP 
    167183    SUBMIT_DIR=$(pwd) 
    168     MPIRUN=/opt/local/bin/mpirun-openmpi-gcc49 -n 2 
     184    MPIRUN=/opt/local/bin/mpirun -n 4 
    169185    ;;  
    170186    ( * ) exit -1 ;; 
     
    184200    ( *ORC* )         oce=orc ; oceDomainType=curvilinear   ;; 
    185201esac 
     202 
    186203case ${ATM} in 
    187204    ( *dynamico*    ) atm=ico ; atmDomainType=unstructured  ;; 
     
    190207esac 
    191208 
    192 case ${OCE} in 
    193     ( ORCA2.3*         ) OcePerio=4 ;; 
    194     ( ORCA1* | eORCA1* ) OcePerio=6 ;; 
    195     ( ORCA025*         ) OcePerio=6 ;; 
     209case ${OCE} in # Periodicity type of ORCA grid 
     210    ( ORCA2.3*             ) OcePerio=4 ;; 
     211    ( ORCA1* | eORCA1*     ) OcePerio=6 ;; 
     212    ( ORCA025* | eORCA025  ) OcePerio=6 ;; 
    196213esac 
    197214# 
     
    210227 
    211228ncks --overwrite --fl_fmt=${FMT_OASIS} --history ${OCE}_coordinates_mask.nc ${OCE}_coordinates_mask_${FMT_OASIS}.nc 
    212 ncks --overwrite --fl_fmt=${FMT_OASIS} --history ${ATM}_grid.nc                    ${ATM}_grid_${FMT_OASIS}.nc 
     229ncks --overwrite --fl_fmt=${FMT_OASIS} --history ${ATM}_grid.nc             ${ATM}_grid_${FMT_OASIS}.nc 
    213230# 
    214 # Creates OCEAN C grid : redundant point removed to compute proper integrals # A passer dans CreateWeights 
     231# Creates OCEAN C grid : redundant points removed to compute proper integrals # A passer dans CreateWeights 
    215232# -------------------------------------------------------------------------------------------------------- 
    216233cat <<EOF >add_c_grid.nco 
     
    296313 
    297314## 
    298 ##  Creates mask on ATM grid 
     315##  Creates ocean fractions on ATM grid 
    299316## =========================================================================== 
    300317cp dia_t${oce}_to_t${atm}_${Suffix}.nc  dia_t${oce}_to_t${atm}_${Suffix}_mask.nc 
     
    342359## Creates mask of coastal OCE points 
    343360## =========================================================================== 
    344 python ComputeNemoCoast.py -n ${OcePerio} -i ${OCE}_coordinates_mask.nc # Creates OceCoastal 
     361ccc_mprun -n 1 python -u ComputeNemoCoast.py -n ${OcePerio} -i ${OCE}_coordinates_mask.nc # Creates OceCoastal 
    345362## 
    346363## Creates mask of coastal ATM points 
     
    357374 
    358375## 
    359 ## Other weigths files 
     376## Other weights files 
    360377## =========================================================================== 
    361378# Loop on commands 
     
    480497NCO="$(ncks --version |& tail -1|sed 's/ncks //')" 
    481498PYTHON_VER=$( python -c "import sys ; print (sys.version.split(' ')[0])" ) 
    482 for InFile in *${oce}_to_*${atm}_*.nc *${atm}_to_*${oce}_*.nc ${ATM}_grid_maskFrom_${OCE}.nc ${ATM}_grid_maskFrom_${OCE}_${FMT_OASIS}.nc; do 
     499for 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 
    483500    ncatted --history \ 
    484501            --attribute uuid,global,d,,                                           \ 
     
    516533done 
    517534## 
    518 ## Update and complete weights file to fit OASIS requested format 
     535## Update and complete weight files to fit OASIS requested format 
    519536## =========================================================================== 
    520537cat <<EOF > add_dim.nco 
     
    567584## Add missing variables in rmp files 
    568585## =========================================================================== 
    569 for rmpFile in rmp_?${atm}_to_[tuv]${oce}_*.nc rmp_[tuv]${oce}_to_t${atm}_*.nc* ; do 
     586for $(ls rmpFile in rmp_?${atm}_to_[tuv]${oce}_*.nc rmp_[tuv]${oce}_to_t${atm}_*.nc* 2>/dev/null) ; do 
    570587    echo ${rmpFile} 
    571588    a_to_o=false ; o_to_a=false 
     
    704721 
    705722## 
    706 ## Creates and save auxiliary files for OASIS 
    707 ## =========================================================================== 
     723## Creates and save auxiliary files for OASIS : grids.nc, areas.nc and masks.nc 
     724## ============================================================================ 
    708725bash ${SUBMIT_DIR}/CreateOasisGrids.bash --oce ${OCE} --atm ${ATM} 
    709726 
     
    713730 
    714731## 
    715 ## Calving 
     732## Calving weights 
    716733## =========================================================================== 
    717734case ${OCE} in 
     735    ( eORCA025 ) 
     736    ccc_mprun -n 1 python -u CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=full     --dir=. 
     737    ;; 
     738     
    718739    ( eORCA1.2 ) 
    719740    cp ${R_IN}/OCE/NEMO/ORCA1_LIM3_PISCES/v3.6_stable/runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc . 
    720     for Type in nosouth iceberg iceshelf ; do 
    721         python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=${Type} --dir=. --isf_icb=runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc 
    722     done 
     741    ccc_mprun -n 1 python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=nosouth  --dir=.  
     742    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 
     743    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 
    723744    ;; 
     745     
    724746    ( * ) 
    725     python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=full --dir=. 
     747    ccc_mprun -n 1 python CalvingWeights.py --oce=${OCE} --atm=${ATM} --type=full    --dir=. 
    726748    ;; 
    727749esac 
     
    780802echo '$HeadURL$ ' >> README.txt 
    781803 
     804# Compute checksums and add thmer to README 
     805# ----------------------------------------- 
    782806cat << EOF >> README.txt 
    783807UUID common to all files : ${UUID} 
    784808 
    785 Files produced, with checksum produce by Unix command shasum (version $(shasum --version)) with default algorithm 
     809Files produced, with checksum produced by Unix command shasum (version $(shasum --version)) with default algorithm 
    786810 
    787811EOF 
Note: See TracChangeset for help on using the changeset viewer.