Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

How to "borrow" (unused) NUDAPT variables in SLUCM for other purposes? #6

Open
sunt05 opened this issue Nov 16, 2023 · 10 comments
Open

Comments

@sunt05
Copy link
Collaborator

sunt05 commented Nov 16, 2023

The following is drafted by @DanLi-BU with light edits by @sunt05


OK, now at least I think I can elaborate my questions clearly. Here they are:

  1. in wrfinput_d04 (attached), there is a variable called BUILD_HEIGHT, which corresponds to HGT_URB2D in the WRF code according to this line in Registry.EM_COMMON:

state real HGT_URB2D ij misc 1 - ir "BUILD_HEIGHT" "AVERAGE BUILDING HEIGHT WEIGHTED BY BUILDING PLAN AREA" "m"

  1. HGT_URB2D/BUILD_HEIGHT is one of the key variables in NUDAPT. If HGT_URB2D is larger than 0, it will trigger using NUDAPT, see https://ral.ucar.edu/sites/default/files/public/product-tool/NUDAPT_44_Documentation.pdf
    This is also corroborated by looking at the function urban_var_init in module_sf_urban.F

  2. Here is the question. If I manually modified BUILD_HEIGHT in wrfinput_d04, say make it 10 m everywhere, and save it as BUILD_HEIGHT in wrfinput_d04. This has NO effect on my WRF results. Well, now that I think about it more, I start to understand why, because WRF actually never has a variable called BUILD_HEIGHT (or never reads in a variable called BUILD_HEIGHT from wrfinput_d0x). The variable name in the code is HGT_URB2D, even though BUILD_HEIGHT is the variable name in wrfinput_d04.

  3. Similarly (but also not exactly), there is a variable called STDH_URB2D in wrfinput_d04 (again see attached), which is also called STDH_URB2D in the WRF code.

state real STDH_URB2D ij misc 1 - ir "STDH_URB2D" "Standard Deviation of Building Height" "m2"

But when I modified it, there was again no effect on the WRF outputs.

  1. Then I searched STDH_URB2D in all .F files and here is what I got:
