| 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 | ''' |
|---|
| 27 | This file creates and writes the dynamic driver netcdf file based on preprepared |
|---|
| 28 | transformed and interpolated wrf and camx files. |
|---|
| 29 | ''' |
|---|
| 30 | |
|---|
| 31 | import os |
|---|
| 32 | import time |
|---|
| 33 | import numpy as np |
|---|
| 34 | import netCDF4 |
|---|
| 35 | from palm_wrf_utils import palm_wrf_gw |
|---|
| 36 | |
|---|
| 37 | |
|---|
| 38 | def palm_dynamic_output(wrf_files, interp_files, camx_interp_fname, dynamic_driver_file, times_sec, |
|---|
| 39 | dimensions, z_levels, z_levels_stag, ztop, |
|---|
| 40 | z_soil_levels, dx, dy, lon_center, lat_center, |
|---|
| 41 | rad_times_proc, rad_values_proc, sma, nested_domain): |
|---|
| 42 | |
|---|
| 43 | print('Processing interpolated files to dynamic driver') |
|---|
| 44 | # dimension of the time coordinate |
|---|
| 45 | dimtimes = len(times_sec) |
|---|
| 46 | # other coordinates |
|---|
| 47 | dimnames = ['z', 'zw', 'zsoil', 'x','xu', 'y', 'yv'] # z height agl in m, zw staggered, zsoil 4 lev from wrf |
|---|
| 48 | dimsize_names = ['zdim' , 'zwdim', 'zsoildim', 'xdim', 'xudim', 'ydim', 'yvdim'] # palm: zw = z - 1 |
|---|
| 49 | x = np.arange(dx, dimensions['xdim']*dx+dx, dx) # clean this |
|---|
| 50 | print('dimension of x:', len(x)) |
|---|
| 51 | y = np.arange(dy, dimensions['ydim']*dy+dy, dy) # clean this |
|---|
| 52 | print('dimension of y:', len(y)) |
|---|
| 53 | # fill values |
|---|
| 54 | fillvalue_float = float(-9999.0) |
|---|
| 55 | # check out file and remove for sure |
|---|
| 56 | try: |
|---|
| 57 | os.remove(dynamic_driver_file) |
|---|
| 58 | except: |
|---|
| 59 | pass |
|---|
| 60 | # create netcdf out file |
|---|
| 61 | print('Driver file:', dynamic_driver_file) |
|---|
| 62 | outfile = netCDF4.Dataset(dynamic_driver_file, "w", format="NETCDF4" ) |
|---|
| 63 | try: |
|---|
| 64 | # time dimension and variable |
|---|
| 65 | outfile.createDimension('time', dimtimes) |
|---|
| 66 | _val_times = outfile.createVariable('time',"f4", ("time")) |
|---|
| 67 | # other dimensions and corresponding variables |
|---|
| 68 | for _dim in zip(dimnames, dimsize_names): #range (len(dimnames)): |
|---|
| 69 | print(_dim[0],_dim[1], dimensions[_dim[1]]) |
|---|
| 70 | outfile.createDimension(_dim[0], dimensions[_dim[1]]) |
|---|
| 71 | _val_z_levels = outfile.createVariable('z',"f4", ("z")) |
|---|
| 72 | _val_z_levels_stag = outfile.createVariable('zw',"f4", ("zw")) |
|---|
| 73 | _val_z_soil_levels = outfile.createVariable('zsoil',"f4", ("zsoil")) |
|---|
| 74 | _val_y = outfile.createVariable('y',"f4", ("y")) |
|---|
| 75 | _val_x = outfile.createVariable('x',"f4", ("x")) |
|---|
| 76 | |
|---|
| 77 | # prepare influx/outflux area sizes |
|---|
| 78 | zstag_all = np.r_[0., z_levels_stag, ztop] |
|---|
| 79 | zwidths = zstag_all[1:] - zstag_all[:-1] |
|---|
| 80 | print('zwidths', zwidths) |
|---|
| 81 | areas_xb = np.zeros((len(z_levels), 1)) |
|---|
| 82 | areas_xb[:,0] = zwidths * dy |
|---|
| 83 | areas_yb = np.zeros((len(z_levels), 1)) |
|---|
| 84 | areas_yb[:,0] = zwidths * dx |
|---|
| 85 | areas_zb = dx*dy |
|---|
| 86 | area_boundaries = (areas_xb.sum()*dimensions['ydim']*2 |
|---|
| 87 | + areas_yb.sum()*dimensions['xdim']*2 |
|---|
| 88 | + areas_zb*dimensions['xdim']*dimensions['ydim']) |
|---|
| 89 | |
|---|
| 90 | # write values for coordinates |
|---|
| 91 | _val_times[:] = times_sec[:] |
|---|
| 92 | _val_z_levels[:] = z_levels[:] |
|---|
| 93 | _val_z_levels_stag[:] = z_levels_stag[:] |
|---|
| 94 | _val_z_soil_levels[:] = z_soil_levels[:] |
|---|
| 95 | _val_y[:] = y[:] |
|---|
| 96 | _val_x[:] = x[:] |
|---|
| 97 | |
|---|
| 98 | # initialization of the variables and setting of init_* variables |
|---|
| 99 | print("Processing initialization from file", interp_files[0]) |
|---|
| 100 | # open corresponding |
|---|
| 101 | infile = netCDF4.Dataset(interp_files[0], "r", format="NETCDF4") |
|---|
| 102 | try: |
|---|
| 103 | # open variables in the input file |
|---|
| 104 | init_atmosphere_pt = infile.variables['init_atmosphere_pt'] |
|---|
| 105 | init_atmosphere_qv = infile.variables['init_atmosphere_qv'] |
|---|
| 106 | init_atmosphere_u = infile.variables['init_atmosphere_u'] |
|---|
| 107 | init_atmosphere_v = infile.variables['init_atmosphere_v'] |
|---|
| 108 | init_atmosphere_w = infile.variables['init_atmosphere_w'] |
|---|
| 109 | init_soil_m = infile.variables['init_soil_m'] |
|---|
| 110 | init_soil_t = infile.variables['init_soil_t'] |
|---|
| 111 | |
|---|
| 112 | # create netcdf structure |
|---|
| 113 | _val_init_atmosphere_pt = outfile.createVariable('init_atmosphere_pt', "f4", ("z", "y", "x"), |
|---|
| 114 | fill_value=fillvalue_float) |
|---|
| 115 | _val_init_atmosphere_pt.setncattr('lod', 2) |
|---|
| 116 | _val_init_atmosphere_qv = outfile.createVariable('init_atmosphere_qv', "f4", ("z", "y", "x"), |
|---|
| 117 | fill_value=fillvalue_float) |
|---|
| 118 | _val_init_atmosphere_qv.setncattr('lod', 2) |
|---|
| 119 | _val_init_atmosphere_u = outfile.createVariable('init_atmosphere_u', "f4", ("z", "y", "xu"), |
|---|
| 120 | fill_value=fillvalue_float) |
|---|
| 121 | _val_init_atmosphere_u.setncattr('lod', 2) |
|---|
| 122 | _val_init_atmosphere_v = outfile.createVariable('init_atmosphere_v', "f4", ("z", "yv", "x"), |
|---|
| 123 | fill_value=fillvalue_float) |
|---|
| 124 | _val_init_atmosphere_v.setncattr('lod', 2) |
|---|
| 125 | _val_init_atmosphere_w = outfile.createVariable('init_atmosphere_w', "f4", ("zw", "y", "x"), |
|---|
| 126 | fill_value=fillvalue_float) |
|---|
| 127 | _val_init_atmosphere_w.setncattr('lod', 2) |
|---|
| 128 | _val_init_soil_t = outfile.createVariable('init_soil_t', "f4", ("zsoil", "y", "x"), fill_value=fillvalue_float) |
|---|
| 129 | _val_init_soil_t.setncattr('lod', 2) |
|---|
| 130 | _val_init_soil_m = outfile.createVariable('init_soil_m', "f4", ("zsoil", "y", "x"), fill_value=fillvalue_float) |
|---|
| 131 | _val_init_soil_m.setncattr('lod', 2) |
|---|
| 132 | # time dependent variables |
|---|
| 133 | if not nested_domain: |
|---|
| 134 | # SURFACE PRESSURE |
|---|
| 135 | _val_surface_forcing_surface_pressure = outfile.createVariable('surface_forcing_surface_pressure', "f4", |
|---|
| 136 | ("time")) |
|---|
| 137 | # BOUNDARY - vertical slices from left, right, south, north, top |
|---|
| 138 | varname = 'pt' |
|---|
| 139 | _val_ls_forcing_pt_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "y"), |
|---|
| 140 | fill_value=fillvalue_float) |
|---|
| 141 | _val_ls_forcing_pt_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "y"), |
|---|
| 142 | fill_value=fillvalue_float) |
|---|
| 143 | _val_ls_forcing_pt_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "x"), |
|---|
| 144 | fill_value=fillvalue_float) |
|---|
| 145 | _val_ls_forcing_pt_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "x"), |
|---|
| 146 | fill_value=fillvalue_float) |
|---|
| 147 | _val_ls_forcing_pt_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "x"), |
|---|
| 148 | fill_value=fillvalue_float) |
|---|
| 149 | |
|---|
| 150 | varname = 'qv' |
|---|
| 151 | _val_ls_forcing_qv_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "y"), |
|---|
| 152 | fill_value=fillvalue_float) |
|---|
| 153 | _val_ls_forcing_qv_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "y"), |
|---|
| 154 | fill_value=fillvalue_float) |
|---|
| 155 | _val_ls_forcing_qv_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "x"), |
|---|
| 156 | fill_value=fillvalue_float) |
|---|
| 157 | _val_ls_forcing_qv_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "x"), |
|---|
| 158 | fill_value=fillvalue_float) |
|---|
| 159 | _val_ls_forcing_qv_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "x"), |
|---|
| 160 | fill_value=fillvalue_float) |
|---|
| 161 | |
|---|
| 162 | varname = 'u' |
|---|
| 163 | _val_ls_forcing_u_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "y"), |
|---|
| 164 | fill_value=fillvalue_float) |
|---|
| 165 | _val_ls_forcing_u_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "y"), |
|---|
| 166 | fill_value=fillvalue_float) |
|---|
| 167 | _val_ls_forcing_u_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "xu"), |
|---|
| 168 | fill_value=fillvalue_float) |
|---|
| 169 | _val_ls_forcing_u_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "xu"), |
|---|
| 170 | fill_value=fillvalue_float) |
|---|
| 171 | _val_ls_forcing_u_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "xu"), |
|---|
| 172 | fill_value=fillvalue_float) |
|---|
| 173 | |
|---|
| 174 | varname = 'v' |
|---|
| 175 | _val_ls_forcing_v_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "z", "yv"), |
|---|
| 176 | fill_value=fillvalue_float) |
|---|
| 177 | _val_ls_forcing_v_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "z", "yv"), |
|---|
| 178 | fill_value=fillvalue_float) |
|---|
| 179 | _val_ls_forcing_v_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "z", "x"), |
|---|
| 180 | fill_value=fillvalue_float) |
|---|
| 181 | _val_ls_forcing_v_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "z", "x"), |
|---|
| 182 | fill_value=fillvalue_float) |
|---|
| 183 | _val_ls_forcing_v_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "yv", "x"), |
|---|
| 184 | fill_value=fillvalue_float) |
|---|
| 185 | |
|---|
| 186 | varname = 'w' |
|---|
| 187 | _val_ls_forcing_w_left = outfile.createVariable('ls_forcing_left_' + varname, "f4", ("time", "zw", "y"), |
|---|
| 188 | fill_value=fillvalue_float) |
|---|
| 189 | _val_ls_forcing_w_right = outfile.createVariable('ls_forcing_right_' + varname, "f4", ("time", "zw", "y"), |
|---|
| 190 | fill_value=fillvalue_float) |
|---|
| 191 | _val_ls_forcing_w_south = outfile.createVariable('ls_forcing_south_' + varname, "f4", ("time", "zw", "x"), |
|---|
| 192 | fill_value=fillvalue_float) |
|---|
| 193 | _val_ls_forcing_w_north = outfile.createVariable('ls_forcing_north_' + varname, "f4", ("time", "zw", "x"), |
|---|
| 194 | fill_value=fillvalue_float) |
|---|
| 195 | _val_ls_forcing_w_top = outfile.createVariable('ls_forcing_top_' + varname, "f4", ("time", "y", "x"), |
|---|
| 196 | fill_value=fillvalue_float) |
|---|
| 197 | |
|---|
| 198 | # geostrophic wind |
|---|
| 199 | _val_ls_forcing_ug = outfile.createVariable('ls_forcing_ug', "f4", ("time", "z"), |
|---|
| 200 | fill_value=fillvalue_float) |
|---|
| 201 | _val_ls_forcing_vg = outfile.createVariable('ls_forcing_vg', "f4", ("time", "z"), |
|---|
| 202 | fill_value=fillvalue_float) |
|---|
| 203 | time.sleep(1) |
|---|
| 204 | |
|---|
| 205 | # write values for initialization variables |
|---|
| 206 | _val_init_atmosphere_pt[:, :, :] = init_atmosphere_pt[0, :, :, :] |
|---|
| 207 | _val_init_atmosphere_qv[:, :, :] = init_atmosphere_qv[0, :, :, :] |
|---|
| 208 | _val_init_atmosphere_u[:, :, :] = init_atmosphere_u[0, :, :, 1:] |
|---|
| 209 | _val_init_atmosphere_v[:, :, :] = init_atmosphere_v[0, :, 1:, :] |
|---|
| 210 | _val_init_atmosphere_w[:, :, :] = init_atmosphere_w[0, :, :, :] |
|---|
| 211 | _val_init_soil_t[:, :, :] = init_soil_t[0, :, :, :] |
|---|
| 212 | for k in range(0,_val_init_soil_m.shape[0]): |
|---|
| 213 | # adjust soil moisture according soil_moisture_adjust field (if exists) |
|---|
| 214 | _val_init_soil_m[k, :, :] = init_soil_m[0, k, :, :] * sma[:, :] |
|---|
| 215 | finally: |
|---|
| 216 | # close interpolated file |
|---|
| 217 | infile.close() |
|---|
| 218 | # |
|---|
| 219 | if not nested_domain: |
|---|
| 220 | # cycle over all included time steps |
|---|
| 221 | for ts in range(0, len(interp_files)): |
|---|
| 222 | print("Processing file",interp_files[ts]) |
|---|
| 223 | # open corresponding interpolated file |
|---|
| 224 | infile = netCDF4.Dataset(interp_files[ts], "r", format="NETCDF4") |
|---|
| 225 | try: |
|---|
| 226 | # open variables in the input file |
|---|
| 227 | init_atmosphere_pt = infile.variables['init_atmosphere_pt'] |
|---|
| 228 | init_atmosphere_qv = infile.variables['init_atmosphere_qv'] |
|---|
| 229 | init_atmosphere_u = infile.variables['init_atmosphere_u'] |
|---|
| 230 | init_atmosphere_v = infile.variables['init_atmosphere_v'] |
|---|
| 231 | init_atmosphere_w = infile.variables['init_atmosphere_w'] |
|---|
| 232 | surface_forcing_surface_pressure = infile.variables['surface_forcing_surface_pressure'] |
|---|
| 233 | |
|---|
| 234 | ################################## |
|---|
| 235 | # write values for time dependent values |
|---|
| 236 | # surface pressure |
|---|
| 237 | _val_surface_forcing_surface_pressure[ts] = np.average(surface_forcing_surface_pressure[:,:,:], axis = (1,2))[0] |
|---|
| 238 | # boundary conditions |
|---|
| 239 | _val_ls_forcing_pt_left[ts, :, :] = init_atmosphere_pt[0, :, :, 0] |
|---|
| 240 | _val_ls_forcing_pt_right[ts, :, :] = init_atmosphere_pt[0, :, :, dimensions['xdim'] - 1] |
|---|
| 241 | _val_ls_forcing_pt_south[ts, :, :] = init_atmosphere_pt[0, :, 0, :] |
|---|
| 242 | _val_ls_forcing_pt_north[ts, :, :] = init_atmosphere_pt[0, :, dimensions['ydim'] - 1, :] |
|---|
| 243 | _val_ls_forcing_pt_top[ts, :, :] = init_atmosphere_pt[0, dimensions['zdim'] - 1, :, :] |
|---|
| 244 | |
|---|
| 245 | _val_ls_forcing_qv_left[ts, :, :] = init_atmosphere_qv[0, :, :, 0] |
|---|
| 246 | _val_ls_forcing_qv_right[ts, :, :] = init_atmosphere_qv[0, :, :, dimensions['xdim'] - 1] |
|---|
| 247 | _val_ls_forcing_qv_south[ts, :, :] = init_atmosphere_qv[0, :, 0, :] |
|---|
| 248 | _val_ls_forcing_qv_north[ts, :, :] = init_atmosphere_qv[0, :, dimensions['ydim'] - 1, :] |
|---|
| 249 | _val_ls_forcing_qv_top[ts, :, :] = init_atmosphere_qv[0, dimensions['zdim'] - 1, :, :] |
|---|
| 250 | |
|---|
| 251 | # Perform mass balancing |
|---|
| 252 | uxleft = init_atmosphere_u[0, :, :, 0] |
|---|
| 253 | uxright = init_atmosphere_u[0, :, :, dimensions['xdim'] - 1] |
|---|
| 254 | vysouth = init_atmosphere_v[0, :, 0, :] |
|---|
| 255 | vynorth = init_atmosphere_v[0, :, dimensions['ydim'] - 1, :] |
|---|
| 256 | wztop = init_atmosphere_w[0, dimensions['zwdim'] - 1, :, :] |
|---|
| 257 | mass_disbalance = ((uxleft * areas_xb).sum() |
|---|
| 258 | - (uxright * areas_xb).sum() |
|---|
| 259 | + (vysouth * areas_yb).sum() |
|---|
| 260 | - (vynorth * areas_yb).sum() |
|---|
| 261 | - (wztop * areas_zb).sum()) |
|---|
| 262 | mass_corr_v = mass_disbalance / area_boundaries |
|---|
| 263 | print('Mass disbalance: {0:8g} m3/s (avg = {1:8g} m/s)'.format( |
|---|
| 264 | mass_disbalance, mass_corr_v)) |
|---|
| 265 | uxleft -= mass_corr_v |
|---|
| 266 | uxright += mass_corr_v |
|---|
| 267 | vysouth -= mass_corr_v |
|---|
| 268 | vynorth += mass_corr_v |
|---|
| 269 | wztop += mass_corr_v |
|---|
| 270 | |
|---|
| 271 | # Verify mass balance (optional) |
|---|
| 272 | #mass_disbalance = ((uxleft * areas_xb).sum() |
|---|
| 273 | # - (uxright * areas_xb).sum() |
|---|
| 274 | # + (vysouth * areas_yb).sum() |
|---|
| 275 | # - (vynorth * areas_yb).sum() |
|---|
| 276 | # - (wztop * areas_zb).sum()) |
|---|
| 277 | #mass_corr_v = mass_disbalance / area_boundaries |
|---|
| 278 | #print('Mass balanced: {0:8g} m3/s (avg = {1:8g} m/s)'.format( |
|---|
| 279 | # mass_disbalance, mass_corr_v)) |
|---|
| 280 | |
|---|
| 281 | _val_ls_forcing_u_left[ts, :, :] = uxleft |
|---|
| 282 | _val_ls_forcing_u_right[ts, :, :] = uxright |
|---|
| 283 | _val_ls_forcing_u_south[ts, :, :] = init_atmosphere_u[0, :, 0, 1:] |
|---|
| 284 | _val_ls_forcing_u_north[ts, :, :] = init_atmosphere_u[0, :, dimensions['ydim'] - 1, 1:] |
|---|
| 285 | _val_ls_forcing_u_top[ts, :, :] = init_atmosphere_u[0, dimensions['zdim'] - 1, :, 1:] |
|---|
| 286 | |
|---|
| 287 | _val_ls_forcing_v_left[ts, :, :] = init_atmosphere_v[0, :, 1:, 0] |
|---|
| 288 | _val_ls_forcing_v_right[ts, :, :] = init_atmosphere_v[0, :, 1:, dimensions['xdim'] - 1] |
|---|
| 289 | _val_ls_forcing_v_south[ts, :, :] = vysouth |
|---|
| 290 | _val_ls_forcing_v_north[ts, :, :] = vynorth |
|---|
| 291 | _val_ls_forcing_v_top[ts, :, :] = init_atmosphere_v[0, dimensions['zdim'] - 1, 1:, :] |
|---|
| 292 | |
|---|
| 293 | _val_ls_forcing_w_left[ts, :, :] = init_atmosphere_w[0, :, :, 0] |
|---|
| 294 | _val_ls_forcing_w_right[ts, :, :] = init_atmosphere_w[0, :, :, dimensions['xdim'] - 1] |
|---|
| 295 | _val_ls_forcing_w_south[ts, :, :] = init_atmosphere_w[0, :, 0, :] |
|---|
| 296 | _val_ls_forcing_w_north[ts, :, :] = init_atmosphere_w[0, :, dimensions['ydim'] - 1, :] |
|---|
| 297 | _val_ls_forcing_w_top[ts, :, :] = wztop |
|---|
| 298 | |
|---|
| 299 | finally: |
|---|
| 300 | # close interpolated file |
|---|
| 301 | infile.close() |
|---|
| 302 | |
|---|
| 303 | # write geostrophic wind |
|---|
| 304 | print('Open wrf file '+wrf_files[ts]) |
|---|
| 305 | nc_wrf = netCDF4.Dataset(wrf_files[ts], 'r') |
|---|
| 306 | try: |
|---|
| 307 | ug, vg = palm_wrf_gw(nc_wrf, lon_center, lat_center, z_levels) |
|---|
| 308 | _val_ls_forcing_ug[ts, :] = ug |
|---|
| 309 | _val_ls_forcing_vg[ts, :] = vg |
|---|
| 310 | finally: |
|---|
| 311 | nc_wrf.close() |
|---|
| 312 | |
|---|
| 313 | # Write chemical boundary conds |
|---|
| 314 | if camx_interp_fname: |
|---|
| 315 | f_camx = netCDF4.Dataset(camx_interp_fname) |
|---|
| 316 | try: |
|---|
| 317 | for vname, vval in f_camx.variables.items(): |
|---|
| 318 | # PALM doesn't support 3D LOD=2 init for chem yet, we have to average the field |
|---|
| 319 | var = outfile.createVariable('init_atmosphere_'+vname, |
|---|
| 320 | 'f4', ('z',), fill_value=fillvalue_float) |
|---|
| 321 | var.units = vval.units |
|---|
| 322 | var.lod = 1 |
|---|
| 323 | var[:] = vval[0,:,:,:].mean(axis=(1,2)) |
|---|
| 324 | |
|---|
| 325 | var = outfile.createVariable('ls_forcing_left_'+vname, |
|---|
| 326 | 'f4', ('time','z','y'), fill_value=fillvalue_float) |
|---|
| 327 | var.units = vval.units |
|---|
| 328 | var[:] = vval[:,:,:,0] |
|---|
| 329 | |
|---|
| 330 | var = outfile.createVariable('ls_forcing_right_'+vname, |
|---|
| 331 | 'f4', ('time','z','y'), fill_value=fillvalue_float) |
|---|
| 332 | var.units = vval.units |
|---|
| 333 | var[:] = vval[:,:,:,-1] |
|---|
| 334 | |
|---|
| 335 | var = outfile.createVariable('ls_forcing_south_'+vname, |
|---|
| 336 | 'f4', ('time','z','x'), fill_value=fillvalue_float) |
|---|
| 337 | var.units = vval.units |
|---|
| 338 | var[:] = vval[:,:,0,:] |
|---|
| 339 | |
|---|
| 340 | var = outfile.createVariable('ls_forcing_north_'+vname, |
|---|
| 341 | 'f4', ('time','z','x'), fill_value=fillvalue_float) |
|---|
| 342 | var.units = vval.units |
|---|
| 343 | var[:] = vval[:,:,-1,:] |
|---|
| 344 | |
|---|
| 345 | var = outfile.createVariable('ls_forcing_top_'+vname, |
|---|
| 346 | 'f4', ('time','y','x'), fill_value=fillvalue_float) |
|---|
| 347 | var.units = vval.units |
|---|
| 348 | var[:] = vval[:,-1,:,:] |
|---|
| 349 | finally: |
|---|
| 350 | f_camx.close() |
|---|
| 351 | |
|---|
| 352 | if len(rad_times_proc) > 0: |
|---|
| 353 | # process radiation inputs |
|---|
| 354 | # radiation time dimension and variable |
|---|
| 355 | outfile.createDimension('time_rad', len(rad_times_proc)) |
|---|
| 356 | _val_times = outfile.createVariable('time_rad',"f4", ("time_rad")) |
|---|
| 357 | _val_times[:] = rad_times_proc[:] |
|---|
| 358 | # radiation variables |
|---|
| 359 | var = outfile.createVariable('rad_sw_in', "f4", ("time_rad"), fill_value=fillvalue_float) |
|---|
| 360 | var.setncattr('lod', 1) |
|---|
| 361 | var.units = 'W/m2' |
|---|
| 362 | var[:] = rad_values_proc[0][:] |
|---|
| 363 | var = outfile.createVariable('rad_lw_in', "f4", ("time_rad"), fill_value=fillvalue_float) |
|---|
| 364 | var.setncattr('lod', 1) |
|---|
| 365 | var.units = 'W/m2' |
|---|
| 366 | var[:] = rad_values_proc[1][:] |
|---|
| 367 | var = outfile.createVariable('rad_sw_in_dif', "f4", ("time_rad"), fill_value=fillvalue_float) |
|---|
| 368 | var.setncattr('lod', 1) |
|---|
| 369 | var.units = 'W/m2' |
|---|
| 370 | var[:] = rad_values_proc[2][:] |
|---|
| 371 | |
|---|
| 372 | finally: |
|---|
| 373 | outfile.close() |
|---|
| 374 | |
|---|
| 375 | |
|---|
| 376 | |
|---|
| 377 | |
|---|
| 378 | |
|---|
