Hi~ I want to ask the difference between helpers/wrf2icar.sh and helpers/wrf/wrf2icar.py, are they independent or what contact in them?

This is really an ongoing challenge, but I'd like to make a concerted effort to get all routines documented and create more detailed doc in readthedocs on the usage of ICAR. ALL are welcome to help with this.

when (some/any) of the vegtype_var, soiltype_var, soil_t_var are defined in the options file, but lsm=0 is selected, the model crashes on initialization in __domain_interface.domain_implementation_MOD_read_land_variables at objects/domain_obj.f90:1622.

Possible Solution

check if lsm /= 0 before allocating vars, this is currently not done.

Steps to Reproduce (for bugs)

1. set lsm=0
   2.define (some/any) of the vegtype_var, soiltype_var, soil_t_var

Your Environment

• Version of the code used: develop branch, but probably all.
• Operating System and version: cori,
• Compiler and version: gnu
• Other relevant information:

Expected Behavior

Current Behavior

Possible Solution

Steps to Reproduce (for bugs)

Your Environment

• Version of the code used:
• Operating System and version:
• Compiler and version:
• Other relevant information:

Hi , recently I'm trying to use wrfout files to run HICAR. I want to know if I want to input wrfout, firstly should I use helpers/wrf/wrf2icar.py to convert wrfout, secondly I put the file to HICAR ,but there are some questions.. so to input wrfout, what is the true steps? thank you for your help

The fixed_dz_advection flag in the namelist is confusing for all but the most experienced users. For certain (most?) applications it should be set to true to produce realistic behaviour, after Dylans modifications of the advection code. Therefore I suggest we hardcode the desired behaviour into either the advection or wind code and remove the flag from the namelist.

Expected Behavior

for anything but linear theory in the wind settings, fixed_dz_advection should be TRUE. (but @gutmann can confirm this).

Both the iterative wind solver and the wind rotation to the grid appear to have minor errors that change the windfield when run in parallel relative to the value when run in serial. Setting a larger halo size (3) fixes this, but with a very large computational cost.

Expected Behavior

wind fields are the same in serial and parallel

Current Behavior

difference between serial and parallel wind fields appear between processes.

Possible Solution

increase size of u and v grids to provide more information outside of immediate boundary, without increasing the halo being passed in all variables.

Steps to Reproduce (for bugs)

1. cafrun -np 36 icar
2. cafrun -np 1 icar
3. difference output files

Your Environment

• Version of the code used: latest develop (pre 2.1 release)
• Operating System and version: linux
• Compiler and version: gfortran 10.1 opencoarrays 2.9.2

Expected Behavior

No warnings produced when running aggregate_parallel_files.py

Current Behavior

The following warning is produced:

