source: palm/trunk/UTIL/WRF_interface/dynamic/palm_dynamic_output.py @ 4843

Last change on this file since 4843 was 4843, checked in by raasch, 3 years ago

local namelist parameter added to switch off the module although the respective module namelist appears in the namelist file, further copyright updates

File size: 21.9 KB
Line 
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'''
27This file creates and writes the dynamic driver netcdf file based on preprepared
28transformed and interpolated wrf and camx files.
29'''
30
31import os
32import time
33import numpy as np
34import netCDF4
35from palm_wrf_utils import palm_wrf_gw
36
37
38def 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
Note: See TracBrowser for help on using the repository browser.