| 1 | #!/usr/bin/env python3 |
|---|
| 2 | # -*- coding: utf-8 -*- |
|---|
| 3 | |
|---|
| 4 | #------------------------------------------------------------------------------# |
|---|
| 5 | # |
|---|
| 6 | # Scripts for processing of WRF and CAMx files to PALM dynamic driver |
|---|
| 7 | # |
|---|
| 8 | # This program is free software: you can redistribute it and/or modify |
|---|
| 9 | # it under the terms of the GNU General Public License as published by |
|---|
| 10 | # the Free Software Foundation, either version 3 of the License, or |
|---|
| 11 | # (at your option) any later version. |
|---|
| 12 | # |
|---|
| 13 | # This program is distributed in the hope that it will be useful, |
|---|
| 14 | # but WITHOUT ANY WARRANTY; without even the implied warranty of |
|---|
| 15 | # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
|---|
| 16 | # GNU General Public License for more details. |
|---|
| 17 | # |
|---|
| 18 | # You should have received a copy of the GNU General Public License |
|---|
| 19 | # along with this program. If not, see <https://www.gnu.org/licenses/>. |
|---|
| 20 | # |
|---|
| 21 | # Copyright 2018-2021 Institute of Computer Science |
|---|
| 22 | # of the Czech Academy of Sciences, Prague |
|---|
| 23 | # Authors: Krystof Eben, Jaroslav Resler, Pavel Krc |
|---|
| 24 | # |
|---|
| 25 | #------------------------------------------------------------------------------# |
|---|
| 26 | '''Scripts for processing of WRF and CAMx files to PALM dynamic driver. |
|---|
| 27 | |
|---|
| 28 | Usage: palm_dynamic -c <config_name> [-w] |
|---|
| 29 | |
|---|
| 30 | Script requires name of the configuration on command line. |
|---|
| 31 | The corresponding case configuration files are placed in subdirectory |
|---|
| 32 | "configuration" and the are named <config_name>.conf. The template |
|---|
| 33 | of the config file is in the file palm_dynamic_defaults.py, the values |
|---|
| 34 | which agree with defaults need not be present in the user config. |
|---|
| 35 | The file palm_dynamic_init.py contains setting and calculation of |
|---|
| 36 | standard initialization values for particular system and can be adjusted. |
|---|
| 37 | The optional parameter -w allows to skip horizontal and vertical |
|---|
| 38 | interpolation in case it is already done. |
|---|
| 39 | ''' |
|---|
| 40 | __version__ = '3.1' |
|---|
| 41 | |
|---|
| 42 | import sys |
|---|
| 43 | import getopt |
|---|
| 44 | import os.path |
|---|
| 45 | from datetime import datetime, timedelta |
|---|
| 46 | import numpy as np |
|---|
| 47 | from pyproj import Proj, transform |
|---|
| 48 | import glob |
|---|
| 49 | |
|---|
| 50 | import netCDF4 |
|---|
| 51 | |
|---|
| 52 | import palm_dynamic_config |
|---|
| 53 | |
|---|
| 54 | ################################## |
|---|
| 55 | # read configuration from command |
|---|
| 56 | # line and parse it |
|---|
| 57 | def print_help(): |
|---|
| 58 | print() |
|---|
| 59 | print(__doc__) |
|---|
| 60 | print('Version: '+__version__) |
|---|
| 61 | |
|---|
| 62 | # set commandline variables |
|---|
| 63 | configname = '' |
|---|
| 64 | wrf_interpolation = True |
|---|
| 65 | if sys.argv[0].endswith('pydevconsole.py'): |
|---|
| 66 | # we are in Pycharm developlent console - set testing data |
|---|
| 67 | configname = 'evropska_12s_basecase' |
|---|
| 68 | else: |
|---|
| 69 | try: |
|---|
| 70 | opts, args = getopt.getopt(sys.argv[1:],"hc:w)",["help=","config=","write="]) |
|---|
| 71 | #print(opts) |
|---|
| 72 | except getopt.GetoptError as err: |
|---|
| 73 | print("Error:", err) |
|---|
| 74 | print_help() |
|---|
| 75 | sys.exit(2) |
|---|
| 76 | # parse options |
|---|
| 77 | for opt, arg in opts: |
|---|
| 78 | if opt in ("-h", "--help"): |
|---|
| 79 | print_help() |
|---|
| 80 | sys.exit(0) |
|---|
| 81 | elif opt in ("-c", "--config"): |
|---|
| 82 | configname = arg |
|---|
| 83 | elif opt in ("-w", "--write"): |
|---|
| 84 | wrf_interpolation = False |
|---|
| 85 | |
|---|
| 86 | if configname == '': |
|---|
| 87 | # script requires config name |
|---|
| 88 | print("This script requires input parameter -c <config_name>") |
|---|
| 89 | print_help() |
|---|
| 90 | exit(2) |
|---|
| 91 | |
|---|
| 92 | # Load config defaults, config file and standard init into the config module |
|---|
| 93 | palm_dynamic_config.configure(configname) |
|---|
| 94 | # Import values (including loaded) from config module into globals |
|---|
| 95 | from palm_dynamic_config import * |
|---|
| 96 | |
|---|
| 97 | # Other modules are imported *AFTER* the config module has been filled with |
|---|
| 98 | # values (so that they can correctly import * from palm_dynamic_config). |
|---|
| 99 | import palm_wrf_utils |
|---|
| 100 | from palm_dynamic_output import palm_dynamic_output |
|---|
| 101 | import palm_dynamic_camx |
|---|
| 102 | |
|---|
| 103 | |
|---|
| 104 | print('Domain name: ', domain) |
|---|
| 105 | print('Resolution name: ', resolution) |
|---|
| 106 | print('Scenario name: ', scenario) |
|---|
| 107 | print('Read domain from static:', grid_from_static) |
|---|
| 108 | if grid_from_static: |
|---|
| 109 | print('STATIC driver input file: ', static_driver_file) |
|---|
| 110 | print('WRF and CAMx file path:', wrf_dir_name) |
|---|
| 111 | print('WRF file mask:', wrf_file_mask) |
|---|
| 112 | print('CAMx file mask:', camx_file_mask) |
|---|
| 113 | print('Simulation start time:', origin_time) |
|---|
| 114 | print('Simulation hours:', simulation_hours) |
|---|
| 115 | |
|---|
| 116 | |
|---|
| 117 | ########################### |
|---|
| 118 | # domain parameters |
|---|
| 119 | |
|---|
| 120 | if grid_from_static: |
|---|
| 121 | # get parameters of the horizontal domain from PALM STATIC driver |
|---|
| 122 | try: |
|---|
| 123 | ncs = netCDF4.Dataset(static_driver_file, "r", format="NETCDF4") |
|---|
| 124 | except Exception as err: |
|---|
| 125 | print("Error opening static driver file:") |
|---|
| 126 | print(static_driver_file) |
|---|
| 127 | print(err) |
|---|
| 128 | sys.exit(1) |
|---|
| 129 | # get horizontal structure of the domain |
|---|
| 130 | nx = ncs.dimensions['x'].size |
|---|
| 131 | ny = ncs.dimensions['y'].size |
|---|
| 132 | dx = ncs.variables['x'][:][1] - ncs.variables['x'][:][0] |
|---|
| 133 | dy = ncs.variables['y'][:][1] - ncs.variables['y'][:][0] |
|---|
| 134 | origin_x = ncs.getncattr('origin_x') |
|---|
| 135 | origin_y = ncs.getncattr('origin_y') |
|---|
| 136 | origin_z = ncs.getncattr('origin_z') |
|---|
| 137 | if origin_time == '': |
|---|
| 138 | # origin_time is not prvided in configuration - read it from static driver |
|---|
| 139 | origin_time = ncs.getncattr('origin_time') |
|---|
| 140 | |
|---|
| 141 | # create vertical structure of the domain |
|---|
| 142 | # calculate dz in case dz is not supplied (dz<=0) |
|---|
| 143 | if dz <= 0.0: |
|---|
| 144 | print("dz = 0.0: set dz = dx") |
|---|
| 145 | dz = dx |
|---|
| 146 | if nz <= 0: |
|---|
| 147 | print("nz > 0 needs to be set in config") |
|---|
| 148 | sys.exit(1) |
|---|
| 149 | |
|---|
| 150 | # read terrain height (relative to origin_z) and |
|---|
| 151 | # calculate and check the height of the surface canopy layer |
|---|
| 152 | if 'zt' in ncs.variables.keys(): |
|---|
| 153 | terrain_rel = ncs.variables['zt'][:] |
|---|
| 154 | else: |
|---|
| 155 | terrain_rel = np.zeros([ny,nx]) |
|---|
| 156 | th = np.ceil(terrain_rel / dz) |
|---|
| 157 | # building height |
|---|
| 158 | if 'buildings_3d' in ncs.variables.keys(): |
|---|
| 159 | print(np.argmax(a != 0, axis=0)) |
|---|
| 160 | bh3 = ncs.variables['buildings_3d'][:] |
|---|
| 161 | # minimum index of nonzeo value along inverted z |
|---|
| 162 | bh = np.argmax(bh3[::-1], axis=0) |
|---|
| 163 | # inversion back and masking grids with no buildings |
|---|
| 164 | bh = bh3.shape[0] - bh |
|---|
| 165 | bh[np.max(bh3, axis=0) == 0] = 0 |
|---|
| 166 | elif 'buildings_2d' in ncs.variables.keys(): |
|---|
| 167 | bh = ncs.variables['buildings_2d'][:] |
|---|
| 168 | bh[bh.mask] = 0 |
|---|
| 169 | bh = np.ceil(bh / dz) |
|---|
| 170 | else: |
|---|
| 171 | bh = np.zeros([ny,nx]) |
|---|
| 172 | # plant canopy height |
|---|
| 173 | if 'lad' in ncs.variables.keys(): |
|---|
| 174 | lad3 = ncs.variables['lad'][:] |
|---|
| 175 | # replace non-zero values with 1 |
|---|
| 176 | lad3[lad3 != 0] = 1 |
|---|
| 177 | # minimum index of nonzeo value along inverted z |
|---|
| 178 | lad = np.argmax(lad3[::-1], axis=0) |
|---|
| 179 | # inversion back and masking grids with no buildings |
|---|
| 180 | lad = lad3.shape[0] - lad |
|---|
| 181 | lad[np.max(lad3, axis=0) == 0] = 0 |
|---|
| 182 | else: |
|---|
| 183 | lad = np.zeros([ny,nx]) |
|---|
| 184 | # calculate maximum of surface canopy layer |
|---|
| 185 | nscl = max(np.amax(th+bh),np.amax(th+lad)) |
|---|
| 186 | # check nz with ncl |
|---|
| 187 | if nz < nscl + nscl_free: |
|---|
| 188 | print('nz has to be higher than ', nscl + nscl_free) |
|---|
| 189 | print('nz=', nz, ', dz=', dz, ', number of scl=', nscl, ',nscl_free=', nscl_free) |
|---|
| 190 | if dz_stretch_level <= (nscl + nscl_free) * dz: |
|---|
| 191 | print('stretching has to start in higher level than ', (nscl + nscl_free) * dz) |
|---|
| 192 | print('dz_stretch_level=', dz_stretch_level, ', nscl=', nscl, ', nscl_free=', nscl_free, ', dz=', dz) |
|---|
| 193 | if 'soil_moisture_adjust' in ncs.variables.keys(): |
|---|
| 194 | soil_moisture_adjust = ncs.variables['soil_moisture_adjust'][:] |
|---|
| 195 | else: |
|---|
| 196 | soil_moisture_adjust = np.ones(shape=(ny, nx), dtype=float) |
|---|
| 197 | # close static driver nc file |
|---|
| 198 | ncs.close() |
|---|
| 199 | else: |
|---|
| 200 | #TODO check all necessary parameters are set and are correct |
|---|
| 201 | # |
|---|
| 202 | # initialize necessary arrays |
|---|
| 203 | try: |
|---|
| 204 | terrain = np.zeros(shape=(ny, nx), dtype=float) |
|---|
| 205 | soil_moisture_adjust = np.ones(shape=(ny, nx), dtype=float) |
|---|
| 206 | except: |
|---|
| 207 | print('domain parameters nx, ny have to be set in configuration') |
|---|
| 208 | sys.exit(1) |
|---|
| 209 | |
|---|
| 210 | # absolute terrain needed for vertical interpolation of wrf data |
|---|
| 211 | terrain = terrain_rel + origin_z |
|---|
| 212 | |
|---|
| 213 | # print domain parameters and check ist existence in caso of setup from config |
|---|
| 214 | try: |
|---|
| 215 | print('Domain parameters:') |
|---|
| 216 | print('nx, ny, nz:', nx, ny, nz) |
|---|
| 217 | print('dx, dy, dz:', dx, dy, dz) |
|---|
| 218 | print('origin_x, origin_y:', origin_x, origin_y) |
|---|
| 219 | print('Base of domain is in level origin_z:', origin_z) |
|---|
| 220 | except: |
|---|
| 221 | print('domain parameters have to be read from static driver or set in configuration') |
|---|
| 222 | sys.exit(1) |
|---|
| 223 | |
|---|
| 224 | # centre of the domain (needed for ug,vg calculation) |
|---|
| 225 | xcent = origin_x + nx * dx / 2.0 |
|---|
| 226 | ycent = origin_y + ny * dy / 2.0 |
|---|
| 227 | # WGS84 projection for transformation to lat-lon |
|---|
| 228 | inproj = Proj('+init='+proj_palm) |
|---|
| 229 | print('inproj', inproj) |
|---|
| 230 | lonlatproj = Proj('+init='+proj_wgs84) |
|---|
| 231 | print('lonlatproj', lonlatproj) |
|---|
| 232 | cent_lon, cent_lat = transform(inproj, lonlatproj, xcent, ycent) |
|---|
| 233 | print('xcent, ycent:',xcent, ycent) |
|---|
| 234 | print('cent_lon, cent_lat:', cent_lon, cent_lat) |
|---|
| 235 | # prepare target grid |
|---|
| 236 | irange = origin_x + dx * (np.arange(nx, dtype='f8') + .5) |
|---|
| 237 | jrange = origin_y + dy * (np.arange(ny, dtype='f8') + .5) |
|---|
| 238 | palm_grid_y, palm_grid_x = np.meshgrid(jrange, irange, indexing='ij') |
|---|
| 239 | palm_grid_lon, palm_grid_lat = transform(inproj, lonlatproj, palm_grid_x, palm_grid_y) |
|---|
| 240 | |
|---|
| 241 | ###################################### |
|---|
| 242 | # build structure of vertical layers |
|---|
| 243 | # remark: |
|---|
| 244 | # PALM input requires nz=ztop in PALM |
|---|
| 245 | # but the output file in PALM has max z higher than z in PARIN. |
|---|
| 246 | # The highest levels in PALM are wrongly initialized !!! |
|---|
| 247 | ##################################### |
|---|
| 248 | # fill out z_levels |
|---|
| 249 | z_levels = np.zeros(nz,dtype=float) |
|---|
| 250 | z_levels_stag = np.zeros(nz-1,dtype=float) |
|---|
| 251 | dzs = dz |
|---|
| 252 | z_levels[0] = dzs/2.0 |
|---|
| 253 | for i in range(nz-1): |
|---|
| 254 | z_levels[i+1] = z_levels[i] + dzs |
|---|
| 255 | z_levels_stag[i] = (z_levels[i+1]+z_levels[i])/2.0 |
|---|
| 256 | dzso = dzs |
|---|
| 257 | if z_levels[i+1] + dzs >= dz_stretch_level: |
|---|
| 258 | dzs = min(dzs * dz_stretch_factor, dz_max) |
|---|
| 259 | ztop = z_levels[-1] + dzs / 2. |
|---|
| 260 | print('z:',z_levels) |
|---|
| 261 | print('zw:',z_levels_stag) |
|---|
| 262 | |
|---|
| 263 | ###################################### |
|---|
| 264 | # get time extent of the PALM simulation |
|---|
| 265 | ##################################### |
|---|
| 266 | # get complete list of wrf files |
|---|
| 267 | wrf_file_list = glob.glob(os.path.join(wrf_dir_name, wrf_file_mask)) |
|---|
| 268 | # get simulation origin and final time as datetime |
|---|
| 269 | start_time = datetime.strptime(origin_time, '%Y-%m-%d %H:%M:%S') |
|---|
| 270 | end_time = start_time + timedelta(hours=simulation_hours) |
|---|
| 271 | end_time_rad = end_time |
|---|
| 272 | print('PALM simulation extent', start_time, end_time, simulation_hours) |
|---|
| 273 | if nested_domain: |
|---|
| 274 | print('Nested domain - process only initialization.') |
|---|
| 275 | print('Set end_time = start_time') |
|---|
| 276 | end_time = start_time |
|---|
| 277 | |
|---|
| 278 | # get wrf times and sort wrf files by time |
|---|
| 279 | print('Analyse WRF files dates:') |
|---|
| 280 | file_times = [] |
|---|
| 281 | for wrf_file in wrf_file_list: |
|---|
| 282 | # get real time from wrf file |
|---|
| 283 | nc_wrf = netCDF4.Dataset(wrf_file, "r", format="NETCDF4") |
|---|
| 284 | ta = nc_wrf.variables['Times'][:] |
|---|
| 285 | t = ta.tobytes().decode("utf-8") |
|---|
| 286 | td = datetime.strptime(t, '%Y-%m-%d_%H:%M:%S') |
|---|
| 287 | print(os.path.basename(wrf_file), ': ', td) |
|---|
| 288 | file_times.append((td,wrf_file)) |
|---|
| 289 | file_times.sort() |
|---|
| 290 | times = [] |
|---|
| 291 | wrf_files = [] |
|---|
| 292 | for tf in file_times: |
|---|
| 293 | if end_time is None or tf[0] <= end_time: |
|---|
| 294 | times.append(tf[0]) |
|---|
| 295 | wrf_files.append(tf[1]) |
|---|
| 296 | #print('PALM output times:') |
|---|
| 297 | #print('\n'.join('{}'.format(t) for t in times)) |
|---|
| 298 | print('PALM output times:', ', '.join('{}'.format(t) for t in times)) |
|---|
| 299 | |
|---|
| 300 | if not times.__contains__(start_time): |
|---|
| 301 | print('WRF files does not contain PALM origin_time timestep - cannot process!') |
|---|
| 302 | exit(1) |
|---|
| 303 | |
|---|
| 304 | if not times.__contains__(end_time): |
|---|
| 305 | print('WRF files does not contain PALM end_time timestep - cannot process!') |
|---|
| 306 | exit(1) |
|---|
| 307 | |
|---|
| 308 | start_index = times.index(start_time) |
|---|
| 309 | end_index = times.index(end_time) |
|---|
| 310 | |
|---|
| 311 | if not nested_domain and end_index-start_index < simulation_hours: |
|---|
| 312 | print('Number of WRF files does not aggre with number of simulation hours') |
|---|
| 313 | exit(1) |
|---|
| 314 | |
|---|
| 315 | print('PALM simulation timestep number are from ', start_index, ' to ', end_index, ' WRF file.') |
|---|
| 316 | |
|---|
| 317 | # create list of processed wrf files |
|---|
| 318 | wrf_files_proc = [] |
|---|
| 319 | for wf in wrf_files[start_index:end_index+1]: |
|---|
| 320 | wrf_files_proc.append(wf) |
|---|
| 321 | # id of process files |
|---|
| 322 | simul_id = domain + "_d" + resolution |
|---|
| 323 | |
|---|
| 324 | ###################################### |
|---|
| 325 | #VERTICAL AND HORIZONTAL INTERPOLATION |
|---|
| 326 | ###################################### |
|---|
| 327 | # get hydro and soil variables contained in wrf files |
|---|
| 328 | shvars = [] |
|---|
| 329 | z_soil_levels = [] |
|---|
| 330 | nc_wrf = netCDF4.Dataset(wrf_files_proc[0], "r", format="NETCDF4") |
|---|
| 331 | try: |
|---|
| 332 | # hydro variables |
|---|
| 333 | shvars = sorted(set(['QVAPOR', 'QCLOUD', 'QRAIN', 'QICE', 'QSNOW', 'QGRAUP']).intersection(nc_wrf.variables.keys())) |
|---|
| 334 | # soil layers |
|---|
| 335 | if 'ZS' in nc_wrf.variables.keys(): |
|---|
| 336 | z_soil_levels = nc_wrf.variables['ZS'][0,:].data.tolist() |
|---|
| 337 | finally: |
|---|
| 338 | nc_wrf.close() |
|---|
| 339 | print('Hydro variables in wrf files:', '+'.join(shvars)) |
|---|
| 340 | |
|---|
| 341 | print('Start of vertical and horizontal interpolation of inputs to the PALM domain.') |
|---|
| 342 | interp_files = [] |
|---|
| 343 | regridder = None |
|---|
| 344 | for wrf_file in wrf_files_proc: |
|---|
| 345 | print ("Input wrf file: ", wrf_file) |
|---|
| 346 | # |
|---|
| 347 | hinterp_file = wrf_file+"_"+simul_id+'.hinterp' |
|---|
| 348 | hinterp_log = wrf_file + "_" + simul_id + '.hinterp.log' |
|---|
| 349 | vinterp_file = wrf_file+"_"+simul_id+'.interp' |
|---|
| 350 | interp_files.append(vinterp_file) |
|---|
| 351 | if wrf_interpolation: |
|---|
| 352 | try: |
|---|
| 353 | os.remove(vinterp_file) |
|---|
| 354 | os.remove(hinterp_file) |
|---|
| 355 | os.remove(hinterp_log) |
|---|
| 356 | except: |
|---|
| 357 | pass |
|---|
| 358 | |
|---|
| 359 | print('Horizontal interpolation...') |
|---|
| 360 | f_wrf = netCDF4.Dataset(wrf_file, 'r') |
|---|
| 361 | try: |
|---|
| 362 | if not regridder: |
|---|
| 363 | trans_wrf = palm_wrf_utils.WRFCoordTransform(f_wrf) |
|---|
| 364 | palm_in_wrf_y, palm_in_wrf_x = trans_wrf.latlon_to_ji( |
|---|
| 365 | palm_grid_lat, palm_grid_lon) |
|---|
| 366 | regridder = palm_wrf_utils.BilinearRegridder( |
|---|
| 367 | palm_in_wrf_x, palm_in_wrf_y, preloaded=True) |
|---|
| 368 | regridder_u = palm_wrf_utils.BilinearRegridder( |
|---|
| 369 | palm_in_wrf_x+.5, palm_in_wrf_y, preloaded=True) |
|---|
| 370 | regridder_v = palm_wrf_utils.BilinearRegridder( |
|---|
| 371 | palm_in_wrf_x, palm_in_wrf_y+.5, preloaded=True) |
|---|
| 372 | |
|---|
| 373 | f_out = netCDF4.Dataset(hinterp_file, 'w', format='NETCDF4') |
|---|
| 374 | try: |
|---|
| 375 | # dimensions |
|---|
| 376 | f_out.createDimension('Time', None) |
|---|
| 377 | for d in 'bottom_top bottom_top_stag soil_layers_stag'.split(): |
|---|
| 378 | f_out.createDimension(d, len(f_wrf.dimensions[d])) |
|---|
| 379 | f_out.createDimension('west_east', len(irange)) |
|---|
| 380 | f_out.createDimension('south_north', len(jrange)) |
|---|
| 381 | |
|---|
| 382 | # copied vars |
|---|
| 383 | for varname in 'PH PHB HGT T W TSLB SMOIS MU MUB P PB PSFC'.split(): |
|---|
| 384 | v_wrf = f_wrf.variables[varname] |
|---|
| 385 | v_out = f_out.createVariable(varname, 'f4', v_wrf.dimensions) |
|---|
| 386 | v_out[:] = regridder.regrid(v_wrf[...,regridder.ys,regridder.xs]) |
|---|
| 387 | |
|---|
| 388 | # U and V have special treatment (unstaggering) |
|---|
| 389 | v_out = f_out.createVariable('U', 'f4', ('Time', 'bottom_top', 'south_north', 'west_east')) |
|---|
| 390 | v_out[:] = regridder_u.regrid(f_wrf.variables['U'][...,regridder_u.ys,regridder_u.xs]) |
|---|
| 391 | v_out = f_out.createVariable('V', 'f4', ('Time', 'bottom_top', 'south_north', 'west_east')) |
|---|
| 392 | v_out[:] = regridder_v.regrid(f_wrf.variables['V'][...,regridder_v.ys,regridder_v.xs]) |
|---|
| 393 | |
|---|
| 394 | # calculated SPECHUM |
|---|
| 395 | v_out = f_out.createVariable('SPECHUM', 'f4', ('Time', 'bottom_top', 'south_north', 'west_east')) |
|---|
| 396 | vdata = regridder.regrid(f_wrf.variables[shvars[0]][...,regridder.ys,regridder.xs]) |
|---|
| 397 | for vname in shvars[1:]: |
|---|
| 398 | vdata += regridder.regrid(f_wrf.variables[vname][...,regridder.ys,regridder.xs]) |
|---|
| 399 | v_out[:] = vdata |
|---|
| 400 | finally: |
|---|
| 401 | f_out.close() |
|---|
| 402 | finally: |
|---|
| 403 | f_wrf.close() |
|---|
| 404 | |
|---|
| 405 | print('Vertical interpolation...') |
|---|
| 406 | palm_wrf_utils.palm_wrf_vertical_interp(hinterp_file, vinterp_file, wrf_file, z_levels, |
|---|
| 407 | z_levels_stag, z_soil_levels, origin_z, terrain, |
|---|
| 408 | wrf_hybrid_levs, vinterp_terrain_smoothing) |
|---|
| 409 | |
|---|
| 410 | if radiation_from_wrf: |
|---|
| 411 | print('Start processing of radiation inputs from the WRF radiation files.') |
|---|
| 412 | # get available times from wrf radiation output files and sort files by time |
|---|
| 413 | # get complete list of wrf radiation files |
|---|
| 414 | rad_file_list = glob.glob(os.path.join(wrf_dir_name, wrf_rad_file_mask)) |
|---|
| 415 | rad_file_times = [] |
|---|
| 416 | print('Analyse WRF radiation files dates:') |
|---|
| 417 | for rad_file in rad_file_list: |
|---|
| 418 | # get real time from wrf radiation file |
|---|
| 419 | nc_rad = netCDF4.Dataset(rad_file, "r", format="NETCDF4") |
|---|
| 420 | ta = nc_rad.variables['Times'][:] |
|---|
| 421 | t = ta.tobytes().decode("utf-8") |
|---|
| 422 | td = datetime.strptime(t, '%Y-%m-%d_%H:%M:%S') |
|---|
| 423 | print(os.path.basename(rad_file), ': ', td) |
|---|
| 424 | rad_file_times.append((td, rad_file)) |
|---|
| 425 | rad_file_times.sort() |
|---|
| 426 | rad_times = [] |
|---|
| 427 | rad_files = [] |
|---|
| 428 | for tf in rad_file_times: |
|---|
| 429 | if end_time_rad is None or tf[0] <= end_time_rad: |
|---|
| 430 | rad_times.append(tf[0]) |
|---|
| 431 | rad_files.append(tf[1]) |
|---|
| 432 | # check list of available times |
|---|
| 433 | if not rad_times.__contains__(start_time): |
|---|
| 434 | print('WRF radiation files does not contain PALM origin_time timestep - cannot process!') |
|---|
| 435 | exit(1) |
|---|
| 436 | if not rad_times.__contains__(end_time_rad): |
|---|
| 437 | print('WRF radiation files does not contain PALM end_time_rad timestep - cannot process!') |
|---|
| 438 | exit(1) |
|---|
| 439 | rad_start_index = rad_times.index(start_time) |
|---|
| 440 | rad_end_index = rad_times.index(end_time_rad) |
|---|
| 441 | print('PALM radiation timestep numbers are from ', rad_start_index, ' to ', rad_end_index, ' WRF radiation file.') |
|---|
| 442 | # get radiation timestep |
|---|
| 443 | rad_timestep = (rad_times[rad_start_index+1] - rad_times[rad_start_index]).total_seconds() |
|---|
| 444 | print('Timestep of the wrf radiation data is ', rad_timestep,' seconds.') |
|---|
| 445 | # check consistency of the series |
|---|
| 446 | for i in range(rad_start_index, rad_end_index): |
|---|
| 447 | if (rad_times[i+1] - rad_times[i]).total_seconds() != rad_timestep: |
|---|
| 448 | print('Inconsistent timeline of radiation inputs! Cannot process.') |
|---|
| 449 | print(rad_times[rad_start_index:rad_end_index]) |
|---|
| 450 | exit(1) |
|---|
| 451 | # create list of processed wrf radiation files and times |
|---|
| 452 | rad_files_proc = [] |
|---|
| 453 | rad_times_proc = [] |
|---|
| 454 | for i in range(rad_start_index, rad_end_index + 1): |
|---|
| 455 | rad_files_proc.append(rad_files[i]) |
|---|
| 456 | rad_times_proc.append((rad_times[i]-rad_times[rad_start_index]).total_seconds()) |
|---|
| 457 | print("Radiation times are ", rad_times_proc) |
|---|
| 458 | # process radiation inputs |
|---|
| 459 | rad_swdown = [] |
|---|
| 460 | rad_lwdown = [] |
|---|
| 461 | rad_swdiff = [] |
|---|
| 462 | nfiles = 0 |
|---|
| 463 | for rad_file in rad_files_proc: |
|---|
| 464 | print ("Input wrf radiation file: ", rad_file) |
|---|
| 465 | ncf = netCDF4.Dataset(rad_file, "r", format="NETCDF4") |
|---|
| 466 | try: |
|---|
| 467 | if nfiles == 0: |
|---|
| 468 | # create list of (i,j) indices used for smoothing of the radiation |
|---|
| 469 | print('Build list of indices for radiation smoothig.') |
|---|
| 470 | palmproj = Proj(init=proj_palm) |
|---|
| 471 | lonlatproj = Proj(init=proj_wgs84) |
|---|
| 472 | rad_ind = [] |
|---|
| 473 | for i in range(ncf.dimensions['west_east'].size): |
|---|
| 474 | for j in range(ncf.dimensions['south_north'].size): |
|---|
| 475 | # transfrom lat-lon to PALM domain in meters |
|---|
| 476 | lonij = ncf.variables['XLONG'][0,j,i] |
|---|
| 477 | latij = ncf.variables['XLAT'][0,j,i] |
|---|
| 478 | xij, yij = transform(lonlatproj, palmproj, lonij, latij) |
|---|
| 479 | if abs(xij-xcent)<=radiation_smoothing_distance and \ |
|---|
| 480 | abs(yij-ycent)<=radiation_smoothing_distance: |
|---|
| 481 | rad_ind.append([i,j]) |
|---|
| 482 | ngrids = len(rad_ind) |
|---|
| 483 | # process wrf radiation file |
|---|
| 484 | nfiles += 1 |
|---|
| 485 | swdown = 0.0 |
|---|
| 486 | lwd = 0.0 |
|---|
| 487 | swddif = 0.0 |
|---|
| 488 | for ij in rad_ind: |
|---|
| 489 | swdown = swdown + ncf.variables['SWDOWN'][0,ij[1],ij[0]] |
|---|
| 490 | lwd = lwd + ncf.variables['GLW'][0,ij[1],ij[0]] |
|---|
| 491 | swddif = swddif + ncf.variables['SWDDIF'][0,ij[1],ij[0]] |
|---|
| 492 | rad_swdown.append(swdown / ngrids) |
|---|
| 493 | rad_lwdown.append(lwd / ngrids) |
|---|
| 494 | rad_swdiff.append(swddif / ngrids) |
|---|
| 495 | finally: |
|---|
| 496 | ncf.close() |
|---|
| 497 | # create list of all radiation values |
|---|
| 498 | rad_values_proc = [rad_swdown, rad_lwdown, rad_swdiff] |
|---|
| 499 | else: |
|---|
| 500 | rad_times_proc = [] |
|---|
| 501 | rad_values_proc = [] |
|---|
| 502 | |
|---|
| 503 | camx_interp_fname = None |
|---|
| 504 | if not nested_domain: |
|---|
| 505 | # process camx files |
|---|
| 506 | camx_file_list = glob.glob(os.path.join(wrf_dir_name, camx_file_mask)) |
|---|
| 507 | if camx_file_list: |
|---|
| 508 | print('Processing CAMx input files: {0}'.format(', '.join(camx_file_list))) |
|---|
| 509 | camx_interp_fname = os.path.join(wrf_dir_name, 'CAMx_{0}_interp.nc'.format(simul_id)) |
|---|
| 510 | |
|---|
| 511 | palm_dynamic_camx.process_files(camx_file_list, camx_interp_fname, |
|---|
| 512 | palm_grid_lat, palm_grid_lon, terrain_rel, z_levels, |
|---|
| 513 | times[start_index:end_index+1], species_names, |
|---|
| 514 | camx_conversions, camx_helpers) |
|---|
| 515 | |
|---|
| 516 | # ===================== CREATE NETCDF DRIVER ============================== |
|---|
| 517 | # calculate relative times from simulation start |
|---|
| 518 | times_sec = [] |
|---|
| 519 | for t in times[start_index:end_index+1]: |
|---|
| 520 | times_sec.append((t-start_time).total_seconds()) |
|---|
| 521 | # collect dimension sizes |
|---|
| 522 | dimensions = {'zdim': nz, 'zwdim': nz-1, 'zsoildim': len(z_soil_levels), 'xdim': nx, 'xudim': nx-1, 'ydim': ny, 'yvdim': ny-1} |
|---|
| 523 | # process interpolated files to dynamic driver |
|---|
| 524 | palm_dynamic_output(wrf_files, interp_files, camx_interp_fname, dynamic_driver_file, times_sec, dimensions, |
|---|
| 525 | z_levels, z_levels_stag, ztop, z_soil_levels, dx, dy, cent_lon, cent_lat, |
|---|
| 526 | rad_times_proc, rad_values_proc, soil_moisture_adjust, nested_domain) |
|---|
| 527 | |
|---|
| 528 | print('Creation of dynamic driver finished.') |
|---|
