[4097] | 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 | ## |
---|
| 16 | import netCDF4, numpy as np |
---|
| 17 | import sys, os, platform, argparse, textwrap, time |
---|
| 18 | from scipy import ndimage |
---|
| 19 | import 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 |
---|
| 33 | parser = argparse.ArgumentParser ( |
---|
| 34 | description = """Compute calving weights""", |
---|
[4159] | 35 | epilog='-------- End of the help message --------') |
---|
[4097] | 36 | |
---|
| 37 | # Adding arguments |
---|
[4159] | 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 input file', type=str, default='./' ) |
---|
| 42 | parser.add_argument ('--repartition_file', help='Data files with iceberg melting' , type=str, default='./runoff-icb_DaiTrenberth_Depoorter_eORCA1_JD.nc' ) |
---|
| 43 | parser.add_argument ('--repartition_var' , help='Variable name for iceshelf' , type=str, default=None) |
---|
[4199] | 44 | parser.add_argument ('--output' , help='output rmp file name', default='rmp_tlmd_to_torc_calving.nc' ) |
---|
| 45 | parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) |
---|
[4097] | 46 | |
---|
| 47 | # Parse command line |
---|
| 48 | myargs = parser.parse_args() |
---|
| 49 | |
---|
| 50 | # Model Names |
---|
[4159] | 51 | src_Name = myargs.atm ; dst_Name = myargs.oce |
---|
| 52 | |
---|
| 53 | # Default vars |
---|
| 54 | if 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 | |
---|
| 61 | else : repartitionVar = myargs.repartition_var |
---|
[4097] | 62 | |
---|
| 63 | # Latitude limits of each calving zone |
---|
| 64 | limit_lat = ( (40.0, 90.0), (-50.0, 40.0), ( -90.0, -50.0) ) |
---|
[4159] | 65 | |
---|
[4097] | 66 | nb_zone = len(limit_lat) |
---|
[4172] | 67 | |
---|
| 68 | if myargs.fmt == 'classic' : FmtNetcdf = 'CLASSIC' |
---|
| 69 | if myargs.fmt == 'netcdf3' : FmtNetcdf = 'CLASSIC' |
---|
| 70 | if myargs.fmt == '64bit' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' |
---|
| 71 | if myargs.fmt == '64bit_data' : FmtNetcdf = 'NETCDF3_64BIT_DATA' |
---|
| 72 | if myargs.fmt == '64bit_offset' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' |
---|
| 73 | if myargs.fmt == 'netcdf4' : FmtNetcdf = 'NETCDF4' |
---|
| 74 | if myargs.fmt == 'netcdf4_classic' : FmtNetcdf = 'NETCDF4_CLASSIC' |
---|
[4097] | 75 | ### ========================================================================================== |
---|
| 76 | |
---|
| 77 | # Model short names |
---|
| 78 | src_name = None ; dst_name = None |
---|
| 79 | if src_Name.count('ICO') != 0 : src_name = 'ico' |
---|
| 80 | if src_Name.count('LMD') != 0 : src_name = 'lmd' |
---|
| 81 | if dst_Name.count('ORCA') != 0 : dst_name = 'orc' |
---|
| 82 | |
---|
| 83 | CplModel = dst_Name + 'x' + src_Name |
---|
| 84 | |
---|
| 85 | print ('Atm names : ' + src_name + ' ' + src_Name ) |
---|
| 86 | print ('Oce names : ' + dst_name + ' ' + dst_Name ) |
---|
| 87 | print (' ') |
---|
| 88 | |
---|
| 89 | # Reading domains characteristics |
---|
| 90 | areas = myargs.dir + '/areas_' + dst_Name + 'x' + src_Name + '.nc' |
---|
| 91 | masks = myargs.dir + '/masks_' + dst_Name + 'x' + src_Name + '.nc' |
---|
| 92 | grids = myargs.dir + '/grids_' + dst_Name + 'x' + src_Name + '.nc' |
---|
| 93 | |
---|
| 94 | print ('Opening ' + areas) |
---|
| 95 | f_areas = netCDF4.Dataset ( areas, "r" ) |
---|
| 96 | print ('Opening ' + masks) |
---|
| 97 | f_masks = netCDF4.Dataset ( masks, "r" ) |
---|
| 98 | print ('Opening ' + grids) |
---|
| 99 | f_grids = netCDF4.Dataset ( grids, "r" ) |
---|
| 100 | |
---|
[4159] | 101 | src_lon = f_grids.variables ['t' + src_name + '.lon'][:] |
---|
| 102 | src_lat = f_grids.variables ['t' + src_name + '.lat'][:] |
---|
| 103 | dst_lon = f_grids.variables ['t' + dst_name + '.lon'][:] |
---|
| 104 | dst_lat = f_grids.variables ['t' + dst_name + '.lat'][:] |
---|
[4097] | 105 | |
---|
[4159] | 106 | src_msk = f_masks.variables ['t' + src_name + '.msk'][:] |
---|
| 107 | dst_msk = f_masks.variables ['t' + dst_name + '.msk'][:] |
---|
[4097] | 108 | dst_mskutil = 1-dst_msk # Reversing the convention : 0 on continent, 1 on ocean |
---|
[4159] | 109 | print ('dst_msk : ' + str(np.sum(dst_msk))) |
---|
| 110 | print ('dst_mskutil : ' + str(np.sum(dst_mskutil))) |
---|
[4097] | 111 | |
---|
[4159] | 112 | # Periodicity masking for NEMO |
---|
[4097] | 113 | if dst_Name == 'ORCA2.3' : nperio_dst = 4 |
---|
| 114 | if dst_Name == 'eORCA1.2' : nperio_dst = 6 |
---|
[4159] | 115 | if dst_Name == 'ORCA025' : nperio_dst = 6 |
---|
[4097] | 116 | if dst_Name == 'eORCA025' : nperio_dst = 6 |
---|
[4159] | 117 | print ('nperio_dst: ' + str(nperio_dst) ) |
---|
[4097] | 118 | dst_mskutil = nemo.lbc_mask (dst_mskutil, nperio=nperio_dst, cd_type='T' ) |
---|
[4159] | 119 | print ('dst_mskutil : ' + str(np.sum(dst_mskutil))) |
---|
[4097] | 120 | |
---|
[4159] | 121 | ## |
---|
| 122 | ## Fill Closed seas and other |
---|
| 123 | ## ================================================== |
---|
| 124 | |
---|
| 125 | # Preparation by closing some straits |
---|
| 126 | # ----------------------------------- |
---|
| 127 | if 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 | |
---|
[4097] | 140 | if dst_Name == 'eORCA1.2' : |
---|
[4159] | 141 | print ('Filling some seas for eORCA1.2') |
---|
[4097] | 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 | |
---|
| 153 | if dst_Name == 'ORCA2.3' : |
---|
[4159] | 154 | print ('Filling some seas for ORCA2.3') |
---|
[4097] | 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 | |
---|
[4159] | 167 | dst_closed = dst_mskutil |
---|
| 168 | |
---|
| 169 | # Fill closed seas with image processing library |
---|
| 170 | # ---------------------------------------------- |
---|
[4097] | 171 | 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' ) |
---|
| 172 | |
---|
| 173 | # Surfaces |
---|
[4159] | 174 | src_srf = f_areas.variables ['t' + src_name + '.srf'] |
---|
| 175 | dst_srf = f_areas.variables ['t' + dst_name + '.srf'] |
---|
[4097] | 176 | dst_srfutil = dst_srf * np.float64 (dst_mskutil) |
---|
| 177 | |
---|
| 178 | dst_srfutil_sum = np.sum( dst_srfutil) |
---|
| 179 | |
---|
[4159] | 180 | src_clo = f_grids.variables ['t' + src_name + '.clo'][:] |
---|
| 181 | src_cla = f_grids.variables ['t' + src_name + '.cla'][:] |
---|
| 182 | dst_clo = f_grids.variables ['t' + dst_name + '.clo'][:] |
---|
| 183 | dst_cla = f_grids.variables ['t' + dst_name + '.cla'][:] |
---|
[4097] | 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 |
---|
[4159] | 188 | orc_index = np.arange (dst_jpj*dst_jpi, dtype=np.int32) + 1 # Fortran indices (starting at one) |
---|
[4097] | 189 | |
---|
| 190 | ### ===== Reading needed data ================================================== |
---|
| 191 | if myargs.type in ['iceberg', 'iceshelf' ]: |
---|
[4159] | 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 ) |
---|
[4097] | 196 | |
---|
| 197 | ## Before loop on basins |
---|
| 198 | remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) |
---|
| 199 | src_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 200 | dst_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 201 | |
---|
| 202 | print (' ') |
---|
| 203 | |
---|
| 204 | ### ===== Starting loop on basins============================================== |
---|
[4159] | 205 | |
---|
| 206 | # Initialise some fields |
---|
[4097] | 207 | remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) |
---|
| 208 | src_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 209 | dst_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 210 | |
---|
[4159] | 211 | basin_msk = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float32) |
---|
| 212 | key_repartition = np.zeros( shape=(nb_zone, dst_jpj, dst_jpi), dtype=np.float64) |
---|
| 213 | |
---|
[4097] | 214 | ## Loop |
---|
| 215 | for 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' ) |
---|
[4159] | 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') |
---|
[4097] | 228 | if myargs.type == 'full' : ok_run = True ; print ('Running general trivial case, uniform repartition' ) |
---|
| 229 | |
---|
| 230 | if ok_run : |
---|
[4159] | 231 | # Calving integral send to one point per latitude bands |
---|
[4097] | 232 | index_src = ((src_grid_size - 1)*n_bas) // (nb_zone-1) + 1 |
---|
| 233 | |
---|
| 234 | # Basin mask |
---|
[4159] | 235 | basin_msk[n_bas] = np.where ( (dst_lat > lat_south ) & (dst_lat <= lat_north ), dst_mskutil, 0 ) |
---|
[4097] | 236 | |
---|
| 237 | # Repartition pattern |
---|
[4159] | 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] |
---|
[4097] | 240 | |
---|
| 241 | # Integral and normalisation |
---|
[4159] | 242 | sum_repartition = np.sum ( key_repartition[n_bas] * dst_srfutil ) |
---|
[4097] | 243 | key_repartition = key_repartition / sum_repartition |
---|
| 244 | |
---|
[4159] | 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 )) ) |
---|
[4097] | 247 | |
---|
| 248 | # Basin surface (modulated by repartition key) |
---|
[4159] | 249 | basin_srfutil = np.sum ( key_repartition[n_bas] * dst_srfutil ) |
---|
[4097] | 250 | |
---|
| 251 | # Weights and links |
---|
| 252 | poids = 1.0 / basin_srfutil |
---|
[4159] | 253 | matrix_local = np.where ( basin_msk[n_bas].ravel() > 0.5, key_repartition[n_bas].ravel()*poids, 0. ) |
---|
[4097] | 254 | matrix_local = matrix_local[matrix_local > 0.0] # Keep only non zero values |
---|
| 255 | num_links = np.count_nonzero (matrix_local) |
---|
[4159] | 256 | # Address on source grid : all links points to the same LMDZ point. |
---|
[4097] | 257 | src_address_local = np.ones(num_links, dtype=np.int32 )*index_src |
---|
[4159] | 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) |
---|
[4097] | 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 | |
---|
[4172] | 270 | |
---|
| 271 | |
---|
[4097] | 272 | ## End of loop |
---|
| 273 | print (' ') |
---|
| 274 | |
---|
| 275 | # |
---|
| 276 | num_links = np.count_nonzero (remap_matrix) |
---|
| 277 | print ( 'Total num_links : {:10d}'.format(num_links) ) |
---|
| 278 | |
---|
[4172] | 279 | # Diag : interpolate uniform field |
---|
| 280 | src_field = np.zeros(shape=(src_grid_size)) |
---|
| 281 | for 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 | |
---|
| 285 | dst_field = np.zeros(shape=(dst_grid_size,)) |
---|
| 286 | for link in np.arange (num_links) : |
---|
| 287 | dst_field [dst_address [link]-1] += remap_matrix[link] * src_field[src_address[link]-1] |
---|
| 288 | |
---|
[4097] | 289 | ### ===== Writing the weights file, for OASIS MCT ========================================== |
---|
| 290 | |
---|
| 291 | # Output file |
---|
[4172] | 292 | calving = myargs.output |
---|
| 293 | f_calving = netCDF4.Dataset ( calving, 'w', format=FmtNetcdf ) |
---|
| 294 | |
---|
[4097] | 295 | print ('Output file: ' + calving ) |
---|
| 296 | |
---|
| 297 | f_calving.Conventions = "CF-1.6" |
---|
| 298 | f_calving.source = "IPSL Earth system model" |
---|
| 299 | f_calving.group = "ICMC IPSL Climate Modelling Center" |
---|
| 300 | f_calving.Institution = "IPSL https.//www.ipsl.fr" |
---|
| 301 | f_calving.Ocean = dst_Name + " https://www.nemo-ocean.eu" |
---|
| 302 | f_calving.Atmosphere = src_Name + " http://lmdz.lmd.jussieu.fr" |
---|
[4159] | 303 | if myargs.type in ['iceberg', 'iceshelf' ]: f_calving.originalFiles = myargs.repartition_file |
---|
[4097] | 304 | f_calving.associatedFiles = grids + " " + areas + " " + masks |
---|
| 305 | f_calving.directory = os.getcwd () |
---|
| 306 | f_calving.description = "Generated with XIOS http://forge.ipsl.jussieu.fr/ioserver and MOSAIX https://forge.ipsl.jussieu.fr/igcmg/browser/TOOLS/MOSAIX" |
---|
| 307 | f_calving.title = calving |
---|
| 308 | f_calving.Program = "Generated by " + sys.argv[0] + " with flags " + str(sys.argv[1:]) |
---|
[4172] | 309 | f_calving.repartitionType = myargs.type |
---|
[4186] | 310 | if myargs.type in [ 'iceberg', 'iceshelf' ] : |
---|
[4172] | 311 | f_calving.repartitionFile = myargs.repartition_file |
---|
| 312 | f_calving.repartitionVar = repartitionVar |
---|
| 313 | f_calving.gridsFile = grids |
---|
| 314 | f_calving.areasFile = areas |
---|
| 315 | f_calving.masksFile = masks |
---|
[4097] | 316 | f_calving.timeStamp = time.asctime() |
---|
| 317 | f_calving.uuid = f_areas.uuid |
---|
| 318 | f_calving.HOSTNAME = platform.node() |
---|
| 319 | #f_calving.LOGNAME = os.getlogin() |
---|
| 320 | f_calving.Python = "Python version " + platform.python_version() |
---|
| 321 | f_calving.OS = platform.system() |
---|
| 322 | f_calving.release = platform.release() |
---|
| 323 | f_calving.hardware = platform.machine() |
---|
| 324 | f_calving.conventions = "SCRIP" |
---|
| 325 | if src_name == 'lmd' : f_calving.source_grid = "curvilinear" |
---|
| 326 | if src_name == 'ico' : f_calving.source_grid = "unstructured" |
---|
| 327 | f_calving.dest_grid = "curvilinear" |
---|
| 328 | f_calving.Model = "IPSL CM6" |
---|
| 329 | f_calving.SVN_Author = "$Author$" |
---|
| 330 | f_calving.SVN_Date = "$Date$" |
---|
| 331 | f_calving.SVN_Revision = "$Revision$" |
---|
| 332 | f_calving.SVN_Id = "$Id$" |
---|
| 333 | f_calving.SVN_HeadURL = "$HeadURL$" |
---|
| 334 | |
---|
[4186] | 335 | d_nb_zone = f_calving.createDimension ('nb_zone' , nb_zone ) |
---|
| 336 | d_num_links = f_calving.createDimension ('num_links' , num_links ) |
---|
| 337 | d_num_wgts = f_calving.createDimension ('num_wgts' , 1 ) |
---|
[4097] | 338 | |
---|
[4159] | 339 | d_src_grid_size = f_calving.createDimension ('src_grid_size' , src_grid_size ) |
---|
| 340 | d_src_grid_corners = f_calving.createDimension ('src_grid_corners', src_clo.shape[0] ) |
---|
[4186] | 341 | d_src_grid_rank = f_calving.createDimension ('src_grid_rank' , 2 ) |
---|
[4097] | 342 | |
---|
[4159] | 343 | d_dst_grid_size = f_calving.createDimension ('dst_grid_size' , dst_grid_size ) |
---|
| 344 | d_dst_grid_corners = f_calving.createDimension ('dst_grid_corners', dst_clo.shape[0] ) |
---|
| 345 | d_dst_grid_rank = f_calving.createDimension ('dst_grid_rank' , 2 ) |
---|
[4097] | 346 | |
---|
| 347 | v_remap_matrix = f_calving.createVariable ( 'remap_matrix', 'f8', ('num_links', 'num_wgts') ) |
---|
| 348 | v_src_address = f_calving.createVariable ( 'src_address' , 'i4', ('num_links',) ) |
---|
| 349 | v_src_address = f_calving.createVariable ( 'dst_address' , 'i4', ('num_links',) ) |
---|
| 350 | |
---|
| 351 | v_remap_matrix[:] = remap_matrix |
---|
| 352 | v_src_address [:] = src_address |
---|
| 353 | v_src_address [:] = src_address |
---|
| 354 | |
---|
| 355 | v_src_grid_dims = f_calving.createVariable ( 'src_grid_dims' , 'i4', ('src_grid_rank',) ) |
---|
| 356 | v_src_grid_center_lon = f_calving.createVariable ( 'src_grid_center_lon', 'i4', ('src_grid_size',) ) |
---|
| 357 | v_src_grid_center_lat = f_calving.createVariable ( 'src_grid_center_lat', 'i4', ('src_grid_size',) ) |
---|
[4159] | 358 | v_src_grid_center_lon.units='degrees_east' ; v_src_grid_center_lon.long_name='Longitude' |
---|
| 359 | v_src_grid_center_lon.long_name='longitude' ; v_src_grid_center_lon.bounds="src_grid_corner_lon" |
---|
| 360 | v_src_grid_center_lat.units='degrees_north' ; v_src_grid_center_lat.long_name='Latitude' |
---|
| 361 | v_src_grid_center_lat.long_name='latitude ' ; v_src_grid_center_lat.bounds="src_grid_corner_lat" |
---|
[4097] | 362 | v_src_grid_corner_lon = f_calving.createVariable ( 'src_grid_corner_lon', 'f8', ('src_grid_size', 'src_grid_corners') ) |
---|
| 363 | v_src_grid_corner_lat = f_calving.createVariable ( 'src_grid_corner_lat', 'f8', ('src_grid_size', 'src_grid_corners') ) |
---|
| 364 | v_src_grid_corner_lon.units="degrees_east" |
---|
| 365 | v_src_grid_corner_lat.units="degrees_north" |
---|
| 366 | v_src_grid_area = f_calving.createVariable ( 'src_grid_area' , 'f8', ('src_grid_size',) ) |
---|
| 367 | v_src_grid_area.long_name="Grid area" ; v_src_grid_area.standard_name="cell_area" ; v_src_grid_area.units="m2" |
---|
| 368 | v_src_grid_imask = f_calving.createVariable ( 'src_grid_imask' , 'f8', ('src_grid_size',) ) |
---|
| 369 | v_src_grid_imask.long_name="Land-sea mask" ; v_src_grid_imask.units="Land:1, Ocean:0" |
---|
| 370 | |
---|
| 371 | v_src_grid_dims [:] = ( src_jpi, src_jpi ) |
---|
| 372 | v_src_grid_center_lon[:] = src_lon[:].ravel() |
---|
| 373 | v_src_grid_center_lat[:] = src_lat[:].ravel() |
---|
[4159] | 374 | v_src_grid_corner_lon[:] = src_clo.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) |
---|
| 375 | v_src_grid_corner_lat[:] = src_cla.reshape( (src_jpi*src_jpj, d_src_grid_corners.__len__()) ) |
---|
[4097] | 376 | v_src_grid_area [:] = src_srf[:].ravel() |
---|
| 377 | v_src_grid_imask [:] = src_msk[:].ravel() |
---|
| 378 | |
---|
| 379 | # -- |
---|
| 380 | |
---|
| 381 | v_dst_grid_dims = f_calving.createVariable ( 'dst_grid_dims' , 'i4', ('dst_grid_rank',) ) |
---|
| 382 | v_dst_grid_center_lon = f_calving.createVariable ( 'dst_grid_center_lon', 'i4', ('dst_grid_size',) ) |
---|
| 383 | v_dst_grid_center_lat = f_calving.createVariable ( 'dst_grid_center_lat', 'i4', ('dst_grid_size',) ) |
---|
[4159] | 384 | v_dst_grid_center_lon.units='degrees_east' ; v_dst_grid_center_lon.long_name='Longitude' |
---|
| 385 | v_dst_grid_center_lon.long_name='longitude' ; v_dst_grid_center_lon.bounds="dst_grid_corner_lon" |
---|
| 386 | v_dst_grid_center_lat.units='degrees_north' ; v_dst_grid_center_lat.long_name='Latitude' |
---|
| 387 | v_dst_grid_center_lat.long_name='latitude' ; v_dst_grid_center_lat.bounds="dst_grid_corner_lat" |
---|
[4097] | 388 | v_dst_grid_corner_lon = f_calving.createVariable ( 'dst_grid_corner_lon', 'f8', ('dst_grid_size', 'dst_grid_corners') ) |
---|
| 389 | v_dst_grid_corner_lat = f_calving.createVariable ( 'dst_grid_corner_lat', 'f8', ('dst_grid_size', 'dst_grid_corners') ) |
---|
| 390 | v_dst_grid_corner_lon.units="degrees_east" |
---|
| 391 | v_dst_grid_corner_lat.units="degrees_north" |
---|
| 392 | v_dst_grid_area = f_calving.createVariable ( 'dst_grid_area' , 'f8', ('dst_grid_size',) ) |
---|
| 393 | v_dst_grid_area.long_name="Grid area" ; v_dst_grid_area.standard_name="cell_area" ; v_dst_grid_area.units="m2" |
---|
| 394 | v_dst_grid_imask = f_calving.createVariable ( 'dst_grid_imask' , 'f8', ('dst_grid_size',) ) |
---|
| 395 | v_dst_grid_imask.long_name="Land-sea mask" ; v_dst_grid_imask.units="Land:1, Ocean:0" |
---|
| 396 | |
---|
[4159] | 397 | v_dst_bassin = f_calving.createVariable ( 'dst_bassin' , 'f8', ('nb_zone', 'dst_grid_size',) ) |
---|
| 398 | v_dst_repartition = f_calving.createVariable ( 'dst_repartition' , 'f8', ('nb_zone', 'dst_grid_size',) ) |
---|
| 399 | |
---|
[4186] | 400 | v_dst_southLimit = f_calving.createVariable ('dst_southLimit', 'f4', ('nb_zone',) ) |
---|
| 401 | v_dst_northLimit = f_calving.createVariable ('dst_northLimit', 'f4', ('nb_zone',) ) |
---|
| 402 | v_dst_southLimit[:] = np.min(limit_lat, axis=(1,) ) |
---|
| 403 | v_dst_northLimit[:] = np.max(limit_lat, axis=(1,) ) |
---|
| 404 | |
---|
[4097] | 405 | v_dst_grid_dims [:] = ( dst_jpi, dst_jpi ) |
---|
| 406 | v_dst_grid_center_lon[:] = dst_lon[:].ravel() |
---|
| 407 | v_dst_grid_center_lat[:] = dst_lat[:].ravel() |
---|
[4159] | 408 | v_dst_grid_corner_lon[:] = dst_clo.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) |
---|
| 409 | v_dst_grid_corner_lat[:] = dst_cla.reshape( (dst_jpi*dst_jpj, d_dst_grid_corners.__len__()) ) |
---|
[4097] | 410 | v_dst_grid_area [:] = dst_srf[:].ravel() |
---|
| 411 | v_dst_grid_imask [:] = dst_msk[:].ravel() |
---|
| 412 | |
---|
[4159] | 413 | v_dst_bassin[:] = basin_msk.reshape ( (nb_zone,dst_grid_size) ) |
---|
| 414 | v_dst_repartition[:] = key_repartition.reshape( (nb_zone,dst_grid_size) ) |
---|
| 415 | |
---|
[4172] | 416 | # Additionnal fields for debugging |
---|
[4159] | 417 | # ================================ |
---|
| 418 | v_dst_grid_imaskutil = f_calving.createVariable ( 'dst_grid_imaskutil' , 'f8', ('dst_grid_size',) ) |
---|
| 419 | v_dst_grid_imaskutil.long_name="Land-sea mask" ; v_dst_grid_imaskutil.units="Land:1, Ocean:0" |
---|
| 420 | v_dst_grid_imaskclose = f_calving.createVariable ( 'dst_grid_imaskclose' , 'f8', ('dst_grid_size',) ) |
---|
| 421 | v_dst_grid_imaskclose.long_name="Land-sea mask" ; v_dst_grid_imaskclose.units="Land:1, Ocean:0" |
---|
| 422 | v_dst_grid_imaskutil [:] = dst_mskutil[:].ravel() |
---|
| 423 | v_dst_grid_imaskclose[:] = dst_closed[:].ravel() |
---|
| 424 | |
---|
[4097] | 425 | v_dst_lon_addressed = f_calving.createVariable ( 'dst_lon_addressed', 'f8', ('num_links',) ) |
---|
| 426 | v_dst_lat_addressed = f_calving.createVariable ( 'dst_lat_addressed', 'f8', ('num_links',) ) |
---|
| 427 | v_dst_lon_addressed.long_name="Longitude" ; v_dst_lon_addressed.standard_name="longitude" ; v_dst_lon_addressed.units="degrees_east" |
---|
[4146] | 428 | v_dst_lat_addressed.long_name="Latitude" ; v_dst_lat_addressed.standard_name="latitude" ; v_dst_lat_addressed.units="degrees_north" |
---|
[4159] | 429 | v_dst_area_addressed = f_calving.createVariable ( 'dst_area_addressed', 'f8', ('num_links',) ) |
---|
| 430 | v_dst_area_addressed.long_name="Grid area" ; v_dst_area_addressed.standard_name="cell_area" ; v_dst_grid_area.units="m2" |
---|
| 431 | v_dst_imask_addressed = f_calving.createVariable ( 'dst_imask_addressed', 'i4', ('num_links',) ) |
---|
| 432 | v_dst_imask_addressed.long_name="Land-sea mask" ; v_dst_imask_addressed.units="Land:1, Ocean:0" |
---|
| 433 | v_dst_imaskutil_addressed = f_calving.createVariable ( 'dst_imaskutil_addressed', 'i4', ('num_links',) ) |
---|
| 434 | v_dst_imaskutil_addressed.long_name="Land-sea mask" ; v_dst_imaskutil_addressed.units="Land:1, Ocean:0" |
---|
[4097] | 435 | |
---|
[4159] | 436 | v_dst_lon_addressed [:] = dst_lon.ravel()[dst_address-1].ravel() |
---|
| 437 | v_dst_lat_addressed [:] = dst_lat.ravel()[dst_address-1].ravel() |
---|
| 438 | v_dst_area_addressed [:] = dst_srf[:].ravel()[dst_address-1].ravel() |
---|
| 439 | v_dst_imask_addressed[:] = dst_msk[:].ravel()[dst_address-1].ravel() |
---|
| 440 | v_dst_imaskutil_addressed[:] = dst_mskutil.ravel()[dst_address-1].ravel() |
---|
| 441 | |
---|
[4172] | 442 | v_src_field = f_calving.createVariable ( 'src_field', 'f8', ('src_grid_size',) ) |
---|
| 443 | v_dst_field = f_calving.createVariable ( 'dst_field', 'f8', ('dst_grid_size',) ) |
---|
| 444 | v_src_field[:] = src_field |
---|
| 445 | v_dst_field[:] = dst_field |
---|
| 446 | |
---|
[4097] | 447 | f_calving.close () |
---|
| 448 | |
---|
| 449 | |
---|
| 450 | ### ===== That's all Folk's !! ============================================== |
---|