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

Last change on this file since 4843 was 4843, checked in by raasch, 4 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.1 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'''Scripts for processing of WRF and CAMx files to PALM dynamic driver.
27
28Usage: palm_dynamic -c <config_name> [-w]
29
30Script requires name of the configuration on command line.
31The corresponding case configuration files are placed in subdirectory
32"configuration" and the are named <config_name>.conf. The template
33of the config file is in the file palm_dynamic_defaults.py, the values
34which agree with defaults need not be present in the user config.
35The file palm_dynamic_init.py contains setting and calculation of
36standard initialization values for particular system and can be adjusted.
37The optional parameter -w allows to skip horizontal and vertical
38interpolation in case it is already done.
39'''
40__version__ = '3.1'
41
42import sys
43import getopt
44import os.path
45from datetime import datetime, timedelta
46import numpy as np
47from pyproj import Proj, transform
48import glob
49
50import netCDF4
51
52import palm_dynamic_config
53
54##################################
55# read configuration from command
56# line and parse it
57def print_help():
58    print()
59    print(__doc__)
60    print('Version: '+__version__)
61
62# set commandline variables
63configname = ''
64wrf_interpolation = True
65if sys.argv[0].endswith('pydevconsole.py'):
66    # we are in Pycharm developlent console - set testing data
67    configname = 'evropska_12s_basecase'
68else:
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
86if 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
93palm_dynamic_config.configure(configname)
94# Import values (including loaded) from config module into globals
95from 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).
99import palm_wrf_utils
100from palm_dynamic_output import palm_dynamic_output
101import palm_dynamic_camx
102
103
104print('Domain name: ', domain)
105print('Resolution name: ', resolution)
106print('Scenario name: ', scenario)
107print('Read domain from static:', grid_from_static)
108if grid_from_static:
109    print('STATIC driver input file: ', static_driver_file)
110print('WRF and CAMx file path:', wrf_dir_name)
111print('WRF file mask:', wrf_file_mask)
112print('CAMx file mask:', camx_file_mask)
113print('Simulation start time:', origin_time)
114print('Simulation hours:', simulation_hours)
115
116
117###########################
118# domain parameters
119
120if 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()
199else:
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
211terrain = terrain_rel + origin_z
212
213# print domain parameters and check ist existence in caso of setup from config
214try:
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)
220except:
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)
225xcent = origin_x + nx * dx / 2.0
226ycent = origin_y + ny * dy / 2.0
227# WGS84 projection for transformation to lat-lon
228inproj = Proj('+init='+proj_palm)
229print('inproj', inproj)
230lonlatproj = Proj('+init='+proj_wgs84)
231print('lonlatproj', lonlatproj)
232cent_lon, cent_lat = transform(inproj, lonlatproj, xcent, ycent)
233print('xcent, ycent:',xcent, ycent)
234print('cent_lon, cent_lat:', cent_lon, cent_lat)
235# prepare target grid
236irange = origin_x + dx * (np.arange(nx, dtype='f8') + .5)
237jrange = origin_y + dy * (np.arange(ny, dtype='f8') + .5)
238palm_grid_y, palm_grid_x = np.meshgrid(jrange, irange, indexing='ij')
239palm_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
249z_levels = np.zeros(nz,dtype=float)
250z_levels_stag = np.zeros(nz-1,dtype=float)
251dzs = dz
252z_levels[0] = dzs/2.0
253for 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)
259ztop = z_levels[-1] + dzs / 2.
260print('z:',z_levels)
261print('zw:',z_levels_stag)
262
263######################################
264# get time extent of the PALM simulation
265#####################################
266# get complete list of wrf files
267wrf_file_list = glob.glob(os.path.join(wrf_dir_name, wrf_file_mask))
268# get simulation origin and final time as datetime
269start_time = datetime.strptime(origin_time, '%Y-%m-%d %H:%M:%S')
270end_time = start_time + timedelta(hours=simulation_hours)
271end_time_rad = end_time
272print('PALM simulation extent', start_time, end_time, simulation_hours)
273if 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
279print('Analyse WRF files dates:')
280file_times = []
281for 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))
289file_times.sort()
290times = []
291wrf_files = []
292for 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))
298print('PALM output times:', ', '.join('{}'.format(t) for t in times))
299
300if not times.__contains__(start_time):
301    print('WRF files does not contain PALM origin_time timestep - cannot process!')
302    exit(1)
303
304if not times.__contains__(end_time):
305    print('WRF files does not contain PALM end_time timestep - cannot process!')
306    exit(1)
307
308start_index = times.index(start_time)
309end_index = times.index(end_time)
310
311if 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
315print('PALM simulation timestep number are from ', start_index, ' to ', end_index, ' WRF file.')
316
317# create list of processed wrf files
318wrf_files_proc = []
319for wf in wrf_files[start_index:end_index+1]:
320    wrf_files_proc.append(wf)
321# id of process files
322simul_id = domain + "_d" + resolution
323
324######################################
325#VERTICAL AND HORIZONTAL INTERPOLATION
326######################################
327# get hydro and soil variables contained in wrf files
328shvars = []
329z_soil_levels = []
330nc_wrf = netCDF4.Dataset(wrf_files_proc[0], "r", format="NETCDF4")
331try:
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()
337finally:
338    nc_wrf.close()
339print('Hydro variables in wrf files:', '+'.join(shvars))
340
341print('Start of vertical and horizontal interpolation of inputs to the PALM domain.')
342interp_files = []
343regridder = None
344for 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
410if 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]
499else:
500    rad_times_proc = []
501    rad_values_proc = []
502
503camx_interp_fname = None
504if 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
518times_sec = []
519for t in times[start_index:end_index+1]:
520    times_sec.append((t-start_time).total_seconds())
521# collect dimension sizes
522dimensions = {'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
524palm_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
528print('Creation of dynamic driver finished.')
Note: See TracBrowser for help on using the repository browser.