[4187] | 1 | # -*- Mode: python -*- |
---|
[6105] | 2 | #!/usr/bin/env python3 |
---|
[4187] | 3 | ### =========================================================================== |
---|
| 4 | ### |
---|
| 5 | ### Compute runoff weights. |
---|
[4199] | 6 | ### For LMDZ only. Not suitable for DYNAMICO |
---|
[4187] | 7 | ### |
---|
| 8 | ### =========================================================================== |
---|
| 9 | ## |
---|
[6190] | 10 | ## MOSAIX is under CeCILL_V2 licence. See "Licence_CeCILL_V2-en.txt" |
---|
| 11 | ## file for an english version of the licence and |
---|
| 12 | ## "Licence_CeCILL_V2-fr.txt" for a french version. |
---|
[4187] | 13 | ## |
---|
[6190] | 14 | ## Permission is hereby granted, free of charge, to any person or |
---|
| 15 | ## organization obtaining a copy of the software and accompanying |
---|
| 16 | ## documentation covered by this license (the "Software") to use, |
---|
| 17 | ## reproduce, display, distribute, execute, and transmit the |
---|
| 18 | ## Software, and to prepare derivative works of the Software, and to |
---|
| 19 | ## permit third-parties to whom the Software is furnished to do so, |
---|
| 20 | ## all subject to the following: |
---|
| 21 | ## |
---|
| 22 | ## Warning, to install, configure, run, use any of MOSAIX software or |
---|
| 23 | ## to read the associated documentation you'll need at least one (1) |
---|
| 24 | ## brain in a reasonably working order. Lack of this implement will |
---|
| 25 | ## void any warranties (either express or implied). Authors assumes |
---|
| 26 | ## no responsability for errors, omissions, data loss, or any other |
---|
| 27 | ## consequences caused directly or indirectly by the usage of his |
---|
| 28 | ## software by incorrectly or partially configured |
---|
| 29 | ## |
---|
| 30 | ## |
---|
[4187] | 31 | # SVN information |
---|
[4188] | 32 | __Author__ = "$Author$" |
---|
| 33 | __Date__ = "$Date$" |
---|
| 34 | __Revision__ = "$Revision$" |
---|
| 35 | __Id__ = "$Id$" |
---|
| 36 | __HeadURL__ = "$HeadURL$" |
---|
[4187] | 37 | __SVN_Date__ = "$SVN_Date: $" |
---|
| 38 | ## |
---|
| 39 | |
---|
| 40 | ## Modules |
---|
[6066] | 41 | import numpy as np, xarray as xr |
---|
[4187] | 42 | import nemo |
---|
| 43 | from scipy import ndimage |
---|
| 44 | import sys, os, platform, argparse, textwrap, time |
---|
| 45 | |
---|
[4259] | 46 | ## Useful constants |
---|
[4187] | 47 | zero = np.dtype('float64').type(0.0) |
---|
| 48 | zone = np.dtype('float64').type(1.0) |
---|
| 49 | epsfrac = np.dtype('float64').type(1.0E-10) |
---|
| 50 | pi = np.pi |
---|
| 51 | rad = pi/np.dtype('float64').type(180.0) # Conversion from degrees to radian |
---|
| 52 | ra = np.dtype('float64').type(6371229.0) # Earth radius |
---|
| 53 | |
---|
| 54 | ## Functions |
---|
| 55 | def geodist (plon1, plat1, plon2, plat2) : |
---|
| 56 | """Distance between two points (on sphere)""" |
---|
| 57 | zs = np.sin (rad*plat1) * np.sin (rad*plat2) + np.cos (rad*plat1) * np.cos (rad*plat2) * np.cos(rad*(plon2-plon1)) |
---|
| 58 | zs = np.maximum (-zone, np.minimum (zone, zs)) |
---|
| 59 | geodist = np.arccos (zs) |
---|
| 60 | return geodist |
---|
| 61 | |
---|
| 62 | ### ===== Reading command line parameters ================================================== |
---|
| 63 | # Creating a parser |
---|
| 64 | parser = argparse.ArgumentParser ( |
---|
| 65 | description = """Compute calving weights""", |
---|
| 66 | epilog='-------- End of the help message --------') |
---|
| 67 | |
---|
| 68 | # Adding arguments |
---|
[6066] | 69 | parser.add_argument ('--oce' , help='oce model name', type=str, default='eORCA1.2', |
---|
[6190] | 70 | choices=['ORCA2.3', 'ORCA1.0', 'ORCA1.1', 'ORCA1_CNRM', 'eORCA1.2', 'eORCA1.4', 'eORCA1.4.2', 'eORCA025', 'eORCA025.1', 'eORCA025.4'] ) |
---|
[4187] | 71 | parser.add_argument ('--atm' , help='atm model name', type=str, default='LMD9695' ) |
---|
| 72 | parser.add_argument ('--atmCoastWidth', help='width of the coastal band in atmosphere (in grid points)', type=int, default=1 ) |
---|
| 73 | parser.add_argument ('--oceCoastWidth', help='width of the coastal band in ocean (in grid points)' , type=int, default=2 ) |
---|
[6066] | 74 | parser.add_argument ('--atmQuantity' , help='Quantity if atm provides quantities (m/s), Surfacic if atm provided flux (m/s/m2)' , type=str, |
---|
| 75 | default='Quantity', choices=['Quantity', 'Surfacic'] ) |
---|
| 76 | parser.add_argument ('--oceQuantity' , help='Quantity if oce requires quantities (ks/s), Surfacic if oce requires flux (m/s/m2)', type=str, |
---|
| 77 | default='Surfacic', choices=['Quantity', 'Surfacic'] ) |
---|
[4187] | 78 | parser.add_argument ('--searchRadius' , help='max distance to connect a land point to an ocean point (in km)', type=float, default=600.0 ) |
---|
| 79 | parser.add_argument ('--grids' , help='grids file name', default='grids.nc' ) |
---|
| 80 | parser.add_argument ('--areas' , help='masks file name', default='areas.nc' ) |
---|
| 81 | parser.add_argument ('--masks' , help='areas file name', default='masks.nc' ) |
---|
[4298] | 82 | parser.add_argument ('--o2a' , help='o2a file name' , default='o2a.nc' ) |
---|
[4199] | 83 | parser.add_argument ('--output', help='output rmp file name', default='rmp_tlmd_to_torc_runoff.nc' ) |
---|
[6066] | 84 | parser.add_argument ('--fmt' , help='NetCDF file format, using nco syntax', default='netcdf4', |
---|
| 85 | choices=['classic', 'netcdf3', '64bit', '64bit_data', '64bit_data', 'netcdf4', 'netcdf4_classsic'] ) |
---|
[6105] | 86 | parser.add_argument ('--ocePerio' , help='periodicity of ocean grid', type=float, default=0, choices=nemo.nperio_valid_range) |
---|
[4187] | 87 | |
---|
| 88 | # Parse command line |
---|
| 89 | myargs = parser.parse_args() |
---|
| 90 | |
---|
| 91 | # |
---|
| 92 | grids = myargs.grids |
---|
| 93 | areas = myargs.areas |
---|
| 94 | masks = myargs.masks |
---|
| 95 | o2a = myargs.o2a |
---|
| 96 | |
---|
| 97 | # Model Names |
---|
| 98 | atm_Name = myargs.atm |
---|
| 99 | oce_Name = myargs.oce |
---|
| 100 | # Width of the coastal band (land points) in the atmopshere |
---|
| 101 | atmCoastWidth = myargs.atmCoastWidth |
---|
| 102 | # Width of the coastal band (ocean points) in the ocean |
---|
| 103 | oceCoastWidth = myargs.oceCoastWidth |
---|
| 104 | searchRadius = myargs.searchRadius * 1000.0 # From km to meters |
---|
| 105 | # Netcdf format |
---|
| 106 | if myargs.fmt == 'classic' : FmtNetcdf = 'CLASSIC' |
---|
| 107 | if myargs.fmt == 'netcdf3' : FmtNetcdf = 'CLASSIC' |
---|
| 108 | if myargs.fmt == '64bit' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' |
---|
| 109 | if myargs.fmt == '64bit_data' : FmtNetcdf = 'NETCDF3_64BIT_DATA' |
---|
| 110 | if myargs.fmt == '64bit_offset' : FmtNetcdf = 'NETCDF3_64BIT_OFFSET' |
---|
| 111 | if myargs.fmt == 'netcdf4' : FmtNetcdf = 'NETCDF4' |
---|
| 112 | if myargs.fmt == 'netcdf4_classic' : FmtNetcdf = 'NETCDF4_CLASSIC' |
---|
| 113 | |
---|
| 114 | # |
---|
| 115 | if atm_Name.find('LMD') >= 0 : atm_n = 'lmd' ; atmDomainType = 'rectilinear' |
---|
| 116 | if atm_Name.find('ICO') >= 0 : atm_n = 'ico' ; atmDomainType = 'unstructured' |
---|
| 117 | |
---|
| 118 | print ('atmQuantity : ' + str (myargs.atmQuantity) ) |
---|
| 119 | print ('oceQuantity : ' + str (myargs.oceQuantity) ) |
---|
[5915] | 120 | |
---|
| 121 | # Ocean grid periodicity |
---|
[6066] | 122 | oce_perio = myargs.ocePerio |
---|
[5915] | 123 | |
---|
[4187] | 124 | ### Read coordinates of all models |
---|
| 125 | ### |
---|
| 126 | |
---|
[6066] | 127 | diaFile = xr.open_dataset ( o2a ) |
---|
| 128 | gridFile = xr.open_dataset ( grids ) |
---|
| 129 | areaFile = xr.open_dataset ( areas ) |
---|
| 130 | maskFile = xr.open_dataset ( masks ) |
---|
[4187] | 131 | |
---|
[6066] | 132 | o2aFrac = diaFile ['OceFrac'].squeeze() |
---|
[4187] | 133 | o2aFrac = np.where ( np.abs(o2aFrac) < 1E10, o2aFrac, 0.0) |
---|
| 134 | |
---|
| 135 | (atm_nvertex, atm_jpj, atm_jpi) = gridFile['t'+atm_n+'.clo'][:].shape |
---|
| 136 | atm_grid_size = atm_jpj*atm_jpi |
---|
| 137 | atm_grid_rank = len(gridFile['t'+atm_n+'.lat'][:].shape) |
---|
| 138 | |
---|
[6066] | 139 | atm_grid_center_lat = gridFile['t'+atm_n+'.lat'].squeeze() |
---|
| 140 | atm_grid_center_lon = gridFile['t'+atm_n+'.lon'].squeeze() |
---|
| 141 | atm_grid_corner_lat = gridFile['t'+atm_n+'.cla'].squeeze() |
---|
| 142 | atm_grid_corner_lon = gridFile['t'+atm_n+'.clo'].squeeze() |
---|
| 143 | atm_grid_area = areaFile['t'+atm_n+'.srf'].squeeze() |
---|
| 144 | atm_grid_imask = 1-maskFile['t'+atm_n+'.msk'][:].squeeze() |
---|
[4187] | 145 | atm_grid_dims = gridFile['t'+atm_n+'.lat'][:].shape |
---|
| 146 | |
---|
[6066] | 147 | if atmDomainType == 'unstructured' : |
---|
| 148 | atm_grid_center_lat = atm_grid_center_lat.rename ({'ycell':'cell'}) |
---|
| 149 | atm_grid_center_lon = atm_grid_center_lon.rename ({'ycell':'cell'}) |
---|
| 150 | atm_grid_corner_lat = atm_grid_corner_lat.rename ({'ycell':'cell'}) |
---|
| 151 | atm_grid_corner_lon = atm_grid_corner_lon.rename ({'ycell':'cell'}) |
---|
| 152 | atm_grid_area = atm_grid_area.rename ({'ycell':'cell'}) |
---|
| 153 | atm_grid_imask = atm_grid_imask.rename ({'ycell':'cell'}) |
---|
| 154 | |
---|
| 155 | if atmDomainType == 'rectilinear' : |
---|
| 156 | atm_grid_center_lat = atm_grid_center_lat.stack (cell=['y', 'x']) |
---|
| 157 | atm_grid_center_lon = atm_grid_center_lon.stack (cell=['y', 'x']) |
---|
| 158 | atm_grid_corner_lat = atm_grid_corner_lat.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'}) |
---|
| 159 | atm_grid_corner_lon = atm_grid_corner_lon.stack (cell=['y', 'x']).rename({'nvertex_lmd':'nvertex'}) |
---|
| 160 | atm_grid_area = atm_grid_area.stack (cell=['y', 'x']) |
---|
| 161 | atm_grid_imask = atm_grid_imask.stack (cell=['y', 'x']) |
---|
| 162 | |
---|
[4187] | 163 | atm_perio = 0 |
---|
| 164 | atm_grid_pmask = atm_grid_imask |
---|
| 165 | atm_address = np.arange(atm_jpj*atm_jpi) |
---|
| 166 | |
---|
| 167 | (oce_nvertex, oce_jpj, oce_jpi) = gridFile['torc.cla'][:].shape ; jpon=oce_jpj*oce_jpj |
---|
| 168 | oce_grid_size = oce_jpj*oce_jpi |
---|
| 169 | oce_grid_rank = len(gridFile['torc.lat'][:].shape) |
---|
| 170 | |
---|
[6066] | 171 | oce_grid_center_lat = gridFile['torc.lat'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 172 | oce_grid_center_lon = gridFile['torc.lon'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 173 | oce_grid_corner_lat = gridFile['torc.cla'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 174 | oce_grid_corner_lon = gridFile['torc.clo'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 175 | oce_grid_area = areaFile['torc.srf'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 176 | oce_grid_imask = 1-maskFile['torc.msk'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
[4187] | 177 | oce_grid_dims = gridFile['torc.lat'][:].shape |
---|
[6066] | 178 | |
---|
[5915] | 179 | if oce_perio == 0 : |
---|
| 180 | if oce_jpi == 182 : oce_perio = 4 # ORCA 2 |
---|
| 181 | if oce_jpi == 362 : oce_perio = 6 # ORCA 1 |
---|
| 182 | if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 |
---|
[6066] | 183 | |
---|
[6105] | 184 | print ("Oce NPERIO parameter : {:}".format(oce_perio)) |
---|
[6066] | 185 | oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel() |
---|
[4187] | 186 | oce_address = np.arange(oce_jpj*oce_jpi) |
---|
| 187 | |
---|
[6042] | 188 | print ("Fill closed sea with image processing library") |
---|
[4187] | 189 | oce_grid_imask2D = np.reshape(oce_grid_pmask,(oce_jpj,oce_jpi)) |
---|
[6066] | 190 | oce_grid_imask2D = nemo.lbc_mask ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), |
---|
| 191 | nperio=oce_perio, cd_type='T', sval=0 ) |
---|
[4187] | 192 | oce_grid_imask = oce_grid_imask2D.ravel() |
---|
| 193 | ## |
---|
[4298] | 194 | print ("Computes an ocean coastal band") |
---|
[4187] | 195 | |
---|
[6066] | 196 | oceLand2D = np.reshape ( np.where (oce_grid_pmask < 0.5, True, False), (oce_jpj, oce_jpi) ) |
---|
| 197 | oceOcean2D = np.reshape ( np.where (oce_grid_pmask > 0.5, True, False), (oce_jpj, oce_jpi) ) |
---|
[4187] | 198 | |
---|
| 199 | NNocean = 1+2*oceCoastWidth |
---|
| 200 | oceOceanFiltered2D = ndimage.uniform_filter(oceOcean2D.astype(float), size=NNocean) |
---|
| 201 | oceCoast2D = np.where (oceOceanFiltered2D<(1.0-0.5/(NNocean**2)),True,False) & oceOcean2D |
---|
[6064] | 202 | oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D,(oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T').ravel() |
---|
[4187] | 203 | |
---|
| 204 | oceOceanFiltered = oceOceanFiltered2D.ravel() |
---|
[6066] | 205 | oceLand = oceLand2D.ravel () |
---|
[4187] | 206 | oceOcean = oceOcean2D.ravel() |
---|
| 207 | oceCoast = oceCoast2D.ravel() |
---|
| 208 | |
---|
[6042] | 209 | print ('Number of points in oceLand : {:8d}'.format (oceLand.sum()) ) |
---|
| 210 | print ('Number of points in oceOcean : {:8d}'.format (oceOcean.sum()) ) |
---|
| 211 | print ('Number of points in oceCoast : {:8d}'.format (oceCoast.sum()) ) |
---|
[4187] | 212 | |
---|
| 213 | # Arrays with coastal points only |
---|
| 214 | oceCoast_grid_center_lon = oce_grid_center_lon[oceCoast] |
---|
| 215 | oceCoast_grid_center_lat = oce_grid_center_lat[oceCoast] |
---|
| 216 | oceCoast_grid_area = oce_grid_area [oceCoast] |
---|
| 217 | oceCoast_grid_imask = oce_grid_imask [oceCoast] |
---|
| 218 | oceCoast_grid_pmask = oce_grid_pmask [oceCoast] |
---|
| 219 | oceCoast_address = oce_address [oceCoast] |
---|
| 220 | |
---|
[4298] | 221 | print ("Computes an atmosphere coastal band " ) |
---|
[4187] | 222 | atmLand = np.where (o2aFrac[:] < epsfrac , True, False) |
---|
| 223 | atmLandFrac = np.where (o2aFrac[:] < zone-epsfrac , True, False) |
---|
| 224 | atmFrac = np.where (o2aFrac[:] > epsfrac , True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) |
---|
| 225 | atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False) |
---|
| 226 | atmOceanFrac = np.where (o2aFrac[:] > epsfrac , True, False) |
---|
| 227 | |
---|
[4298] | 228 | ## For LMDZ only !! |
---|
[4187] | 229 | if atmDomainType == 'rectilinear' : |
---|
[6042] | 230 | print ("Extend coastal band " ) |
---|
[4187] | 231 | NNatm = 1+2*atmCoastWidth |
---|
| 232 | atmLand2D = np.reshape ( atmLand, ( atm_jpj, atm_jpi) ) |
---|
| 233 | |
---|
| 234 | atmLandFiltered2D = ndimage.uniform_filter(atmLand2D.astype(float), size=NNatm) |
---|
| 235 | atmCoast2D = np.where (atmLandFiltered2D<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac |
---|
| 236 | |
---|
| 237 | atmLandFiltered = atmLandFiltered2D.ravel() |
---|
| 238 | atmCoast = atmCoast2D.ravel() |
---|
| 239 | |
---|
[6042] | 240 | print ('Number of points in atmLand : {:8d}'.format (atmLand.sum()) ) |
---|
| 241 | print ('Number of points in atmOcean : {:8d}'.format (atmOcean.sum()) ) |
---|
| 242 | print ('Number of points in atmCoast : {:8d}'.format (atmCoast.sum()) ) |
---|
[4187] | 243 | |
---|
| 244 | else : |
---|
| 245 | atmCoast = atmFrac |
---|
| 246 | |
---|
| 247 | |
---|
| 248 | # Arrays with coastal points only |
---|
| 249 | atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast] |
---|
| 250 | atmCoast_grid_center_lat = atm_grid_center_lat[atmCoast] |
---|
| 251 | atmCoast_grid_area = atm_grid_area [atmCoast] |
---|
| 252 | atmCoast_grid_imask = atm_grid_imask [atmCoast] |
---|
| 253 | atmCoast_grid_pmask = atm_grid_pmask [atmCoast] |
---|
| 254 | atmCoast_address = atm_address [atmCoast] |
---|
| 255 | |
---|
| 256 | # Initialisations before the loop |
---|
| 257 | remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) |
---|
| 258 | atm_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 259 | oce_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 260 | |
---|
| 261 | ## Loop on atmosphere coastal points |
---|
[6042] | 262 | if searchRadius > 0. : |
---|
| 263 | print ("Loop on atmosphere coastal points") |
---|
| 264 | for ja in np.arange(len(atmCoast_grid_pmask)) : |
---|
| 265 | z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) |
---|
| 266 | z_mask = np.where (z_dist*ra < searchRadius, True, False) |
---|
| 267 | num_links = int(z_mask.sum()) |
---|
| 268 | if num_links == 0 : continue |
---|
[6066] | 269 | z_area = oceCoast_grid_area[z_mask].sum().values |
---|
[6042] | 270 | poids = np.ones ((num_links),dtype=np.float64) / z_area |
---|
| 271 | if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja] |
---|
| 272 | if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] |
---|
| 273 | if ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print |
---|
[6066] | 274 | print ( 'ja:{:8d}, num_links:{:8d}, z_area:{:8.4e}, atm area:{:8.4e}, weights sum:{:8.4e} ' |
---|
| 275 | .format(ja, num_links, z_area, atm_grid_area[ja].values, poids.sum() ) ) |
---|
[6042] | 276 | # |
---|
| 277 | matrix_local = poids |
---|
| 278 | atm_address_local = np.ones(num_links, dtype=np.int32 ) * atmCoast_address[ja] |
---|
| 279 | # Address on destination grid |
---|
| 280 | oce_address_local = oceCoast_address[z_mask] |
---|
| 281 | # Append to global arrays |
---|
| 282 | remap_matrix = np.append ( remap_matrix, matrix_local ) |
---|
| 283 | atm_address = np.append ( atm_address , atm_address_local ) |
---|
| 284 | oce_address = np.append ( oce_address , oce_address_local ) |
---|
[4187] | 285 | |
---|
[6042] | 286 | print ('End of loop') |
---|
[4187] | 287 | |
---|
| 288 | num_links = remap_matrix.shape[0] |
---|
| 289 | |
---|
[6042] | 290 | print ("Write output file") |
---|
[4187] | 291 | runoff = myargs.output |
---|
| 292 | print ('Output file: ' + runoff ) |
---|
| 293 | |
---|
| 294 | |
---|
[6066] | 295 | remap_matrix = xr.DataArray ( np.reshape(remap_matrix, (num_links, 1)), dims = ['num_links', 'num_wgts'] ) |
---|
[4187] | 296 | |
---|
[6066] | 297 | # OASIS uses Fortran style indexing, starting at one |
---|
| 298 | src_address = xr.DataArray ( atm_address.astype(np.int32)+1, dims = ['num_links'], |
---|
| 299 | attrs={"convention": "Fortran style addressing, starting at 1"}) |
---|
| 300 | dst_address = xr.DataArray ( oce_address.astype(np.int32)+1, dims = ['num_links'], |
---|
| 301 | attrs={"convention": "Fortran style addressing, starting at 1"}) |
---|
[4187] | 302 | |
---|
[6066] | 303 | src_grid_dims = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32), dims = ['src_grid_rank',] ) |
---|
| 304 | src_grid_center_lon = xr.DataArray (atm_grid_center_lon.values , dims = ['src_grid_size',] ) |
---|
| 305 | src_grid_center_lat = xr.DataArray (atm_grid_center_lat.values , dims = ['src_grid_size',] ) |
---|
[4187] | 306 | |
---|
[6066] | 307 | src_grid_center_lon.attrs['units']='degrees_east' ; src_grid_center_lon.attrs['long_name']='Longitude' |
---|
| 308 | src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" |
---|
| 309 | src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' |
---|
| 310 | src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" |
---|
| 311 | |
---|
| 312 | src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) |
---|
| 313 | src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) |
---|
| 314 | src_grid_corner_lon.attrs['units']="degrees_east" |
---|
| 315 | src_grid_corner_lat.attrs['units']="degrees_north" |
---|
[4187] | 316 | |
---|
[6066] | 317 | src_grid_area = xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] ) |
---|
| 318 | src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" |
---|
[4187] | 319 | |
---|
[6066] | 320 | src_grid_imask = xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] ) |
---|
| 321 | src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 322 | |
---|
[6066] | 323 | src_grid_pmask = xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] ) |
---|
| 324 | src_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; src_grid_pmask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 325 | |
---|
| 326 | # -- |
---|
[6066] | 327 | dst_grid_dims = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32), dims = ['dst_grid_rank',] ) |
---|
| 328 | dst_grid_center_lon = xr.DataArray (oce_grid_center_lon.values, dims = ['dst_grid_size',] ) |
---|
| 329 | dst_grid_center_lat = xr.DataArray (oce_grid_center_lat.values, dims = ['dst_grid_size',] ) |
---|
[4187] | 330 | |
---|
[6066] | 331 | dst_grid_center_lon.attrs['units']='degrees_east' ; dst_grid_center_lon.attrs['long_name']='Longitude' |
---|
| 332 | dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" |
---|
| 333 | dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' |
---|
| 334 | dst_grid_center_lat.attrs['long_name']='latitude ' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" |
---|
[4187] | 335 | |
---|
[6066] | 336 | dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) |
---|
| 337 | dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) |
---|
| 338 | dst_grid_corner_lon.attrs['units']="degrees_east" |
---|
| 339 | dst_grid_corner_lat.attrs['units']="degrees_north" |
---|
[4187] | 340 | |
---|
[6066] | 341 | dst_grid_area = xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] ) |
---|
| 342 | dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" |
---|
[4187] | 343 | |
---|
[6066] | 344 | dst_grid_imask = xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] ) |
---|
| 345 | dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 346 | |
---|
[6066] | 347 | dst_grid_pmask = xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] ) |
---|
| 348 | dst_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; dst_grid_pmask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 349 | |
---|
[6066] | 350 | src_lon_addressed = xr.DataArray (atm_grid_center_lon.values[atm_address] , dims = ['num_links'] ) |
---|
| 351 | src_lat_addressed = xr.DataArray (atm_grid_center_lat.values[atm_address] , dims = ['num_links'] ) |
---|
| 352 | src_area_addressed = xr.DataArray (atm_grid_area .values[atm_address] , dims = ['num_links'] ) |
---|
| 353 | src_imask_addressed = xr.DataArray (1-atm_grid_imask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 354 | src_pmask_addressed = xr.DataArray (1-atm_grid_pmask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) |
---|
[4187] | 355 | |
---|
[6066] | 356 | dst_lon_addressed = xr.DataArray (oce_grid_center_lon.values[atm_address], dims = ['num_links'] ) |
---|
| 357 | dst_lat_addressed = xr.DataArray (oce_grid_center_lat.values[oce_address], dims = ['num_links'] ) |
---|
| 358 | dst_area_addressed = xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 359 | dst_imask_addressed = xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 360 | dst_pmask_addressed = xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 361 | |
---|
| 362 | src_lon_addressed.attrs['long_name']="Longitude" ; src_lon_addressed.attrs['standard_name']="longitude" ; src_lon_addressed.attrs['units']="degrees_east" |
---|
| 363 | src_lat_addressed.attrs['long_name']="Latitude" ; src_lat_addressed.attrs['standard_name']="latitude" ; src_lat_addressed.attrs['units']="degrees_north" |
---|
| 364 | |
---|
| 365 | dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" |
---|
| 366 | dst_lat_addressed.attrs['long_name']="Latitude" ; dst_lat_addressed.attrs['standard_name']="latitude" ; dst_lat_addressed.attrs['units']="degrees_north" |
---|
| 367 | |
---|
[4187] | 368 | if atmDomainType == 'rectilinear' : |
---|
[6066] | 369 | atmLand = xr.DataArray ( atmLand.ravel() , dims = ['src_grid_size',] ) |
---|
| 370 | atmLandFiltered = xr.DataArray ( atmLandFrac.ravel() , dims = ['src_grid_size',] ) |
---|
| 371 | atmLandFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) |
---|
| 372 | atmFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) |
---|
| 373 | atmOcean = xr.DataArray ( atmOcean.ravel() , dims = ['src_grid_size',] ) |
---|
| 374 | atmOceanFrac = xr.DataArray ( atmOceanFrac.ravel(), dims = ['src_grid_size',] ) |
---|
[4187] | 375 | |
---|
[6066] | 376 | atmCoast = xr.DataArray (atmCoast.astype(np.int32) , dims = ['src_grid_size',]) |
---|
| 377 | oceLand = xr.DataArray (oceLand.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
| 378 | oceOcean = xr.DataArray (oceOcean.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
| 379 | oceOceanFiltered = xr.DataArray (oceOceanFiltered.astype(np.float32), dims = ['dst_grid_size',]) |
---|
| 380 | oceCoast = xr.DataArray (oceCoast.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
[4187] | 381 | |
---|
| 382 | |
---|
[6066] | 383 | f_runoff = xr.Dataset ( { |
---|
[6090] | 384 | 'remap_matrix' : remap_matrix, |
---|
[6066] | 385 | 'src_address' : src_address, |
---|
| 386 | 'dst_address' : dst_address, |
---|
| 387 | 'src_grid_dims' : src_grid_dims, |
---|
| 388 | 'src_grid_center_lon' : src_grid_center_lon, |
---|
| 389 | 'src_grid_center_lat' : src_grid_center_lat, |
---|
| 390 | 'src_grid_corner_lon' : src_grid_corner_lon, |
---|
| 391 | 'src_grid_corner_lat' : src_grid_corner_lat, |
---|
| 392 | 'src_grid_area' : src_grid_area, |
---|
| 393 | 'src_grid_area' : src_grid_area, |
---|
| 394 | 'src_grid_pmask' : src_grid_pmask, |
---|
| 395 | 'dst_grid_dims' : dst_grid_dims, |
---|
| 396 | 'dst_grid_center_lon' : dst_grid_center_lon, |
---|
| 397 | 'st_grid_center_lat' : dst_grid_center_lat, |
---|
| 398 | 'dst_grid_corner_lon' : dst_grid_corner_lon, |
---|
| 399 | 'dst_grid_corner_lat' : dst_grid_corner_lat, |
---|
| 400 | 'dst_grid_area' : dst_grid_area, |
---|
| 401 | 'dst_grid_imask' : dst_grid_imask, |
---|
| 402 | 'dst_grid_pmask' : dst_grid_pmask, |
---|
| 403 | 'src_lon_addressed' : src_lon_addressed, |
---|
| 404 | 'src_lat_addressed' : src_lat_addressed, |
---|
| 405 | 'src_area_addressed' : src_area_addressed, |
---|
| 406 | 'dst_lon_addressed' : dst_lon_addressed, |
---|
| 407 | 'dst_lat_addressed' : dst_lat_addressed, |
---|
| 408 | 'dst_area_addressed' : dst_area_addressed, |
---|
| 409 | 'dst_imask_addressed' : dst_imask_addressed, |
---|
| 410 | 'dst_pmask_addressed' : dst_pmask_addressed, |
---|
| 411 | 'atmCoast' : atmCoast, |
---|
| 412 | 'oceLand' : oceLand, |
---|
| 413 | 'oceOcean' : oceOcean, |
---|
| 414 | 'oceOceanFiltered' : oceOceanFiltered, |
---|
| 415 | 'oceCoast' : oceCoast |
---|
| 416 | } ) |
---|
[4187] | 417 | |
---|
[6066] | 418 | f_runoff.attrs['Conventions'] = "CF-1.6" |
---|
| 419 | f_runoff.attrs['source'] = "IPSL Earth system model" |
---|
| 420 | f_runoff.attrs['group'] = "ICMC IPSL Climate Modelling Center" |
---|
| 421 | f_runoff.attrs['Institution'] = "IPSL https.//www.ipsl.fr" |
---|
| 422 | f_runoff.attrs['Ocean'] = oce_Name + " https://www.nemo-ocean.eu" |
---|
| 423 | f_runoff.attrs['Atmosphere'] = atm_Name + " http://lmdz.lmd.jussieu.fr" |
---|
| 424 | f_runoff.attrs['associatedFiles'] = grids + " " + areas + " " + masks |
---|
| 425 | f_runoff.attrs['description'] = "Generated with RunoffWeights.py" |
---|
| 426 | f_runoff.attrs['title'] = runoff |
---|
[6090] | 427 | f_runoff.attrs['Program'] = "Generated by " + sys.argv[0] + " with flags " + ' '.join (sys.argv[1:]) |
---|
[6066] | 428 | f_runoff.attrs['atmCoastWidth'] = "{:d} grid points".format(atmCoastWidth) |
---|
| 429 | f_runoff.attrs['oceCoastWidth'] = "{:d} grid points".format(oceCoastWidth) |
---|
| 430 | f_runoff.attrs['searchRadius'] = "{:.0f} km".format(searchRadius/1000.) |
---|
| 431 | f_runoff.attrs['atmQuantity'] = myargs.atmQuantity |
---|
| 432 | f_runoff.attrs['oceQuantity'] = myargs.oceQuantity |
---|
| 433 | f_runoff.attrs['gridsFile'] = grids |
---|
| 434 | f_runoff.attrs['areasFile'] = areas |
---|
| 435 | f_runoff.attrs['masksFile'] = masks |
---|
| 436 | f_runoff.attrs['o2aFile'] = o2a |
---|
[6090] | 437 | f_runoff.attrs['timeStamp'] = time.asctime () |
---|
| 438 | try : f_calving.attrs['directory'] = os.getcwd () |
---|
| 439 | except : pass |
---|
| 440 | try : f_runoff.attrs['HOSTNAME'] = platform.node () |
---|
| 441 | except : pass |
---|
| 442 | try : f_runoff.attrs['LOGNAME'] = os.getlogin () |
---|
| 443 | except : pass |
---|
| 444 | try : f_runoff.attrs['Python'] = "Python version " + platform.python_version () |
---|
| 445 | except : pass |
---|
| 446 | try : f_runoff.attrs['OS'] = platform.system () |
---|
| 447 | except : pass |
---|
| 448 | try : f_runoff.attrs['release'] = platform.release () |
---|
| 449 | except : pass |
---|
| 450 | try : f_runoff.attrs['hardware'] = platform.machine () |
---|
| 451 | except : pass |
---|
[6066] | 452 | f_runoff.attrs['conventions'] = "SCRIP" |
---|
| 453 | f_runoff.attrs['source_grid'] = "curvilinear" |
---|
| 454 | f_runoff.attrs['dest_grid'] = "curvilinear" |
---|
| 455 | f_runoff.attrs['Model'] = "IPSL CM6" |
---|
| 456 | f_runoff.attrs['SVN_Author'] = "$Author$" |
---|
| 457 | f_runoff.attrs['SVN_Date'] = "$Date$" |
---|
| 458 | f_runoff.attrs['SVN_Revision'] = "$Revision$" |
---|
| 459 | f_runoff.attrs['SVN_Id'] = "$Id$" |
---|
| 460 | f_runoff.attrs['SVN_HeadURL'] = "$HeadURL$" |
---|
| 461 | |
---|
| 462 | f_runoff.to_netcdf ( runoff, mode='w', format=FmtNetcdf ) |
---|
[6042] | 463 | f_runoff.close () |
---|
[4187] | 464 | |
---|
[6066] | 465 | ## |
---|
| 466 | |
---|
[6042] | 467 | print ('That''s all folks !') |
---|
[4187] | 468 | ## ====================================================================================== |
---|