Adding spatially explicit albedo for SLUCM

[Peter9192]

In this project we want to make it possible to pass spatially explicit fields for albedo and emissivity into the model.

This issue is an inventory of the steps needed to make that work.

• Add new variables to input data in WPS. See e.g. https://www.youtube.com/watch?v=xibGSLSaj6w. Basically we need to add a new entry in the geogrid table, and make sure corresponding the source data is available in the static data folder. As a shortcut for initial testing, we could manually add a new variable to the geo_em or wrfinput files after running geogrid/metgrid.
• Add variables to registry. WRF automatically generates code to initialize all model state variables, etc. The state variables for the SLUCM are registered here:

  WRF/Registry/Registry.EM_COMMON

  Lines 819 to 835 in 0a11865

  	# urban model variables
  	state real DZR l em - Z r "DZR" "THICKNESSES OF ROOF LAYERS" "m"
  	state real DZB l em - Z r "DZB" "THICKNESSES OF WALL LAYERS" "m"
  	state real DZG l em - Z r "DZG" "THICKNESSES OF ROAD LAYERS" "m"
  	state real URB_PARAM i{urb}j misc 1 - i1 "URB_PARAM" "NUDAPT_NBSD Urban Parameters" "parameter"
  	state real LP_URB2D ij misc 1 - i01r "BUILD_AREA_FRACTION" "BUILDING PLAN AREA DENSITY" "dimensionless"
  	state real HI_URB2D i{uhi}j misc 1 Z ir "HEIGHT_HISTOGRAMS" "DISTRIBUTION OF BUILDING HEIGHTS" "dimensionless"
  	state real LB_URB2D ij misc 1 - ir "BUILD_SURF_RATIO" "BUILDING SURFACE AREA TO PLAN AREA RATIO" "dimensionless"
  	state real HGT_URB2D ij misc 1 - ir "BUILD_HEIGHT" "AVERAGE BUILDING HEIGHT WEIGHTED BY BUILDING PLAN AREA" "m"
  	state real MH_URB2D ij misc 1 - i01r "MH_URB2D" "Mean Building Height" "m"
  	state real STDH_URB2D ij misc 1 - ir "STDH_URB2D" "Standard Deviation of Building Height" "m2"
  	state real LF_URB2D i{udr}j misc 1 Z ir "LF_URB2D" "Frontal Area Index" "dimensionless"
  	state real ZD_URB2D ij misc 1 - i1 "ZD_URB2D" "Zero-plane Displacement" "m"
  	state real Z0_URB2D ij misc 1 - i01r "Z0_URB2D" "Roughness length for momentum" "m"
  	state real LF_URB2D_S ij misc 1 - i01r "LF_URB2D_S" "Frontal area index (no wind directional dependency)" ""
  	# AHE with month and hour dimension flattened to one dimension, Jan = (0:23), Feb = (24:47)
  	state real AHE i{m_hr}j misc 1 - i01r "AHE" "Anthropogenic heat emission" "W m-2"

  Evidently, albedo and emissivity are missing, so they need to be added. See https://www2.mmm.ucar.edu/wrf/users/wrf_users_guide/build/html/software.html#registry. If I understand correctly, if we set this to be an 'input' variable, it will automatically read it from the wrfinput.
