Update WRF namelist translation

README.md

@@ -1,21 +1,50 @@
-# ERFtools
-A collection of scripts for facilitating the usage of ERF.
+# ERF Tools
+A collection of Python-based modules and scripts for facilitating the usage of
+the Energy Research and Forecasting (ERF) model.
+
+## Installation
+
+Clone and create a conda (or mamba) environment:
+```shell
+git clone https://github.com/erf-model/erftools.git
+cd erftools
+conda create -n erftools python=3.9
+conda activate erftools
+```
+
+Then, install with:
+```shell
+pip install -e . # editable install
+```
+or install with all optional dependencies:
+```shell
+pip install -e .[all]
+```
 
 ## Examples
 
 ### Converting a WRF namelist into ERF inputs
 ```python
 from erftools.preprocessing import WRFInputDeck
 wrf = WRFInputDeck('namelist.input')
+wrf.write_inputfile('inputs')
+
+# optional: extract data from wrfinput
 wrf.process_initial_conditions('wrfinput_d01',
                                landuse_table_path='/Users/equon/WRF/run/LANDUSE.TBL',
                                write_hgt='terrain_height.txt',
                                write_z0='roughness_height.txt')
-wrf.write_inputfile('inputs')
+```
+
+The namelist conversion may also be accomplished from the command line:
+```bash
+wrf_namelist_to_erf namelist.input inputs
 ```
 
 ### Postprocessing data logs
