[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 |
---|
[6244] | 165 | |
---|
| 166 | |
---|
[4187] | 167 | atm_address = np.arange(atm_jpj*atm_jpi) |
---|
| 168 | |
---|
| 169 | (oce_nvertex, oce_jpj, oce_jpi) = gridFile['torc.cla'][:].shape ; jpon=oce_jpj*oce_jpj |
---|
| 170 | oce_grid_size = oce_jpj*oce_jpi |
---|
| 171 | oce_grid_rank = len(gridFile['torc.lat'][:].shape) |
---|
| 172 | |
---|
[6066] | 173 | oce_grid_center_lat = gridFile['torc.lat'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 174 | oce_grid_center_lon = gridFile['torc.lon'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 175 | oce_grid_corner_lat = gridFile['torc.cla'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 176 | oce_grid_corner_lon = gridFile['torc.clo'].squeeze().stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 177 | oce_grid_area = areaFile['torc.srf'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
| 178 | oce_grid_imask = 1-maskFile['torc.msk'].stack(oce_grid_size=['y_grid_T', 'x_grid_T']) |
---|
[4187] | 179 | oce_grid_dims = gridFile['torc.lat'][:].shape |
---|
[6066] | 180 | |
---|
[6244] | 181 | # imask : includes periodicity point |
---|
| 182 | # pmask : exclude periodicity points |
---|
| 183 | |
---|
| 184 | #if oce_perio == 0 : |
---|
| 185 | # if oce_jpi == 182 : oce_perio = 4 # ORCA 2 |
---|
| 186 | # if oce_jpi == 362 : oce_perio = 6 # ORCA 1 |
---|
| 187 | # if oce_jpi == 1442 : oce_perio = 6 # ORCA 025 |
---|
| 188 | oce_preio = nemo.__guessNperio__ (oce_jpj, oce_jpi, nperio=oce_perio) |
---|
| 189 | |
---|
[6105] | 190 | print ("Oce NPERIO parameter : {:}".format(oce_perio)) |
---|
[6244] | 191 | oce_grid_imask = nemo.lbc (np.reshape(oce_grid_imask.values, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T' ).ravel() |
---|
| 192 | oce_grid_pmask = nemo.lbc_mask (np.reshape(oce_grid_imask , (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0).ravel() |
---|
[4187] | 193 | oce_address = np.arange(oce_jpj*oce_jpi) |
---|
| 194 | |
---|
[6244] | 195 | print ("Fill closed seas with image processing library") |
---|
| 196 | oce_grid_imask2D = np.reshape (oce_grid_imask, (oce_jpj,oce_jpi)) |
---|
| 197 | oce_grid_imask2D = 1-nemo.fill_closed_seas (1-oce_grid_imask2D, nperio=oce_perio, cd_type='T') |
---|
| 198 | #oce_grid_imask2D = nemo.lbc ( 1-ndimage.binary_fill_holes (1-nemo.lbc(oce_grid_imask2D, nperio=oce_perio, cd_type='T')), |
---|
| 199 | # nperio=oce_perio, cd_type='T' ) |
---|
| 200 | oce_grid_imask = oce_grid_imask2D.ravel () |
---|
| 201 | |
---|
[4187] | 202 | ## |
---|
[6244] | 203 | print ("Computes an ocean coastal band with image processing library") |
---|
| 204 | oceLand2D = np.reshape ( np.where (oce_grid_imask < 0.5, True, False), (oce_jpj, oce_jpi) ) |
---|
| 205 | oceOcean2D = np.reshape ( np.where (oce_grid_imask > 0.5, True, False), (oce_jpj, oce_jpi) ) |
---|
[4187] | 206 | |
---|
| 207 | NNocean = 1+2*oceCoastWidth |
---|
[6244] | 208 | oceOceanFiltered2D = ndimage.uniform_filter (oceOcean2D.astype(float), size=NNocean) |
---|
[4187] | 209 | oceCoast2D = np.where (oceOceanFiltered2D<(1.0-0.5/(NNocean**2)),True,False) & oceOcean2D |
---|
[6244] | 210 | oceCoast2D = nemo.lbc_mask (np.reshape(oceCoast2D, (oce_jpj,oce_jpi)), nperio=oce_perio, cd_type='T', sval=0.).ravel() |
---|
[4187] | 211 | |
---|
| 212 | oceOceanFiltered = oceOceanFiltered2D.ravel() |
---|
[6066] | 213 | oceLand = oceLand2D.ravel () |
---|
[4187] | 214 | oceOcean = oceOcean2D.ravel() |
---|
| 215 | oceCoast = oceCoast2D.ravel() |
---|
| 216 | |
---|
[6042] | 217 | print ('Number of points in oceLand : {:8d}'.format (oceLand.sum()) ) |
---|
| 218 | print ('Number of points in oceOcean : {:8d}'.format (oceOcean.sum()) ) |
---|
| 219 | print ('Number of points in oceCoast : {:8d}'.format (oceCoast.sum()) ) |
---|
[4187] | 220 | |
---|
| 221 | # Arrays with coastal points only |
---|
| 222 | oceCoast_grid_center_lon = oce_grid_center_lon[oceCoast] |
---|
| 223 | oceCoast_grid_center_lat = oce_grid_center_lat[oceCoast] |
---|
| 224 | oceCoast_grid_area = oce_grid_area [oceCoast] |
---|
| 225 | oceCoast_grid_imask = oce_grid_imask [oceCoast] |
---|
| 226 | oceCoast_grid_pmask = oce_grid_pmask [oceCoast] |
---|
| 227 | oceCoast_address = oce_address [oceCoast] |
---|
| 228 | |
---|
[6244] | 229 | print ("Computes an atmosphere coastal band" ) |
---|
[4187] | 230 | atmLand = np.where (o2aFrac[:] < epsfrac , True, False) |
---|
| 231 | atmLandFrac = np.where (o2aFrac[:] < zone-epsfrac , True, False) |
---|
| 232 | atmFrac = np.where (o2aFrac[:] > epsfrac , True, False) & np.where (o2aFrac[:] < (zone-epsfrac), True, False) |
---|
| 233 | atmOcean = np.where (o2aFrac[:] < (zone-epsfrac), True, False) |
---|
| 234 | atmOceanFrac = np.where (o2aFrac[:] > epsfrac , True, False) |
---|
| 235 | |
---|
[4298] | 236 | ## For LMDZ only !! |
---|
[4187] | 237 | if atmDomainType == 'rectilinear' : |
---|
[6042] | 238 | print ("Extend coastal band " ) |
---|
[4187] | 239 | NNatm = 1+2*atmCoastWidth |
---|
| 240 | atmLand2D = np.reshape ( atmLand, ( atm_jpj, atm_jpi) ) |
---|
| 241 | |
---|
| 242 | atmLandFiltered2D = ndimage.uniform_filter(atmLand2D.astype(float), size=NNatm) |
---|
| 243 | atmCoast2D = np.where (atmLandFiltered2D<(1.0-0.5/(NNatm**2)),True,False) & atmLandFrac |
---|
| 244 | |
---|
| 245 | atmLandFiltered = atmLandFiltered2D.ravel() |
---|
| 246 | atmCoast = atmCoast2D.ravel() |
---|
| 247 | |
---|
[6042] | 248 | print ('Number of points in atmLand : {:8d}'.format (atmLand.sum()) ) |
---|
| 249 | print ('Number of points in atmOcean : {:8d}'.format (atmOcean.sum()) ) |
---|
| 250 | print ('Number of points in atmCoast : {:8d}'.format (atmCoast.sum()) ) |
---|
[4187] | 251 | |
---|
| 252 | else : |
---|
| 253 | atmCoast = atmFrac |
---|
| 254 | |
---|
| 255 | |
---|
| 256 | # Arrays with coastal points only |
---|
| 257 | atmCoast_grid_center_lon = atm_grid_center_lon[atmCoast] |
---|
| 258 | atmCoast_grid_center_lat = atm_grid_center_lat[atmCoast] |
---|
| 259 | atmCoast_grid_area = atm_grid_area [atmCoast] |
---|
| 260 | atmCoast_grid_imask = atm_grid_imask [atmCoast] |
---|
| 261 | atmCoast_grid_pmask = atm_grid_pmask [atmCoast] |
---|
| 262 | atmCoast_address = atm_address [atmCoast] |
---|
| 263 | |
---|
| 264 | # Initialisations before the loop |
---|
| 265 | remap_matrix = np.empty ( shape=(0), dtype=np.float64 ) |
---|
| 266 | atm_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 267 | oce_address = np.empty ( shape=(0), dtype=np.int32 ) |
---|
| 268 | |
---|
| 269 | ## Loop on atmosphere coastal points |
---|
[6042] | 270 | if searchRadius > 0. : |
---|
| 271 | print ("Loop on atmosphere coastal points") |
---|
[6244] | 272 | for ja in np.arange (len(atmCoast_grid_pmask)) : |
---|
[6042] | 273 | z_dist = geodist ( atmCoast_grid_center_lon[ja], atmCoast_grid_center_lat[ja], oceCoast_grid_center_lon, oceCoast_grid_center_lat) |
---|
| 274 | z_mask = np.where (z_dist*ra < searchRadius, True, False) |
---|
[6244] | 275 | num_links = int (z_mask.sum()) |
---|
[6042] | 276 | if num_links == 0 : continue |
---|
[6066] | 277 | z_area = oceCoast_grid_area[z_mask].sum().values |
---|
[6042] | 278 | poids = np.ones ((num_links),dtype=np.float64) / z_area |
---|
| 279 | if myargs.atmQuantity == 'Surfacic' : poids = poids * atm_grid_area[ja] |
---|
| 280 | if myargs.oceQuantity == 'Quantity' : poids = poids * oceCoast_grid_area[z_mask] |
---|
| 281 | if ja % (len(atmCoast_grid_pmask)//50) == 0 : # Control print |
---|
[6066] | 282 | print ( 'ja:{:8d}, num_links:{:8d}, z_area:{:8.4e}, atm area:{:8.4e}, weights sum:{:8.4e} ' |
---|
| 283 | .format(ja, num_links, z_area, atm_grid_area[ja].values, poids.sum() ) ) |
---|
[6042] | 284 | # |
---|
| 285 | matrix_local = poids |
---|
[6244] | 286 | atm_address_local = np.ones (num_links, dtype=np.int32 ) * atmCoast_address[ja] |
---|
[6042] | 287 | # Address on destination grid |
---|
[6244] | 288 | oce_address_local = oceCoast_address [z_mask] |
---|
[6042] | 289 | # Append to global arrays |
---|
| 290 | remap_matrix = np.append ( remap_matrix, matrix_local ) |
---|
| 291 | atm_address = np.append ( atm_address , atm_address_local ) |
---|
| 292 | oce_address = np.append ( oce_address , oce_address_local ) |
---|
[4187] | 293 | |
---|
[6042] | 294 | print ('End of loop') |
---|
[4187] | 295 | |
---|
| 296 | num_links = remap_matrix.shape[0] |
---|
| 297 | |
---|
[6042] | 298 | print ("Write output file") |
---|
[4187] | 299 | runoff = myargs.output |
---|
| 300 | print ('Output file: ' + runoff ) |
---|
| 301 | |
---|
| 302 | |
---|
[6066] | 303 | remap_matrix = xr.DataArray ( np.reshape(remap_matrix, (num_links, 1)), dims = ['num_links', 'num_wgts'] ) |
---|
[4187] | 304 | |
---|
[6066] | 305 | # OASIS uses Fortran style indexing, starting at one |
---|
| 306 | src_address = xr.DataArray ( atm_address.astype(np.int32)+1, dims = ['num_links'], |
---|
| 307 | attrs={"convention": "Fortran style addressing, starting at 1"}) |
---|
| 308 | dst_address = xr.DataArray ( oce_address.astype(np.int32)+1, dims = ['num_links'], |
---|
| 309 | attrs={"convention": "Fortran style addressing, starting at 1"}) |
---|
[4187] | 310 | |
---|
[6066] | 311 | src_grid_dims = xr.DataArray (np.array(atm_grid_dims, dtype=np.int32), dims = ['src_grid_rank',] ) |
---|
| 312 | src_grid_center_lon = xr.DataArray (atm_grid_center_lon.values , dims = ['src_grid_size',] ) |
---|
| 313 | src_grid_center_lat = xr.DataArray (atm_grid_center_lat.values , dims = ['src_grid_size',] ) |
---|
[4187] | 314 | |
---|
[6066] | 315 | src_grid_center_lon.attrs['units']='degrees_east' ; src_grid_center_lon.attrs['long_name']='Longitude' |
---|
| 316 | src_grid_center_lon.attrs['long_name']='longitude' ; src_grid_center_lon.attrs['bounds']="src_grid_corner_lon" |
---|
| 317 | src_grid_center_lat.attrs['units']='degrees_north' ; src_grid_center_lat.attrs['long_name']='Latitude' |
---|
| 318 | src_grid_center_lat.attrs['long_name']='latitude ' ; src_grid_center_lat.attrs['bounds']="src_grid_corner_lat" |
---|
| 319 | |
---|
| 320 | src_grid_corner_lon = xr.DataArray (atm_grid_corner_lon.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) |
---|
| 321 | src_grid_corner_lat = xr.DataArray (atm_grid_corner_lat.values.transpose(), dims = [ 'src_grid_size', 'src_grid_corners'] ) |
---|
| 322 | src_grid_corner_lon.attrs['units']="degrees_east" |
---|
| 323 | src_grid_corner_lat.attrs['units']="degrees_north" |
---|
[4187] | 324 | |
---|
[6066] | 325 | src_grid_area = xr.DataArray (atm_grid_area.values, dims = ['src_grid_size',] ) |
---|
| 326 | src_grid_area.attrs['long_name']="Grid area" ; src_grid_area.attrs['standard_name']="cell_area" ; src_grid_area.attrs['units']="m2" |
---|
[4187] | 327 | |
---|
[6066] | 328 | src_grid_imask = xr.DataArray (atm_grid_imask.values, dims = ['src_grid_size',] ) |
---|
| 329 | src_grid_imask.attrs['long_name']="Land-sea mask" ; src_grid_imask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 330 | |
---|
[6066] | 331 | src_grid_pmask = xr.DataArray (atm_grid_pmask.values, dims = ['src_grid_size',] ) |
---|
| 332 | src_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; src_grid_pmask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 333 | |
---|
| 334 | # -- |
---|
[6066] | 335 | dst_grid_dims = xr.DataArray (np.array(oce_grid_dims, dtype=np.int32), dims = ['dst_grid_rank',] ) |
---|
| 336 | dst_grid_center_lon = xr.DataArray (oce_grid_center_lon.values, dims = ['dst_grid_size',] ) |
---|
| 337 | dst_grid_center_lat = xr.DataArray (oce_grid_center_lat.values, dims = ['dst_grid_size',] ) |
---|
[4187] | 338 | |
---|
[6066] | 339 | dst_grid_center_lon.attrs['units']='degrees_east' ; dst_grid_center_lon.attrs['long_name']='Longitude' |
---|
| 340 | dst_grid_center_lon.attrs['long_name']='longitude' ; dst_grid_center_lon.attrs['bounds']="dst_grid_corner_lon" |
---|
| 341 | dst_grid_center_lat.attrs['units']='degrees_north' ; dst_grid_center_lat.attrs['long_name']='Latitude' |
---|
| 342 | dst_grid_center_lat.attrs['long_name']='latitude ' ; dst_grid_center_lat.attrs['bounds']="dst_grid_corner_lat" |
---|
[4187] | 343 | |
---|
[6066] | 344 | dst_grid_corner_lon = xr.DataArray (np.transpose(oce_grid_corner_lon.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) |
---|
| 345 | dst_grid_corner_lat = xr.DataArray (np.transpose(oce_grid_corner_lat.values), dims = [ 'dst_grid_size', 'dst_grid_corners'] ) |
---|
| 346 | dst_grid_corner_lon.attrs['units']="degrees_east" |
---|
| 347 | dst_grid_corner_lat.attrs['units']="degrees_north" |
---|
[4187] | 348 | |
---|
[6066] | 349 | dst_grid_area = xr.DataArray (oce_grid_area.values, dims = ['dst_grid_size',] ) |
---|
| 350 | dst_grid_area.attrs['long_name']="Grid area" ; dst_grid_area.attrs['standard_name']="cell_area" ; dst_grid_area.attrs['units']="m2" |
---|
[4187] | 351 | |
---|
[6066] | 352 | dst_grid_imask = xr.DataArray (oce_grid_imask.astype(np.int32), dims = ['dst_grid_size',] ) |
---|
| 353 | dst_grid_imask.attrs['long_name']="Land-sea mask" ; dst_grid_imask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 354 | |
---|
[6066] | 355 | dst_grid_pmask = xr.DataArray (oce_grid_pmask, dims = ['dst_grid_size',] ) |
---|
| 356 | dst_grid_pmask.attrs['long_name']="Land-sea mask (periodicity removed)" ; dst_grid_pmask.attrs['units']="Land:1, Ocean:0" |
---|
[4187] | 357 | |
---|
[6066] | 358 | src_lon_addressed = xr.DataArray (atm_grid_center_lon.values[atm_address] , dims = ['num_links'] ) |
---|
| 359 | src_lat_addressed = xr.DataArray (atm_grid_center_lat.values[atm_address] , dims = ['num_links'] ) |
---|
| 360 | src_area_addressed = xr.DataArray (atm_grid_area .values[atm_address] , dims = ['num_links'] ) |
---|
| 361 | src_imask_addressed = xr.DataArray (1-atm_grid_imask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 362 | src_pmask_addressed = xr.DataArray (1-atm_grid_pmask .values[atm_address].astype(np.int32) , dims = ['num_links'] ) |
---|
[4187] | 363 | |
---|
[6066] | 364 | dst_lon_addressed = xr.DataArray (oce_grid_center_lon.values[atm_address], dims = ['num_links'] ) |
---|
| 365 | dst_lat_addressed = xr.DataArray (oce_grid_center_lat.values[oce_address], dims = ['num_links'] ) |
---|
| 366 | dst_area_addressed = xr.DataArray (oce_grid_area.values[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 367 | dst_imask_addressed = xr.DataArray (1-oce_grid_imask[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 368 | dst_pmask_addressed = xr.DataArray (1-oce_grid_pmask[oce_address].astype(np.int32) , dims = ['num_links'] ) |
---|
| 369 | |
---|
| 370 | src_lon_addressed.attrs['long_name']="Longitude" ; src_lon_addressed.attrs['standard_name']="longitude" ; src_lon_addressed.attrs['units']="degrees_east" |
---|
| 371 | src_lat_addressed.attrs['long_name']="Latitude" ; src_lat_addressed.attrs['standard_name']="latitude" ; src_lat_addressed.attrs['units']="degrees_north" |
---|
| 372 | |
---|
| 373 | dst_lon_addressed.attrs['long_name']="Longitude" ; dst_lon_addressed.attrs['standard_name']="longitude" ; dst_lon_addressed.attrs['units']="degrees_east" |
---|
| 374 | dst_lat_addressed.attrs['long_name']="Latitude" ; dst_lat_addressed.attrs['standard_name']="latitude" ; dst_lat_addressed.attrs['units']="degrees_north" |
---|
| 375 | |
---|
[4187] | 376 | if atmDomainType == 'rectilinear' : |
---|
[6066] | 377 | atmLand = xr.DataArray ( atmLand.ravel() , dims = ['src_grid_size',] ) |
---|
| 378 | atmLandFiltered = xr.DataArray ( atmLandFrac.ravel() , dims = ['src_grid_size',] ) |
---|
| 379 | atmLandFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) |
---|
| 380 | atmFrac = xr.DataArray ( atmFrac.ravel() , dims = ['src_grid_size',] ) |
---|
| 381 | atmOcean = xr.DataArray ( atmOcean.ravel() , dims = ['src_grid_size',] ) |
---|
| 382 | atmOceanFrac = xr.DataArray ( atmOceanFrac.ravel(), dims = ['src_grid_size',] ) |
---|
[4187] | 383 | |
---|
[6066] | 384 | atmCoast = xr.DataArray (atmCoast.astype(np.int32) , dims = ['src_grid_size',]) |
---|
| 385 | oceLand = xr.DataArray (oceLand.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
| 386 | oceOcean = xr.DataArray (oceOcean.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
| 387 | oceOceanFiltered = xr.DataArray (oceOceanFiltered.astype(np.float32), dims = ['dst_grid_size',]) |
---|
| 388 | oceCoast = xr.DataArray (oceCoast.astype(np.int32) , dims = ['dst_grid_size',]) |
---|
[4187] | 389 | |
---|
| 390 | |
---|
[6066] | 391 | f_runoff = xr.Dataset ( { |
---|
[6090] | 392 | 'remap_matrix' : remap_matrix, |
---|
[6066] | 393 | 'src_address' : src_address, |
---|
| 394 | 'dst_address' : dst_address, |
---|
| 395 | 'src_grid_dims' : src_grid_dims, |
---|
| 396 | 'src_grid_center_lon' : src_grid_center_lon, |
---|
| 397 | 'src_grid_center_lat' : src_grid_center_lat, |
---|
| 398 | 'src_grid_corner_lon' : src_grid_corner_lon, |
---|
| 399 | 'src_grid_corner_lat' : src_grid_corner_lat, |
---|
| 400 | 'src_grid_area' : src_grid_area, |
---|
| 401 | 'src_grid_area' : src_grid_area, |
---|
[6314] | 402 | 'src_grid_imask' : src_grid_imask, |
---|
[6066] | 403 | 'src_grid_pmask' : src_grid_pmask, |
---|
| 404 | 'dst_grid_dims' : dst_grid_dims, |
---|
| 405 | 'dst_grid_center_lon' : dst_grid_center_lon, |
---|
[6314] | 406 | 'dst_grid_center_lat' : dst_grid_center_lat, |
---|
[6066] | 407 | 'dst_grid_corner_lon' : dst_grid_corner_lon, |
---|
| 408 | 'dst_grid_corner_lat' : dst_grid_corner_lat, |
---|
| 409 | 'dst_grid_area' : dst_grid_area, |
---|
| 410 | 'dst_grid_imask' : dst_grid_imask, |
---|
| 411 | 'dst_grid_pmask' : dst_grid_pmask, |
---|
| 412 | 'src_lon_addressed' : src_lon_addressed, |
---|
| 413 | 'src_lat_addressed' : src_lat_addressed, |
---|
| 414 | 'src_area_addressed' : src_area_addressed, |
---|
| 415 | 'dst_lon_addressed' : dst_lon_addressed, |
---|
| 416 | 'dst_lat_addressed' : dst_lat_addressed, |
---|
| 417 | 'dst_area_addressed' : dst_area_addressed, |
---|
| 418 | 'dst_imask_addressed' : dst_imask_addressed, |
---|
| 419 | 'dst_pmask_addressed' : dst_pmask_addressed, |
---|
| 420 | 'atmCoast' : atmCoast, |
---|
| 421 | 'oceLand' : oceLand, |
---|
| 422 | 'oceOcean' : oceOcean, |
---|
| 423 | 'oceOceanFiltered' : oceOceanFiltered, |
---|
| 424 | 'oceCoast' : oceCoast |
---|
| 425 | } ) |
---|
[4187] | 426 | |
---|
[6066] | 427 | f_runoff.attrs['Conventions'] = "CF-1.6" |
---|
| 428 | f_runoff.attrs['source'] = "IPSL Earth system model" |
---|
| 429 | f_runoff.attrs['group'] = "ICMC IPSL Climate Modelling Center" |
---|
| 430 | f_runoff.attrs['Institution'] = "IPSL https.//www.ipsl.fr" |
---|
| 431 | f_runoff.attrs['Ocean'] = oce_Name + " https://www.nemo-ocean.eu" |
---|
| 432 | f_runoff.attrs['Atmosphere'] = atm_Name + " http://lmdz.lmd.jussieu.fr" |
---|
| 433 | f_runoff.attrs['associatedFiles'] = grids + " " + areas + " " + masks |
---|
| 434 | f_runoff.attrs['description'] = "Generated with RunoffWeights.py" |
---|
| 435 | f_runoff.attrs['title'] = runoff |
---|
[6090] | 436 | f_runoff.attrs['Program'] = "Generated by " + sys.argv[0] + " with flags " + ' '.join (sys.argv[1:]) |
---|
[6066] | 437 | f_runoff.attrs['atmCoastWidth'] = "{:d} grid points".format(atmCoastWidth) |
---|
| 438 | f_runoff.attrs['oceCoastWidth'] = "{:d} grid points".format(oceCoastWidth) |
---|
| 439 | f_runoff.attrs['searchRadius'] = "{:.0f} km".format(searchRadius/1000.) |
---|
| 440 | f_runoff.attrs['atmQuantity'] = myargs.atmQuantity |
---|
| 441 | f_runoff.attrs['oceQuantity'] = myargs.oceQuantity |
---|
| 442 | f_runoff.attrs['gridsFile'] = grids |
---|
| 443 | f_runoff.attrs['areasFile'] = areas |
---|
| 444 | f_runoff.attrs['masksFile'] = masks |
---|
| 445 | f_runoff.attrs['o2aFile'] = o2a |
---|
[6090] | 446 | f_runoff.attrs['timeStamp'] = time.asctime () |
---|
| 447 | try : f_calving.attrs['directory'] = os.getcwd () |
---|
| 448 | except : pass |
---|
| 449 | try : f_runoff.attrs['HOSTNAME'] = platform.node () |
---|
| 450 | except : pass |
---|
| 451 | try : f_runoff.attrs['LOGNAME'] = os.getlogin () |
---|
| 452 | except : pass |
---|
| 453 | try : f_runoff.attrs['Python'] = "Python version " + platform.python_version () |
---|
| 454 | except : pass |
---|
| 455 | try : f_runoff.attrs['OS'] = platform.system () |
---|
| 456 | except : pass |
---|
| 457 | try : f_runoff.attrs['release'] = platform.release () |
---|
| 458 | except : pass |
---|
| 459 | try : f_runoff.attrs['hardware'] = platform.machine () |
---|
| 460 | except : pass |
---|
[6066] | 461 | f_runoff.attrs['conventions'] = "SCRIP" |
---|
| 462 | f_runoff.attrs['source_grid'] = "curvilinear" |
---|
| 463 | f_runoff.attrs['dest_grid'] = "curvilinear" |
---|
| 464 | f_runoff.attrs['Model'] = "IPSL CM6" |
---|
| 465 | f_runoff.attrs['SVN_Author'] = "$Author$" |
---|
| 466 | f_runoff.attrs['SVN_Date'] = "$Date$" |
---|
| 467 | f_runoff.attrs['SVN_Revision'] = "$Revision$" |
---|
| 468 | f_runoff.attrs['SVN_Id'] = "$Id$" |
---|
| 469 | f_runoff.attrs['SVN_HeadURL'] = "$HeadURL$" |
---|
| 470 | |
---|
| 471 | f_runoff.to_netcdf ( runoff, mode='w', format=FmtNetcdf ) |
---|
[6042] | 472 | f_runoff.close () |
---|
[4187] | 473 | |
---|
[6066] | 474 | ## |
---|
| 475 | |
---|
[6042] | 476 | print ('That''s all folks !') |
---|
[4187] | 477 | ## ====================================================================================== |
---|