• Use the new state variables in the code. The urban canopy model code lives here, and e.g. the albedo is used here:

  WRF/phys/module_sf_urban.F

  Lines 914 to 919 in 0a11865

  	SR1=SX*(1.-ALBR)
  	SGR1=SX*(1.-ALBV)
  	SG1=SX*VFGS*(1.-ALBG)
  	SB1=SX*VFWS*(1.-ALBB)
  	SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
  	SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)

  
  NB variables are defined here:

  WRF/phys/module_sf_urban.F

  Lines 233 to 238 in 0a11865

  	! ALBR [-] : surface albedo of roof
  	! ALBB [-] : surface albedo of building wall
  	! ALBG [-] : surface albedo of road
  	! EPSR [-] : surface emissivity of roof
  	! EPSB [-] : surface emissivity of building wall
  	! EPSG [-] : surface emissivity of road

  and are currently initialized as fixed values:
  

  WRF/phys/module_sf_urban.F

  Lines 486 to 487 in 0a11865

  	REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG
  	REAL :: EPSR, EPSB, EPSG, Z0R, Z0B, Z0G, Z0HB, Z0HG

  
  The URBPARAM table is read in subroutine urban_param_init:

  WRF/phys/module_sf_urban.F

  Lines 2099 to 2473 in 0a11865

  	if(USE_WUDAPT_LCZ.eq.0)then !AndreaLCZ
  	OPEN (UNIT=11, &
  	FILE='URBPARM.TBL', &
  	ACCESS='SEQUENTIAL', &
  	STATUS='OLD', &
  	ACTION='READ', &
  	POSITION='REWIND', &
  	IOSTAT=IOSTATUS)
  	IF (IOSTATUS > 0) THEN
  	FATAL_ERROR('Error opening URBPARM.TBL. Make sure URBPARM.TBL (found in run/) is linked to the running directory.')
  	ENDIF
  	else
  	OPEN (UNIT=11, &
  	FILE='URBPARM_LCZ.TBL', &
  	ACCESS='SEQUENTIAL', &
  	STATUS='OLD', &
  	ACTION='READ', &
  	POSITION='REWIND', &
  	IOSTAT=IOSTATUS)
  	IF (IOSTATUS > 0) THEN
  	FATAL_ERROR('Error opening URBPARM_LCZ.TBL. Make sure URBPARM_LCZ.TBL (found in run/) is linked to the running directory.')
  	ENDIF
  	endif
  	READLOOP : do
  	read(11,'(A512)', iostat=iostatus) string
  	if (iostatus /= 0) exit READLOOP
  	if (string(1:1) == "#") cycle READLOOP
  	if (trim(string) == "") cycle READLOOP
  	indx = index(string,":")
  	if (indx <= 0) cycle READLOOP
  	name = trim(adjustl(string(1:indx-1)))
  	! Here are the variables we expect to be defined in the URBPARM.TBL:
  	if (name == "Number of urban categories") then
  	read(string(indx+1:),*) icate
  	IF (.not. ALLOCATED(ZR_TBL)) then
  	ALLOCATE( ZR_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ZR_TBL in urban_param_init')
  	ALLOCATE( SIGMA_ZED_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating SIGMA_ZED_TBL in urban_param_init')
  	ALLOCATE( Z0C_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0C_TBL in urban_param_init')
  	ALLOCATE( Z0HC_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HC_TBL in urban_param_init')
  	ALLOCATE( ZDC_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ZDC_TBL in urban_param_init')
  	ALLOCATE( SVF_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating SVF_TBL in urban_param_init')
  	ALLOCATE( R_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating R_TBL in urban_param_init')
  	ALLOCATE( RW_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating RW_TBL in urban_param_init')
  	ALLOCATE( HGT_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating HGT_TBL in urban_param_init')
  	ALLOCATE( AH_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating AH_TBL in urban_param_init')
  	ALLOCATE( ALH_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ALH_TBL in urban_param_init')
  	ALLOCATE( BETR_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating BETR_TBL in urban_param_init')
  	ALLOCATE( BETB_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating BETB_TBL in urban_param_init')
  	ALLOCATE( BETG_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating BETG_TBL in urban_param_init')
  	ALLOCATE( CAPR_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPR_TBL in urban_param_init')
  	ALLOCATE( CAPB_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPB_TBL in urban_param_init')
  	ALLOCATE( CAPG_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPG_TBL in urban_param_init')
  	ALLOCATE( AKSR_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSR_TBL in urban_param_init')
  	ALLOCATE( AKSB_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSB_TBL in urban_param_init')
  	ALLOCATE( AKSG_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSG_TBL in urban_param_init')
  	ALLOCATE( ALBR_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBR_TBL in urban_param_init')
  	ALLOCATE( ALBB_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBB_TBL in urban_param_init')
  	ALLOCATE( ALBG_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBG_TBL in urban_param_init')
  	ALLOCATE( EPSR_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSR_TBL in urban_param_init')
  	ALLOCATE( EPSB_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSB_TBL in urban_param_init')
  	ALLOCATE( EPSG_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSG_TBL in urban_param_init')
  	ALLOCATE( Z0R_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0R_TBL in urban_param_init')
  	ALLOCATE( Z0B_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0B_TBL in urban_param_init')
  	ALLOCATE( Z0G_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0G_TBL in urban_param_init')
  	ALLOCATE( AKANDA_URBAN_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating AKANDA_URBAN_TBL in urban_param_init')
  	ALLOCATE( Z0HB_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HB_TBL in urban_param_init')
  	ALLOCATE( Z0HG_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HG_TBL in urban_param_init')
  	ALLOCATE( TRLEND_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating TRLEND_TBL in urban_param_init')
  	ALLOCATE( TBLEND_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating TBLEND_TBL in urban_param_init')
  	ALLOCATE( TGLEND_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating TGLEND_TBL in urban_param_init')
  	ALLOCATE( FRC_URB_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating FRC_URB_TBL in urban_param_init')
  	! ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
  	! if(allocate_status /= 0) FATAL_ERROR('Error allocating ROOF_WIDTH in urban_param_init')
  	! ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
  	! if(allocate_status /= 0) FATAL_ERROR('Error allocating ROAD_WIDTH in urban_param_init')
  	!for BEP
  	ALLOCATE( NUMDIR_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating NUMDIR_TBL in urban_param_init')
  	ALLOCATE( STREET_DIRECTION_TBL(MAXDIRS , ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating STREET_DIRECTION_TBL in urban_param_init')
  	ALLOCATE( STREET_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating STREET_WIDTH_TBL in urban_param_init')
  	ALLOCATE( BUILDING_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating BUILDING_WIDTH_TBL in urban_param_init')
  	ALLOCATE( NUMHGT_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating NUMHGT_TBL in urban_param_init')
  	ALLOCATE( HEIGHT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating HEIGHT_BIN_TBL in urban_param_init')
  	ALLOCATE( HPERCENT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating HPERCENT_BIN_TBL in urban_param_init')
  	ALLOCATE( COP_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating COP_TBL in urban_param_init')
  	ALLOCATE( BLDAC_FRC_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating BLDAC_FRC_TBL in urban_param_init')
  	ALLOCATE( COOLED_FRC_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating COOLED_FRC_TBL in urban_param_init')
  	ALLOCATE( PWIN_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating PWIN_TBL in urban_param_init')
  	ALLOCATE( BETA_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating BETA_TBL in urban_param_init')
  	ALLOCATE( SW_COND_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating SW_COND_TBL in urban_param_init')
  	ALLOCATE( TIME_ON_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating TIME_ON_TBL in urban_param_init')
  	ALLOCATE( TIME_OFF_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating TIME_OFF_TBL in urban_param_init')
  	ALLOCATE( TARGTEMP_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating TARGTEMP_TBL in urban_param_init')
  	ALLOCATE( GAPTEMP_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating GAPTEMP_TBL in urban_param_init')
  	ALLOCATE( TARGHUM_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating TARGHUM_TBL in urban_param_init')
  	ALLOCATE( GAPHUM_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating GAPHUM_TBL in urban_param_init')
  	ALLOCATE( PERFLO_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating PERFLO_TBL in urban_param_init')
  	ALLOCATE( HSESF_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating HSESF_TBL in urban_param_init')
  	ALLOCATE( PV_FRAC_ROOF_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating PV_FRAC_ROOF_TBL in urban_param_init')
  	ALLOCATE( GR_FRAC_ROOF_TBL(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating GR_FRAC_ROOF_TBL in urban_param_init')
  	endif
  	numdir_tbl = 0
  	street_direction_tbl = -1.E36
  	street_width_tbl = 0
  	building_width_tbl = 0
  	numhgt_tbl = 0
  	height_bin_tbl = -1.E36
  	hpercent_bin_tbl = -1.E36
  	!end BEP
  	else if (name == "ZR") then
  	read(string(indx+1:),*) zr_tbl(1:icate)
  	else if (name == "SIGMA_ZED") then
  	read(string(indx+1:),*) sigma_zed_tbl(1:icate)
  	else if (name == "ROOF_WIDTH") then
  	ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ROOF_WIDTH in urban_param_init')
  	read(string(indx+1:),*) roof_width(1:icate)
  	else if (name == "ROAD_WIDTH") then
  	ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
  	if(allocate_status /= 0) FATAL_ERROR('Error allocating ROAD_WIDTH in urban_param_init')
  	read(string(indx+1:),*) road_width(1:icate)
  	else if (name == "AH") then
  	read(string(indx+1:),*) ah_tbl(1:icate)
  	else if (name == "ALH") then
  	read(string(indx+1:),*) alh_tbl(1:icate)
  	else if (name == "FRC_URB") then
  	read(string(indx+1:),*) frc_urb_tbl(1:icate)
  	else if (name == "CAPR") then
  	read(string(indx+1:),*) capr_tbl(1:icate)
  	! Convert CAPR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
  	capr_tbl = capr_tbl * ( 1.0 / 4.1868 ) * 1.E-6
  	else if (name == "CAPB") then
  	read(string(indx+1:),*) capb_tbl(1:icate)
  	! Convert CABR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
  	capb_tbl = capb_tbl * ( 1.0 / 4.1868 ) * 1.E-6
  	else if (name == "CAPG") then
  	read(string(indx+1:),*) capg_tbl(1:icate)
  	! Convert CABG_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
  	capg_tbl = capg_tbl * ( 1.0 / 4.1868 ) * 1.E-6
  	else if (name == "AKSR") then
  	read(string(indx+1:),*) aksr_tbl(1:icate)
  	! Convert AKSR_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
  	AKSR_TBL = AKSR_TBL * ( 1.0 / 4.1868 ) * 1.E-2
  	else if (name == "AKSB") then
  	read(string(indx+1:),*) aksb_tbl(1:icate)
  	! Convert AKSB_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
  	AKSB_TBL = AKSB_TBL * ( 1.0 / 4.1868 ) * 1.E-2
  	else if (name == "AKSG") then
  	read(string(indx+1:),*) aksg_tbl(1:icate)
  	! Convert AKSG_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
  	AKSG_TBL = AKSG_TBL * ( 1.0 / 4.1868 ) * 1.E-2
  	else if (name == "ALBR") then
  	read(string(indx+1:),*) albr_tbl(1:icate)
  	else if (name == "ALBB") then
  	read(string(indx+1:),*) albb_tbl(1:icate)
  	else if (name == "ALBG") then
  	read(string(indx+1:),*) albg_tbl(1:icate)
  	else if (name == "EPSR") then
  	read(string(indx+1:),*) epsr_tbl(1:icate)
  	else if (name == "EPSB") then
  	read(string(indx+1:),*) epsb_tbl(1:icate)
  	else if (name == "EPSG") then
  	read(string(indx+1:),*) epsg_tbl(1:icate)
  	else if (name == "AKANDA_URBAN") then
  	read(string(indx+1:),*) akanda_urban_tbl(1:icate)
  	else if (name == "Z0B") then
  	read(string(indx+1:),*) z0b_tbl(1:icate)
  	else if (name == "Z0G") then
  	read(string(indx+1:),*) z0g_tbl(1:icate)
  	else if (name == "DDZR") then
  	read(string(indx+1:),*) dzr(1:num_roof_layers)
  	! Convert thicknesses from m to cm
  	dzr = dzr * 100.0
  	else if (name == "DDZB") then
  	read(string(indx+1:),*) dzb(1:num_wall_layers)
  	! Convert thicknesses from m to cm
  	dzb = dzb * 100.0
  	else if (name == "DDZG") then
  	read(string(indx+1:),*) dzg(1:num_road_layers)
  	! Convert thicknesses from m to cm
  	dzg = dzg * 100.0
  	else if (name == "BOUNDR") then
  	read(string(indx+1:),*) boundr_data
  	else if (name == "BOUNDB") then
  	read(string(indx+1:),*) boundb_data
  	else if (name == "BOUNDG") then
  	read(string(indx+1:),*) boundg_data
  	else if (name == "TRLEND") then
  	read(string(indx+1:),*) trlend_tbl(1:icate)
  	else if (name == "TBLEND") then
  	read(string(indx+1:),*) tblend_tbl(1:icate)
  	else if (name == "TGLEND") then
  	read(string(indx+1:),*) tglend_tbl(1:icate)
  	else if (name == "CH_SCHEME") then
  	read(string(indx+1:),*) ch_scheme_data
  	else if (name == "TS_SCHEME") then
  	read(string(indx+1:),*) ts_scheme_data
  	else if (name == "AHOPTION") then
  	read(string(indx+1:),*) ahoption
  	else if (name == "AHDIUPRF") then
  	read(string(indx+1:),*) ahdiuprf(1:24)
  	else if (name == "ALHOPTION") then
  	read(string(indx+1:),*) alhoption
  	else if (name == "ALHSEASON") then
  	read(string(indx+1:),*) alhseason(1:4)
  	else if (name == "ALHDIUPRF") then
  	read(string(indx+1:),*) alhdiuprf(1:48)
  	else if (name == "PORIMP") then
  	read(string(indx+1:),*) porimp(1:3)
  	else if (name == "DENGIMP") then
  	read(string(indx+1:),*) dengimp(1:3)
  	else if (name == "IMP_SCHEME") then
  	read(string(indx+1:),*) imp_scheme
  	else if (name == "IRI_SCHEME") then
  	read(string(indx+1:),*) iri_scheme
  	else if (name == "OASIS") then
  	read(string(indx+1:),*) oasis
  	else if (name == "GROPTION") then
  	read(string(indx+1:),*) groption
  	else if (name == "FGR") then
  	read(string(indx+1:),*) fgr
  	else if (name == "DZGR") then
  	read(string(indx+1:),*) dzgr(1:4)
  	!for BEP
  	else if (name == "STREET PARAMETERS") then
  	STREETLOOP : do
  	read(11,'(A512)', iostat=iostatus) string
  	if (string(1:1) == "#") cycle STREETLOOP
  	if (trim(string) == "") cycle STREETLOOP
  	if (string == "END STREET PARAMETERS") exit STREETLOOP
  	read(string, *) k ! , dirst, ws, bs
  	numdir_tbl(k) = numdir_tbl(k) + 1
  	read(string, *) k, street_direction_tbl(numdir_tbl(k),k), &
  	street_width_tbl(numdir_tbl(k),k), &
  	building_width_tbl(numdir_tbl(k),k)
  	enddo STREETLOOP
  	else if (name == "BUILDING HEIGHTS") then
  	read(string(indx+1:),*) k
  	HEIGHTLOOP : do
  	read(11,'(A512)', iostat=iostatus) string
  	if (string(1:1) == "#") cycle HEIGHTLOOP
  	if (trim(string) == "") cycle HEIGHTLOOP
  	if (string == "END BUILDING HEIGHTS") exit HEIGHTLOOP
  	read(string,*) dummy_hgt, dummy_pct
  	numhgt_tbl(k) = numhgt_tbl(k) + 1
  	height_bin_tbl(numhgt_tbl(k), k) = dummy_hgt
  	hpercent_bin_tbl(numhgt_tbl(k),k) = dummy_pct
  	enddo HEIGHTLOOP
  	pctsum = sum ( hpercent_bin_tbl(:,k) , mask=(hpercent_bin_tbl(:,k)>-1.E25 ) )
  	if ( pctsum /= 100.) then
  	write (*,'(//,"Building height percentages for category ", I2, " must sum to 100.0")') k
  	write (*,'("Currently, they sum to ", F6.2,/)') pctsum
  	FATAL_ERROR('pctsum is not equal to 100.')
  	endif
  	else if ( name == "Z0R") then
  	read(string(indx+1:),*) Z0R_tbl(1:icate)
  	else if ( name == "COP") then
  	read(string(indx+1:),*) cop_tbl(1:icate)
  	else if ( name == "BLDAC_FRC") then
  	read(string(indx+1:),*) bldac_frc_tbl(1:icate)
  	else if ( name == "COOLED_FRC") then
  	read(string(indx+1:),*) cooled_frc_tbl(1:icate)
  	else if ( name == "PWIN") then
  	read(string(indx+1:),*) pwin_tbl(1:icate)
  	else if ( name == "BETA") then
  	read(string(indx+1:),*) beta_tbl(1:icate)
  	else if ( name == "SW_COND") then
  	read(string(indx+1:),*) sw_cond_tbl(1:icate)
  	else if ( name == "TIME_ON") then
  	read(string(indx+1:),*) time_on_tbl(1:icate)
  	else if ( name == "TIME_OFF") then
  	read(string(indx+1:),*) time_off_tbl(1:icate)
  	else if ( name == "TARGTEMP") then
  	read(string(indx+1:),*) targtemp_tbl(1:icate)
  	else if ( name == "GAPTEMP") then
  	read(string(indx+1:),*) gaptemp_tbl(1:icate)
  	else if ( name == "TARGHUM") then
  	read(string(indx+1:),*) targhum_tbl(1:icate)
  	else if ( name == "GAPHUM") then
  	read(string(indx+1:),*) gaphum_tbl(1:icate)
  	else if ( name == "PERFLO") then
  	read(string(indx+1:),*) perflo_tbl(1:icate)
  	else if (name == "HSEQUIP") then
  	read(string(indx+1:),*) hsequip_tbl(1:24)
  	else if (name == "HSEQUIP_SCALE_FACTOR") then
  	read(string(indx+1:),*) hsesf_tbl(1:icate)
  	else if (name == "IRHO") then
  	read(string(indx+1:),*) IRHO_TBL(1:24)
  	else if ( name == "PV_FRAC_ROOF") then
  	read(string(indx+1:),*) pv_frac_roof_tbl(1:icate)
  	else if ( name == "GR_FRAC_ROOF") then
  	read(string(indx+1:),*) gr_frac_roof_tbl(1:icate)
  	else if (name == "GR_FLAG") then
  	read(string(indx+1:),*) gr_flag_tbl
  	else if (name == "GR_TYPE") then
  	read(string(indx+1:),*) gr_type_tbl
  	!end BEP
  	else
  	FATAL_ERROR('URBPARM.TBL: Unrecognized NAME = "'//trim(name)//'" in Subr URBAN_PARAM_INIT')
  	endif
  	enddo READLOOP
  	CLOSE(11)

  
  And the values for each grid point are then inferred in subroutine read_param:

  WRF/phys/module_sf_urban.F

  Lines 1998 to 2000 in 0a11865

  	ALBR= ALBR_TBL(UTYPE)
  	ALBB= ALBB_TBL(UTYPE)
  	ALBG= ALBG_TBL(UTYPE)

  
  To change this, I think the easiest option is to change the relevant lines (e.g. 914-919) to use a different variable, which could be added as one of the input (and output!?) variables to the model, similar to e.g. nudapt params here:

  WRF/phys/module_sf_urban.F

  Lines 392 to 402 in 0a11865

  	!-------------------------------------------------------------------------------
  	! I: NUDAPT Input Parameters
  	!-------------------------------------------------------------------------------
  	REAL, INTENT(INOUT) :: mh_urb ! mean building height [m]
  	REAL, INTENT(INOUT) :: stdh_urb ! standard deviation of building height [m]
  	REAL, INTENT(INOUT) :: hgt_urb ! area weighted mean building height [m]
  	REAL, INTENT(INOUT) :: lp_urb ! plan area fraction [-]
  	REAL, INTENT(INOUT) :: frc_urb ! urban fraction [-]
  	REAL, INTENT(INOUT) :: lb_urb ! building surface to plan area ratio [-]
  	REAL, INTENT(INOUT), DIMENSION(4) :: lf_urb ! frontal area index [-]
  	REAL, INTENT(INOUT) :: zo_check ! check for printing ZOC

  
  These variables then also need to be added as input to the subroutine, not sure whether they are passed directly to the urban subroutine, similar to here:

  WRF/phys/module_sf_urban.F

  Line 319 in 0a11865

  U10,V10,TH2,Q2,UST,mh_urb,stdh_urb,lf_urb, & ! O

  or to urban_var_init as in here

  WRF/phys/module_sf_urban.F

  Line 2592 in 0a11865

  HGT_URB2D,MH_URB2D,STDH_URB2D, & ! inout

  
  They also need to be passed in where the corresponding subroutines are called, e.g. here

  WRF/phys/module_physics_init.F

  Line 3331 in 0a11865

  HGT_URB2D,MH_URB2D,STDH_URB2D, & !urban

so far so good. I'm gonna see how far I get with this.

TODO: compare Jisk/wudapt/nudapt approaches

[Peter9192]

I made a call graph to see how data flows through wrf into the urban scheme.

flowchart TD
    subgraph init
    start_domain --> |grid| start_domain_em
    start_domain_em --> |urb2d| phy_init
    phy_init --> |urb2d| bl_init
    bl_init --> |urb2d| urban_var_init
    bl_init --> |urb2d| urban_var_init
    end 

    subgraph physics
    solve_em --> |grid| first_rk_step_part1
    first_rk_step_part1 --> |urb2d| surface_driver

    surface_driver ==> |urb2d| clmdrv
    clmdrv --> |urb0d init | urban

    surface_driver ==> |urb2d| lsm
    lsm --> |urb0d | urban

    surface_driver ==> |urb2d| lsm_mosaic
    lsm_mosaic --> |urb0d | urban
    lsm_mosaic --> |urb0d | urban

    surface_driver --> |urb2d| noahmp_urban 
    noahmp_urban -->  |urb0d| urban
    end

Loading

Notice that:

• subroutine urban_var_init gets 2d arrays
• subroutine urban gets scalar values only, i.e. one grid point at a time
• all subroutines that call urban have a block that maps from 2d to scalar just before the call (and for some vars also after it)
• high in the tree, WRF simply passes on the grid. The grid contains all state variables, so we don't have to worry about passing on parameters down the call tree until the individual variables are first dereferenced from the grid. Therefore I excluded the upper part of the call tree.
• subroutine clmdrv sets most nudapt params to 0 (initial values??). The other callers of urban set it to actual values from the input array. They seem to be copies of each other.
• some routines have multiple calls to urban or urban_var_init. Again, they seem to be copies of each other.

[Peter9192]

More notes:

subroutine urban_var_init:
* initialize all 2d arrays with zeros(?) or some other relevant value

subroutine urban_param_init:
* Called by bl_init
* Parse urban parameters from table into local arrays
* table is parsed once and local arrays are stored in module scope of urban scheme
* if use_wudapt_lcz: use URBPARM.TBL
* if use_wudapt_lcz: use URBPARM_LCZ.TBL
* local arrays look like: HGT_TBL(ICATE), i.e. they all have _TBL and hold 3 or 10 categories
* some parameters are calculated based on values from the table, e.g. SVF_TBL (sky view factor for each urban category)

subroutine read_param:
* sets value of urban parameters from table
* get param from local table arrays by urban type
* e.g. SVF = SVF_TBL(UTYPE)

subroutine urban:
* call read_param: initially get all params from table
* nudapt modification: if mh_urb > 0,
* (re)calculate street and building width based on nudapt inputs
* overwrite urbparam table values with those from nudapt
* ZR = mh_urb (mean building height)
* R = lp_urb (plan area fraction)
* RW = 1.0-R (1 - plan area fraction)
* HGT = mh_urb/(BW+SW) (normalized height)
* SIGMA_ZED = stdh_urb (standard deviation of building height)
* lambda_f = lf_urb (frontal area index; something to do with wind dir)
* Lots of geometry calculations based on Macdonalds
ZDC, Z0C, lambda_fr, Z0R, Z0HC
* Calculate sky view factor VFWS, SVF