-Data logs are output with the `erf.data_log` param and can include time histories of surface conditions and planar averaged profiles (e.g., for idealized LES simulations)
+Data logs are output with the `erf.data_log` param and can include time
+histories of surface conditions and planar averaged profiles (e.g., for
+idealized LES simulations)
 ```python
 from erftools.postprocessing import DataLog
 log = DataLog(f'{simdir}/surf_hist.dat',
@@ -31,6 +60,8 @@ print(log.ds) # data are stored in an xarray dataset
 
 Some notes and recommendations:
 
-* An aspirational goal is to contribute code that can be used as in the examples above, with clear, intuitive naming and following PEP-8 style as a set of guidelines rather than gospel.
-* To avoid duplication, model constants are defined in `erf.constants`, which should replicate `ERF/Source/ERF_Constants.H`.
-* In the same vein, equation of state evaluations are defined in `erf.EOS`, which should replicate `ERF/Source/Utils/ERF_EOS.H`. 
+* An aspirational goal is to contribute code that can be used as in the examples above, with clear, intuitive naming.
+* To avoid duplication, model constants are defined in `erftools.constants`, which should replicate `ERF/Source/ERF_Constants.H`.
+* In the same vein, equation of state evaluations are defined in `erftools.utils.EOS`, which should replicate `ERF/Source/Utils/ERF_EOS.H`.
+* Other utilities for calculating/deriving/diagnosing quantities of interest are also in `erftools.utils.*`
+* Please follow PEP-8 style--as a set of guidelines rather than gospel--to facilitate code usage and maintenance by the community.

erftools/parms.py → erftools/erfparms.py

@@ -16,12 +16,23 @@ def check_unknown_params(data_dict, dataclass_type):
     provided_params = set(data_dict.keys())
     unknown_params = list(provided_params - known_params)
 
+    # ignore MOST settings
+    unknown_params = [param for param in unknown_params
+                      if not param.startswith('most.')]
+
+    # ignore refinement indicators
     if 'refinement_indicators' in data_dict:
-        boxes = data_dict['refinement_indicators'].split()
-        for box in boxes:
+        refine_names = data_dict['refinement_indicators']
+        if isinstance(refine_names, str):
+            refine_names = [refine_names]
+        for box in refine_names:
             unknown_params = [param for param in unknown_params
                               if not param.startswith(f'{box}.')]
 
+    # ignore nc_init_file_*
+    unknown_params = [param for param in unknown_params
+                      if not param.startswith('nc_init_file_')]
+
     if unknown_params:
         warnings.warn(f'Non-standard {dataclass_type.__name__} ignored: {unknown_params}')
 
@@ -33,7 +44,6 @@ class AMRParms:
     max_level: int = 0
     ref_ratio: Union[int,List[int]] = 2
     ref_ratio_vect: List[int] = field(default_factory=list)
-    regrid_int: int = -1
 
     v: int = 0  # verbosity
 
@@ -125,6 +135,7 @@ class ERFParms:
 
     # Refinement
     refinement_indicators: Union[str,List[str]] = field(default_factory=list)
+    regrid_int: int = -1
 
     # Grid Stretching
     grid_stretching_ratio: float = 1.
@@ -240,9 +251,6 @@ class ERFParms:
     rayleigh_dampcoef: float = 0.2
     rayleigh_zdamp: float = 500.
 
-    # BCs
-    use_explicit_most: bool = False
-
     # Initialization
     init_type: str = 'None'
     init_sounding_ideal: bool = False
@@ -263,6 +271,27 @@ class ERFParms:
 
     # Radiation
     radiation_model: str = 'None'
+    rad_freq_in_steps: int = 1
+    rad_write_fluxes: bool = False
+    #nswbands: int = 14
+    #nlwbands: int = 16
+    #nswgpts: int = 224
+    #nlwgpts: int = 256
+    co2vmr: float = 0.00036 # from DevTests/Radiation
+    o3vmr: Union[float,List[float]] = 1.887e-7 # default
+    n2ovmr: float = 3.2e-7 # from DevTests/Radiation
+    covmr: float = 1.5e-7 # from DevTests/Radiation
+    ch4vmr: float = 1.7e-6 # from DevTests/Radiation
+    o2vmr: float = 0.209 # default
+    n2vmr: float = 0.7906 # default
+    rrtmgp_file_path: str = './'
+    #rrtmgp_coeffs_sw: str = 'rrtmgp-data-sw-g224-2018-12-04.nc'
+    #rrtmgp_coeffs_lw: str = 'rrtmgp-data-lw-g256-2018-12-04.nc'
+    #rrtmgp_cloud_optics_sw: str = 'rrtmgp-cloud-optics-coeffs-sw.nc'
+    #rrtmgp_cloud_optics_lw: str = 'rrtmgp-cloud-optics-coeffs-lw.nc'
+
+    # Land Surface
+    land_surface_model: str = 'None'
 
     def __post_init__(self):
         if self.anelastic:
@@ -321,15 +350,15 @@ def __post_init__(self):
 
         les_types = self.les_type if isinstance(self.les_type,list) else [self.les_type]
         for turbmodel in les_types:
-            assert turbmodel in ['None','Smagorinsky','Deardorff'], \
+            assert turbmodel in ['None','Smagorinsky','Smagorinsky2D','Deardorff'], \
                     f'Unexpected erf.les_type={turbmodel}'
         if any([turbmodel == 'Smagorinsky' for turbmodel in les_types]):
             smag_Cs = self.Cs if isinstance(self.Cs,list) else len(les_types)*[self.Cs]
             assert all([Cs > 0 for Cs in smag_Cs]), 'Need to specify valid Smagorinsky Cs'
 
         pbl_types = self.pbl_type if isinstance(self.pbl_type,list) else [self.pbl_type]
         for pblscheme in pbl_types:
-            assert pblscheme in ['None','MYNN25','YSU'], \
+            assert pblscheme in ['None','MYJ','MYNN25','MRF','YSU'], \
                     f'Unexpected erf.pbl_type={pblscheme}'
 
         assert self.abl_driver_type in \

erftools/inputs.py

@@ -1,10 +1,14 @@
 import sys
 import warnings
+from datetime import datetime
+from dateutil import parser
+
+import click
 import contextlib
 import numpy as np
 
 # parameters parsed by AMReX's ParmParse
-from .parms import AMRParms, GeometryParms, ERFParms, check_unknown_params
+from .erfparms import AMRParms, GeometryParms, ERFParms, check_unknown_params
 
 
 def parmparse(prefix, ppdata):
@@ -37,6 +41,14 @@ def list_to_str(mylist,dtype=None):
 def strs_to_str(mylist):
     return ' '.join([s.strip('"').strip("'") for s in mylist])
 
+def to_datetime(dt):
+    # ensure that we have a datetime object, parse if needed
+    if isinstance(dt, datetime):
+        return dt
+    elif isinstance(dt, list):
+        assert len(dt) == 2, 'expected date string + time string'
+        dt = ' '.join(dt)
+    return parser.parse(dt)
 
 class ERFInputs(object):
     """Input data container with validation and output"""
@@ -50,8 +62,13 @@ def __init__(self,inpfile=None,**ppdata):
         self.max_step = int(ppdata.get('max_step',-1))
         self.start_time = float(ppdata.get('start_time',0.))
         self.stop_time = float(ppdata.get('stop_time',1e34))
-        self.start_date = ppdata.get('start_date',None)
-        self.stop_date = ppdata.get('stop_date',None)
+
+        self.start_datetime = ppdata.get('start_datetime',None)
+        if self.start_datetime is not None:
+            self.start_datetime = to_datetime(self.start_datetime)
+        self.stop_datetime = ppdata.get('stop_datetime',None)
+        if self.stop_datetime is not None:
+            self.stop_datetime = to_datetime(self.stop_datetime)
 
         # read amr, geometry, and erf inputs
         amrparms = parmparse('amr',ppdata)
@@ -128,7 +145,7 @@ def read_bcs(self,ppdata):
             assert 'zhi.type' in ppdata.keys()
             self.zlo = parmparse('zlo',ppdata)
             self.zhi = parmparse('zhi',ppdata)
-            if self.zlo['type'] == 'MOST':
+            if self.zlo['type'] == 'surface_layer':
                 self.most = parmparse('erf.most',ppdata)
 
     def read_refinement(self,ppdata):
@@ -164,7 +181,9 @@ def read_refinement(self,ppdata):
                     self.refine[box][test] = thresh
 
     def validate(self):
-        # additional validation that depends on different parmparse types
+        """Additional validation (outside of what is performed in
+        erfparms.*) that depends on different parmparse tables
+        """
         if self.erf.terrain_type.lower() != 'none':
             if self.erf.terrain_z_levels:
                 nz = self.amr.n_cell[2]
@@ -178,19 +197,22 @@ def validate(self):
                               f' will be inactive with'
                               f' amr.max_level={self.amr.max_level}')
             if 'field_name' in refineparams and \
-                    self.amr.regrid_int <= 0:
+                    self.erf.regrid_int <= 0:
                 warnings.warn(f'{box} dynamic refinement will be inactive'
-                              f' with amr.regrid_int={self.amr.regrid_int}')
+                              f' with erf.regrid_int={self.erf.regrid_int}')
+
+        if isinstance(self.erf.o3vmr, list):
+            assert len(self.erf.o3vmr) == self.amr.n_cell[2]
 
     def write(self,fpath=None):
         with open_file_or_stdout(fpath) as f:
             f.write('# ------------------------------- INPUTS TO ERF -------------------------------\n')
             f.write('# written by erftools.inputs (https://github.com/erf-model/erftools)\n')
-            if self.start_date and self.stop_date:
+            if self.start_datetime and self.stop_datetime:
                 f.write(f"""
 max_step       = {self.max_step}
-start_datetime = {self.start_date}  # epoch time: {self.start_time} s
-stop_datetime  = {self.stop_date}  # epoch time: {self.stop_time} s
+start_datetime = {self.start_datetime}  # epoch time: {self.start_time} s
+stop_datetime  = {self.stop_datetime}  # epoch time: {self.stop_time} s
 """)
             else:
                 f.write(f"""
@@ -238,14 +260,15 @@ def write(self,fpath=None):
             if hasattr(self,'zlo'):
                 f.write('\n')
                 write_bc(f,'zlo',self.zlo)
-                if self.zlo['type'] == 'MOST':
+                if self.zlo['type'] == 'surface_layer':
                     for key,val in self.most.items():
                         f.write(f'erf.most.{key} = {val}\n')
                 write_bc(f,'zhi',self.zhi)
 
             ########################################
             f.write('\n# TIME STEP CONTROL\n')
-            f.write(f'erf.substepping_type = {self.erf.substepping_type}\n')
+            if self.erf.substepping_type != 'implicit':
+                f.write(f'erf.substepping_type = {self.erf.substepping_type}\n')
             if self.erf.fixed_dt > 0:
                 f.write(f'erf.fixed_dt         = {self.erf.fixed_dt}\n')
             else:
@@ -270,7 +293,7 @@ def write(self,fpath=None):
             ########################################
             f.write('\n# REFINEMENT / REGRIDDING\n')
             f.write(f'amr.max_level = {self.amr.max_level}\n')
-            f.write(f'amr.regrid_int = {self.amr.regrid_int}\n')
+            f.write(f'erf.regrid_int = {self.erf.regrid_int}\n')
             if len(self.erf.refinement_indicators) > 0:
                 if len(self.amr.ref_ratio_vect) > 0:
                     f.write('amr.ref_ratio_vect = '
@@ -337,14 +360,15 @@ def write(self,fpath=None):
             f.write('\n# PLOTFILES\n')
             if self.erf.plotfile_type != 'amrex':
                 f.write(f'erf.plotfile_type = {self.erf.plotfile_type}\n')
-            if (self.erf.plot_int_1 > 0) or (self.erf.plot_per_1 > 0):
-                f.write(f'erf.plot_file_1 = {self.erf.plot_file_1}\n')
-                if self.erf.plot_per_1 > 0:
-                    f.write(f'erf.plot_per_1  = {self.erf.plot_per_1}\n')
-                else:
-                    f.write(f'erf.plot_int_1  = {self.erf.plot_int_1}\n')
-                f.write('erf.plot_vars_1 = '
-                        f"{strs_to_str(self.erf.plot_vars_1)}\n")
+
+            f.write(f'erf.plot_file_1 = {self.erf.plot_file_1}\n')
+            if self.erf.plot_per_1 > 0:
+                f.write(f'erf.plot_per_1  = {self.erf.plot_per_1}\n')
+            else:
+                f.write(f'erf.plot_int_1  = {self.erf.plot_int_1}\n')
+            f.write('erf.plot_vars_1 = '
+                    f"{strs_to_str(self.erf.plot_vars_1)}\n")
+
             if (self.erf.plot_int_2 > 0) or (self.erf.plot_per_2 > 0):
                 f.write(f'erf.plot_file_2 = {self.erf.plot_file_2}\n')
                 if self.erf.plot_per_2 > 0:
@@ -404,13 +428,14 @@ def write(self,fpath=None):
             lestypes = self.erf.les_type \
                     if isinstance(self.erf.les_type,list) \
                     else [self.erf.les_type]
-            have_smag = ('Smagorinsky' in lestypes)
+            have_smag = ('Smagorinsky' in lestypes) or ('Smagorinsky2D' in lestypes)
             have_dear = ('Deardorff' in lestypes)
             if have_smag or have_dear:
                 f.write(f'\nerf.les_type  = {strs_to_str(lestypes)}\n')
             if have_smag:
                 Cs_list = self.erf.Cs if isinstance(self.erf.Cs,list) else [self.erf.Cs]
-                f.write(f'erf.Cs        = {list_to_str(Cs_list)}\n')
+                comment = ' # Cs=0.25 is WRF default' if (0.25 in Cs_list) else ''
+                f.write(f'erf.Cs        = {list_to_str(Cs_list)}{comment}\n')
             if have_dear:
                 Ck_list = self.erf.Ck if isinstance(self.erf.Ck,list) else [self.erf.Ck]
                 f.write(f'erf.Ck        = {list_to_str(Ck_list)}\n')
@@ -445,7 +470,28 @@ def write(self,fpath=None):
             ########################################
             if self.erf.radiation_model != 'None':
                 f.write(f"""\n# RADIATION
-erf.radiation_model = {self.erf.radiation_model}
+erf.radiation_model   = {self.erf.radiation_model}
+erf.rad_freq_in_steps = {self.erf.rad_freq_in_steps}
+erf.rad_write_fluxes  = {bool_to_str(self.erf.rad_write_fluxes)}
+erf.co2vmr            = {self.erf.co2vmr}
+""")
+                if isinstance(self.erf.o3vmr, (list,np.ndarray)):
+                    f.write(f'erf.o3vmr             = {list_to_str(self.erf.o3vmr)}')
+                else:
+                    f.write(f'erf.o3vmr             = {self.erf.o3vmr}')
+                f.write(f"""
+erf.n2ovmr            = {self.erf.n2ovmr}
+erf.covmr             = {self.erf.covmr}
+erf.ch4vmr            = {self.erf.ch4vmr}
+erf.o2vmr             = {self.erf.o2vmr}
+erf.n2vmr             = {self.erf.n2vmr}
+erf.rrtmgp_file_path  = {self.erf.rrtmgp_file_path}
+""")
+
+            ########################################
+            if self.erf.land_surface_model != 'None':
+                f.write(f"""\n# LAND SURFACE
+erf.land_surface_model = {self.erf.land_surface_model}
 """)
 
             ########################################
@@ -516,12 +562,15 @@ def write(self,fpath=None):
 erf.init_sounding_ideal = {bool_to_str(self.erf.init_sounding_ideal)}
 """)
             elif self.erf.init_type.lower() == 'wrfinput':
-                f.write(f"""
-erf.init_type      = wrfinput
-erf.use_real_bcs   = {bool_to_str(self.erf.use_real_bcs)}
-erf.nc_init_file_0 = {self.erf.nc_init_file_0}""")
+                f.write('\nerf.init_type      = wrfinput\n')
+                nc_init_files = [
+                    param for param in dir(self.erf)
+                    if param.startswith('nc_init_file_')]
+                for initfile in nc_init_files:
+                    f.write(f'erf.{initfile} = {getattr(self.erf,initfile)}\n')
                 if self.erf.use_real_bcs:
                     f.write(f"""
+erf.use_real_bcs   = {bool_to_str(self.erf.use_real_bcs)}
 erf.nc_bdy_file    = {self.erf.nc_bdy_file}
 erf.real_width     = {self.erf.real_width}
 erf.real_set_width = {self.erf.real_set_width}
@@ -557,3 +606,10 @@ def write_bc(f, bcname, bcdict):
         else:
             assert isinstance(val, list)
             f.write(f"{bcname}.{key} = {list_to_str(val,float)}\n")
+
+
+@click.command()
+@click.argument('input_file', type=click.Path(exists=True, readable=True))
+def check_erf_input(input_file):
+    """Validate ERF input file"""
+    ERFInputs(input_file)