./WRF//phys/module_sf_urban.F:2832:                             HGT_URB2D,MH_URB2D,STDH_URB2D,                & ! inout
./WRF//phys/module_sf_urban.F:2928:   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: STDH_URB2D
./WRF//phys/module_sf_urban.F:2976:               IF (HGT_URB2D(I,J)>0. .or. STDH_URB2D(I,J)<0.) THEN
./WRF//phys/module_sf_urban.F:2981:                  WRITE(mesg,*) 'STDH_URB2D',STDH_URB2D(I,J)
./WRF//phys/module_sf_urban.F:2988:                     STDH_URB2D(I,J)=0.
./WRF//phys/module_sf_urban.F:3013:               IF (HGT_URB2D(I,J)>0. .or. STDH_URB2D(I,J)<0. ) THEN
./WRF//phys/module_sf_urban.F:3018:                  WRITE(mesg,*) 'STDH_URB2D',STDH_URB2D(I,J)
./WRF//phys/module_sf_urban.F:3025:                     STDH_URB2D(I,J)=0.
./WRF//phys/module_sf_urban.F:3050:              IF (HGT_URB2D(I,J)>0. .or. STDH_URB2D(I,J)<0. ) THEN
./WRF//phys/module_sf_urban.F:3055:                  WRITE(mesg,*) 'STDH_URB2D',STDH_URB2D(I,J)
./WRF//phys/module_sf_urban.F:3062:                     STDH_URB2D(I,J)=0.
./WRF//phys/module_sf_urban.F:3087:              IF (HGT_URB2D(I,J)>0. .or. STDH_URB2D(I,J)<0. ) THEN
./WRF//phys/module_sf_urban.F:3092:                  WRITE(mesg,*) 'STDH_URB2D',STDH_URB2D(I,J)
./WRF//phys/module_sf_urban.F:3099:                     STDH_URB2D(I,J)=0.
./WRF//phys/module_sf_urban.F:3130:                   STDH_URB2D(I,J)=0.
./WRF//phys/module_physics_init.F:156:                         HGT_URB2D,MH_URB2D,STDH_URB2D,          & !Optional multi-layer urban
./WRF//phys/module_physics_init.F:666:   REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: STDH_URB2D !SLUCM
./WRF//phys/module_physics_init.F:1419:                HGT_URB2D,MH_URB2D,STDH_URB2D,                  & !Optional multi-layer urban
./WRF//phys/module_physics_init.F:2429:                HGT_URB2D,MH_URB2D,STDH_URB2D,                  & !Optional multi-layer urban
./WRF//phys/module_physics_init.F:2778:    REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: STDH_URB2D !SLUCM
./WRF//phys/module_physics_init.F:3137:                              HGT_URB2D,MH_URB2D,STDH_URB2D,                   & !urban
./WRF//phys/module_physics_init.F:3254:                              HGT_URB2D,MH_URB2D,STDH_URB2D,                   & !urban
./WRF//phys/module_sf_noahmpdrv.F:2938:    STDH_URB = STDH_URB2D(I,J)
./WRF//dyn_em/module_initialize_real.F:3023:               grid%STDH_URB2D(i,j)  = grid%URB_PARAM(i,93,j)
./WRF//dyn_em/module_first_rk_step_part1.F:707:     &        ,MH_URB2D=grid%mh_urb2d,STDH_URB2D=grid%stdh_urb2d          & !SLUCM
./WRF//dyn_em/start_em.F:1137:                      grid%MH_URB2D,grid%STDH_URB2D,grid%LF_URB2D,                     & !SLUCM
  1. My understanding is that STDH_URB2D is created after run real.exe (the line of code highlighted earlier, also see here, but there was never a place where the WRF model "reads in" the variable STDH_URB2D from wrfinput_d04. So it makes sense that modifying it does not affect the simulation results.

  2. But then when there is indeed valid NUDAPT data (as we did in previous weeks), how does real.exe "pass" STDH_URB2D, as well as HGT_URB2D, to wrf.exe through wrfinput_d0x? Does this have anything to do with this grid%? Interesting...

@sunt05
Copy link
Collaborator Author

sunt05 commented Nov 16, 2023

The interesting exploration by @DanLi-BU does reminds me of the coupling work we did a while ago for WRF-SUEWS.

  1. But then when there is indeed valid NUDAPT data (as we did in previous weeks), how does real.exe "pass" STDH_URB2D, as well as HGT_URB2D, to wrf.exe through wrfinput_d0x? Does this have anything to do with this grid%? Interesting...

And it's right - the custom variables in wrfinput are indeed loaded via the grid%xx derived type - see here

So we will need to tweak the related part in module_first_rk_step_part1.F to get this working (hopefully...).

@sunt05 sunt05 changed the title How to *borrow* (unused) NUDAPT variables in SLUCM for other purposes? How to _borrow_ (unused) NUDAPT variables in SLUCM for other purposes? Nov 16, 2023
@sunt05 sunt05 changed the title How to _borrow_ (unused) NUDAPT variables in SLUCM for other purposes? How to "borrow" (unused) NUDAPT variables in SLUCM for other purposes? Nov 16, 2023
@DanLi-BU
Copy link
Owner

This is the modified wrfinput_d04. I changed STDH_URB2D but there are no effects.

wrfinput_d04_mod.nc.zip

@DanLi-BU

This comment was marked as resolved.

@sunt05
Copy link
Collaborator Author

sunt05 commented Nov 16, 2023

Please upload all input files so I can run this case at my side.

@DanLi-BU

This comment was marked as resolved.

@DanLi-BU

This comment was marked as resolved.

@DanLi-BU
Copy link
Owner

DanLi-BU commented Dec 7, 2023

This is the code I used to modify the wrfinputs

wrfinput_d04_modification.zip

@sunt05
Copy link
Collaborator Author

sunt05 commented Dec 7, 2023

Now we've resolved an issue with the STDH_URB2D variable, which lost its metadata during a Python-based modification process - just wanted to write down the interesting debugging journey that might help fellow WRF users/coders.


🔍 Initial Search & Dead Ends

  • Started by searching for where STDH_URB2D was loaded from within the WRF input files.
  • Realised variables like this aren't always in plain FORTRAN source code; they could be injected via Perl scripts.
  • Checked all occurrences of STDH_URB2D, but no clear reason why it showed zeros instead of expected values after running WRF.

💡 Breakthrough with Compilation Files

  • Turned our attention to intermediate .f90 files generated during compilation.

  • Used VSCode and its Fortran plugin, which made locating references easier.

  • Found key insights in share/mediation_wrfmain.F:

    CALL open_r_dataset ( fid, TRIM(inpname) , grid , config_flags , "DATASET=INPUT", ierr )

  • Further details were hidden behind macros here:

    #include "fine_stream_input.inc"

🛠️ Debugging Steps Taken

Implemented some debug lines using the existing code structure:

WRITE ( message , FMT = '("processing wrfinput file (stream 0) for domain ",I8)' ) grid%id
CALL wrf_debug              ( 100 , message )

Remember to adjust the debug_level in your namelist.input!

Added checks right after reading the NetCDF file revealed unexpected all-zero values for grid%stdh_urb2d.

🕵️‍♂️ Comparing Variables
Checked another urban-related variable, frc_urb2d, which loaded correctly – indicating an isolated issue with just one variable.

Solution Identified
A few lines of quick xarray code pinpointed that while modifying values, we didn't preserve attributes (attrs).
The fix? Modify only the necessary variable values rather than replacing entire objects.

@sunt05
Copy link
Collaborator Author

sunt05 commented Dec 7, 2023

This is the code I used to modify the wrfinputs

wrfinput_d04_modification.zip

Below is the core part of the zipped ipynb:

Click to expand!

This is the content that will be hidden until the section is expanded. You can include text, images, code, etc.

import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np


# List of file names to modify
file_names = ['d01', 'd02', 'd03', 'd04']

# Base path for the files
base_path = "wrfinput_"

# Iterate over each file, modify, and save
for file_name in file_names:
    # Open the dataset
    ds_path = base_path + file_name
    ds = xr.open_dataset(ds_path)

    # Your existing code for modifying stdh_urb2d and hgt_urb2d
    stdh_urb2d = ds['STDH_URB2D']
    hgt_urb2d = ds['MH_URB2D']
    frc_urb2d = ds['FRC_URB2D']
    lu_index = ds['LU_INDEX']

    # Define the condition for LU_INDEX
    lu_condition = lu_index.isin([23, 24, 25, 26])

    # Update values based on the conditions
    stdh_urb2d1 = xr.where((frc_urb2d > 0) & lu_condition, -4.5, stdh_urb2d)
    hgt_urb2d = xr.where((frc_urb2d > 0) & lu_condition, 15, hgt_urb2d)

    # Assign the modified variables back to the dataset
    ds['STDH_URB2D'] = stdh_urb2d1
    ds['MH_URB2D'] = hgt_urb2d

    # Save the modified dataset to a new file
    modified_ds_path = base_path + file_name + "_mod.nc"
    ds.to_netcdf(modified_ds_path)
    print(f"Saved modified dataset to {modified_ds_path}")

Updated version:

Click to expand!
# List of file names to modify
file_names = ["d01", "d02", "d03", "d04"]

# Base path for the files
base_path = "wrfinput_"

# Iterate over each file, modify, and save
for file_name in file_names:
    # Open the dataset
    ds_path = base_path + file_name
    ds = xr.open_dataset(ds_path)

    # Auxiliary variables
    frc_urb2d = ds["FRC_URB2D"]
    lu_index = ds["LU_INDEX"]

    # Define the dictionary of variables with their new values
    dict_var_mod = {
        "STDH_URB2D": -4.5,
        "MH_URB2D": 15,
    }

    # Define the condition for LU_INDEX
    lu_condition = lu_index.isin([23, 24, 25, 26])

    # Update values based on the conditions
    for var, val in dict_var_mod.items():
        ds[var] = ds[var].where((frc_urb2d > 0) & lu_condition, val)

    # Save the modified dataset to a new file
    modified_ds_path = base_path + file_name + "_mod.nc"
    ds.to_netcdf(modified_ds_path)
    print(f"Saved modified dataset to {modified_ds_path}")

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

When branches are created from issues, their pull requests are automatically linked.

2 participants