/glade/u/apps/opt/conda/envs/npl-2024a/lib/python3.11/site-packages/xarray/core/utils.py:494:
   FutureWarning: The return type of `Dataset.dims` will be changed to return a set of dimension 
   names in future, in order to be more consistent with `DataArray.dims`. To access a mapping 
   from dimension names to lengths, please use `Dataset.sizes`.
  warnings.warn(

Possible Solution

Switch to Dataset.sizes like the future warning recommends. It seems to be occurring from the following section, but the warning gets cut off so it might be a different place.

        if len(dims) == 1:
            nt = d.dims[dims[0]]
            data = np.zeros((nt))
        if len(dims) == 2:
            data = np.zeros((ny + y_off, nx + x_off))
        if len(dims) == 3:
            data = np.zeros((d.dims[dims[0]], ny + y_off, nx + x_off))
        if len(dims) == 4:
            nt = d.dims[dims[0]]
            nz = d.dims[dims[1]]
            data = np.zeros((nt, nz, ny + y_off, nx + x_off))

Steps to Reproduce (for bugs)

1. Run helpers/aggregate_parallel_files.py on ICAR output files.

Your Environment

• Version of the code: develop branch
• Operating System and version: Derecho
• Compiler and version: Python 3.11.6
• Other relevant information: Conda environment npl-20241

Convert IO code to loop through a list instead of explicitly adding code for every new variable. Requires additional data structures to track metadata for all variables, but not a huge lift.

Expected Behavior

gen_ideal_test.py should be restored to previous version to be truly a generalized test case. Or a new method that is more generalized and easier to produce tests with should be created.

Current Behavior

gen_ideal_test.py is an ideal test case for the Sleve work in #128

Possible Solution

gen_ideal_test.py should be moved to gen_sleve_ideal_test.py and the previous gen_ideal_test.py should be stored or a new idea test case that is more generalized and easier to produce tests with should be created.

Using the linear wind removal feature requires that the input data be on atmospheric model levels (or better, constant height above the terrain.

It would be nice if this worked when input forcing data were on pressure levels. This requires that either input data be initially interpolated to a common vertical structure, or that the look up table be generated for arbitrary heights. The latter is hard, particularly for points below topography, and time varying heights...

Hi Ethan

Compiling the latest icar-v2 (with fftw-3.3.8, libsz-2.1.1, hdf5-1.12.0, netcdf-c-4.7.4, netcdf-fortran-4.5.3), with the gfortran-10 compiler and CAF_MODE=single flag, I get the following error:

10 | submodule(domain_interface) domain_implementation
| 1
......
16 | use io_routines, only : io_read, io_write
| 2
Error: Symbol ‘io_read’ at (1) conflicts with the symbol at (2)

Could it be that io_read is already defined, so it creates this conflict here?
(1 ist at domain_interface, 2 at io_read)

Best wishes
Sebastian

Specify an option for Po to be used in converting potential temperature to temperature. Right now this is hard coded at 100000Pa, but if another dataset uses something else (e.g. 102000) this will result in the wrong temperature being computed internally.

The issue occurs when there are vertical ICAR levels at an elevation lower than the lowest vertical level of the forcing. E.g. the lowest forcing level is defined at 500m while there are vertical ICAR levels defined at 25m, 87.5m, 175m and 325m => four levels below the lowest forcing level.

ICAR does not extrapolate quantities linearly to these lower levels but rather keeps them constant. The following plot shows the extrapolation behavior of ICAR so far (ICAR old), the expected behavior (ICAR new) and the forcing data (forcing)

This particularly leads to artefacts when calculating the Brunt-Väisälä frequency. For ICAR old this results in a layer of low atmospheric stability in the lowest levels below the first vertical forcing level:

The proposed fix for this issue is to allow vertical extrapolation downward via linear interpolation.

Newly added convection schemes (BMJ, NSAS) require PBL height as an input. The current 'simple' PBL scheme does not compute this.
Currently the PBL height input to the above-mentioned convection schemes is fixed at an arbitrary value (hard-coded in cu_driver.f90).

First, we should see how sensitive these schemes actually are to this hard-coded PBL value by varying the value and analyzing the output. If this is negligible, the need for PBL improvement may be small. If it turns out that this is significant, then two possible ways forward exist:

1. include a more sophisticated scheme that calculates PBL height. YSU is a logical choice, since it is already partially implemented.
2. Expanding the existing scheme based on Hong and Pan 1996 could be expanded to calculate the PBL height.

Travis CI no longer builds since it migrated to travis.com.

Expected Behavior

Some sort of CI should run successfully on all PRs prior to merging.

Current Behavior

Nothing happens

Possible Solution

Setup github actions to run CI builds within github!

Steps to Reproduce (for bugs)

1. Create PR to github.com/NCAR/icar
2. There's no step 2

Many helper scripts are in disrepair, no longer work, or at least don't work with py3. Need someone to clean up, make them more generally useful, and "modern"

Apparently, double precision is not enough to store time accurately. "easy" switch: modify the time object (and maybe the time delta object though that isn't as necessary) to use quad-precision internally, while maintaining the same external interface. Also, make sure rounding to the nearest second for output files works correctly.

I've tried assigning the absolute and relative path for LUT_filename in the namelist file, but when read_LUT = True and write_LUT = False, ICAR still creates a linear theory LUT. I'm using the same hi res file (domain) every time. I also tried commenting write_LUT, but it still created the LUT.

I had problem with compiling ICAR with gfortran, and I have tried gfortran (GCC) of version 4.4.7 and version 4.9.4, but both give error either like:
f951: error: unrecognized command line option "-fno-protect-parens"
make: *** [build/icar_constants.o] Error 1

or like:
"collect2: error: ld returned 1 exit status"
make: *** [build/icar_constants.o] Error 1

I'm trying to run an ideal case to make sure ICAR is installed properly. After creating a python2.7 environment and installing pynio and pyngl, the script 'create_geo_testfiles.py' works (creates files named lo.nc and hi.nc). The gen_init_ideal.py and gen_ideal.py scripts fail due to the following error:

NameError: a global name '_nio_addvar' is not defined

It appears _nio_addvar is not present in the mygis.py module.

In order to make gen_ideal.py and gen_init_ideal.py work (until the problem described above) I had to change the dimension variables from floats to integers (nz, nx, ny, udims, vdims, and dims) because python would complain about using float values as indices.

Interpolating across a large gradient (e.g. 25->0) can produce small negative values at what should be the 0 edge. This may be a nearly unavoidable rounding/truncation error, but needs to be investigated as it can create issues for precipitation.

This is probably identical to NCAR/GARD#63 as they use the same regridding code.

question1: in Topography.py, in function "genHill": j = (np.arange(self.nx) - self.ny/2) / self.ny * np.pi *2, is this nx should be ny?

question2: in icar, if all variables's levels must be same? in wrfout data , vertical velocity w's levels are u's levels plus 1,
when running icar, there will be some errors.

question3: also in Topography.py, is Schar's advection test's code correct?, in Schar's article in 2004, I don't found the same formula ,especially "hgt = (
hill_height *
self.h_x(ig, lmbda) * self.h_x_star( ig, a)
)" , that make me very confused

hope your reply, thank you !

I'm tracking down whether this is a problem with the input data or the code itself. ICAR was built with MODE=debugslow on Intel compilers. About 30% of the way through the first time step an error forrtl: error(73): floating divide by zero occurs. The domain is situated over Idaho.

Physics Options:

pbl = 2 ! Simple
lsm = 3 ! Noah
water = 2 ! simple
mp = 2 ! simple
rad = 2 ! simple
conv = 1 ! Tiedke Scheme
adv = 1 ! upwind
wind = 1 ! linear theory

The traceback provides the following (edited) information:

Routine             Line               Source
module_cu_tiedtke    2923    cu_tiedtke.f90
module_cu_tiedtke    1211    cu_tiedtke.f90
module_cu_tiedtke     648     cu_tiedtke.f90
module_cu_tiedtke     443     cu_tiedtke.f90
convection_mp_con    153     cu_driver.f90
time_step_mp_step     477     time_step.f90
MAIN__                         90     driver.f90

ICAR is producing some weird sharp horizontal features that appear to be driving by similarly sharp linear features in the vertical wind field. This likely stems from an large jump in linear theory LUT values for the wind perturbation over a single step in the direction axis of the LUT. This step change occurs in both U and V fields and would result in large w values.

Problem found when running the 4km headwaters case with ERAi forcing for Sept. 2013.
LUT values [ndir=24, nspd=8. nsqr=8] from around dir=10 there is a step change, then around dir=20 it flips back again...

When using WRF output (27 km as low res, 9 km as high res) as forcing data for ICAR, I get the following error:

Initializing Boundaries
Interpreting geopotential height as residing between model layers
Setting up domain geographic Look Up Tables
ERROR finding triangle
39 39 38 38
61 60 61 60
33.12846 238.2691
33.35105 238.3338
33.14337 238.4846
33.22448 238.0857
33.01720 238.2369
ERROR: Point not located in bounding triangle
238.2691 33.12846
Triangle vertices
238.3338 33.35105
238.4846 33.14337
238.2852 33.18402
-0.3326202 2.0896307E-04 1.332411
Point in poly returns: F
Triangulation is broken

Permit advection algorithm to be specified separately for each advected scalar (e.g. upwind for a/b, MPDATA w/o FCT for c and MPDATA w/FCT for d...) and make Advection code loop through a list instead of specifying each variable by hand to permit arbitrary new species to be advected.

Add MPI or coarray fortran parallelization. This is a fairly large issue, but is tractable.

The 2.0 implementation of output variables uses objects, and could be extended to permit output variables to "compute" themselves prior to output, e.g. winds could rotate to east-west/N-S coordinates, accumulated_precip could add a bucket, time averaging or min/max statistics could be added.

It is not clear how some of this would work for variables that are output to multiple files (e.g. restart and standard output files), would separate variables need to be instantiated in each output object?

Hi~ I want to know how to modify the model top height, and in icar_options.nml, what's the usage of "dz_levels" ,could I modify the number of "dz_levels" and if could, according to what to modify, thank you for your reply~

Expected Behavior

If icar executable exits with STOP the github action action should fail.

Current Behavior

Looking at a recent github actions report and clicking on Execute Test Run, one can see that the action is reporting success when it should be a failure.

Section of Execute Test Run output in build_and_test_ICAR shows

 Initializing Options
STOP Options file does not exist. 
 Using options file = icar_options.nm

which reports

build_and_test_ICAR
succeeded 15 days ago in 12m 23

Possible Solution

Unsure, working on it now.

Steps to Reproduce (for bugs)

1. Looked at Github Actions report, unsure behavior when testing locally.

Expected Behavior

With iterative winds, vertical motion should smoothly decay to 0 at the top of the model

Current Behavior

Currently with a flat z height below the model top, vertical motion decays to 0 at that height, then everything above it often has a large magnitude.

Possible Solution

investigate wind solver and z coordinate calculation, update iterative winds to always use model top.

Steps to Reproduce (for bugs)

1. set flat_z_height to -5
2. set wind = 3
3. may need to use sleeve coordinate?
4. run a domain with some topography

Your Environment

• all environments

Expected Behavior

Builds correctly.

Current Behavior

We have gotten a report that when compiling with the new CCE 17 Fortran compiler, the following error is being reported:

ftn-1917 ftn: ERROR default$mapper$exchangeable_t$exchangeable_interface_,
                File = objects/exchangeable_h.f90, Line = 1, Column = 1

Since the type for base_ptr has a coarray ultimate component, it must be a
nonpointer nonallocatable scalar, must not be a coarray and must not be a
function result.

Steps to Reproduce (for bugs)

1. Load cce 17
2. Compile

Your Environment

• Version of the code used: develop branch
• Cluster: Perlmutter
• Compiler and version: Cray Compiler Environment 17

Permit the vertical coordinate system to vary spatially so that the atm grid can be smoothed over topography instead of strictly terrain following. This shouldn't be difficult but requires minor changes in many places, most notably in the advection code.

There have been issues reported by various people with ICAR not being able to "find" the correct grid cells when interpolating the forcing data to the ICAR grid. This started when moving to the stricter triangulation method instead of assuming a rectilinear forcing grid. It also seems more likely when using higher resolution ICAR grids.

Need to work on making this interpolation routine a little more robust, and maybe falling back to something simpler with a warning printed instead of stopping completely when it may just be a floating point precision / rounding issue.

When I run icar，There was a issue. I don't know how to solve it. It is like this:

Initializing Options
Using options file = short_icar_options.nml
Model version: 1.0.1

Initializing Domain

Initializing Boundaries
Interpreting geopotential height as residing between model layers
Setting up domain geographic Look Up Tables
ERROR: Failed to find point 97.3592224 24.5479889 in a quadrant near:
97.6912231 24.6180420
1 1
109 109

The helper scripts are doing some imports from python modules I can't find. The culprits are
stat_down
mygis
load_data

Is there a reason why the saturation vapor pressure formulas are different in the LSM driver (https://github.com/NCAR/icar/blob/master/physics/lsm_driver.f90#L110) and the the simple radiation code (https://github.com/NCAR/icar/blob/master/physics/ra_simple.f90#L88)?

If not, the ra_simple version could be updated for consistency.