From 6a0787f80e8818a3d9fe5054efa5c1b00dd09af4 Mon Sep 17 00:00:00 2001 From: dmey Date: Sun, 5 Jan 2020 15:15:20 +0100 Subject: [PATCH] add wrf-teb integration --- Registry/Registry.EM_COMMON | 84 +- Registry/registry.dimspec | 6 + dyn_em/module_first_rk_step_part1.F | 44 + dyn_em/module_initialize_real.F | 2 +- dyn_em/start_em.F | 18 + phys/module_physics_init.F | 217 +++- phys/module_sf_noahdrv.F | 197 +++- phys/module_sf_psychrolib.F | 1430 +++++++++++++++++++++++++++ phys/module_sf_teb.F | 639 ++++++++++++ phys/module_sf_urban.F | 1074 +++++++++++++++++++- phys/module_surface_driver.F | 171 +++- run/URBPARM_TEB.TBL | 886 +++++++++++++++++ 12 files changed, 4736 insertions(+), 32 deletions(-) create mode 100644 phys/module_sf_psychrolib.F create mode 100644 phys/module_sf_teb.F create mode 100644 run/URBPARM_TEB.TBL diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 35200d3e9f..b24e632572 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -857,13 +857,13 @@ state real SMR_URB3D ilj misc 1 Z rd=(interp_m state real TRL_URB3D ilj misc 1 Z rd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "TRL_URB" "ROOF LAYER TEMPERATURE" "K" state real TBL_URB3D ilj misc 1 Z rd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "TBL_URB" "WALL LAYER TEMPERATURE" "K" state real TGL_URB3D ilj misc 1 Z rd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "TGL_URB" "ROAD LAYER TEMPERATURE" "K" -state real SH_URB2D ij misc 1 - rd=(interp_fcnm)u=(copy_fcnm) "SH_URB" "SENSIBLE HEAT FLUX FROM URBAN SFC" "W m{-2}" -state real LH_URB2D ij misc 1 - rd=(interp_fcnm)u=(copy_fcnm) "LH_URB" "LATENT HEAT FLUX FROM URBAN SFC" "W m{-2}" -state real G_URB2D ij misc 1 - rd=(interp_fcnm)u=(copy_fcnm) "G_URB" "GROUND HEAT FLUX INTO URBAN" "W m{-2}" -state real RN_URB2D ij misc 1 - rd=(interp_fcnm)u=(copy_fcnm) "RN_URB" "NET RADIATION ON URBAN SFC" "W m{-2}" -state real TS_URB2D ij misc 1 - rd=(interp_fcnm)u=(copy_fcnm) "TS_URB" "SKIN TEMPERATURE" "K" -state real FRC_URB2D ij misc 1 - i012rd=(interp_fcnm)u=(copy_fcnm) "FRC_URB2D" "URBAN FRACTION" "dimensionless" -state integer UTYPE_URB2D ij misc 1 - rd=(interp_fcnm)u=(copy_fcnm) "UTYPE_URB" "URBAN TYPE" "dimensionless" +state real SH_URB2D ij misc 1 - rhd=(interp_fcnm)u=(copy_fcnm) "SH_URB" "SENSIBLE HEAT FLUX FROM URBAN SFC" "W m{-2}" +state real LH_URB2D ij misc 1 - rhd=(interp_fcnm)u=(copy_fcnm) "LH_URB" "LATENT HEAT FLUX FROM URBAN SFC" "W m{-2}" +state real G_URB2D ij misc 1 - rhd=(interp_fcnm)u=(copy_fcnm) "G_URB" "GROUND HEAT FLUX INTO URBAN" "W m{-2}" +state real RN_URB2D ij misc 1 - rhd=(interp_fcnm)u=(copy_fcnm) "RN_URB" "NET RADIATION ON URBAN SFC" "W m{-2}" +state real TS_URB2D ij misc 1 - rhd=(interp_fcnm)u=(copy_fcnm) "TS_URB" "SKIN TEMPERATURE" "K" +state real FRC_URB2D ij misc 1 - i012rhd=(interp_fcnm)u=(copy_fcnm) "FRC_URB2D" "URBAN FRACTION" "dimensionless" +state integer UTYPE_URB2D ij misc 1 - rhd=(interp_fcnm)u=(copy_fcnm) "UTYPE_URB" "URBAN TYPE" "dimensionless" state real TRB_URB4D i{umap1}j misc 1 Z r "TRB_URB4D" "ROOF LAYER TEMPERATURE" "K" state real TW1_URB4D i{umap2}j misc 1 Z r "TW1_URB4D" "WALL LAYER TEMPERATURE" "K" state real TW2_URB4D i{umap2}j misc 1 Z r "TW2_URB4D" "WALL LAYER TEMPERATURE" "K" @@ -892,6 +892,56 @@ state real CHC_SFCDIF ij misc 1 - r "C state real CMGR_SFCDIF ij misc 1 - r "CMGR_SFCDIF" "" "" state real CHGR_SFCDIF ij misc 1 - r "CHGR_SFCDIF" "" "" +# Town Energy Balance (TEB) +state real TEB_TI_BLD ij misc 1 - r "TEB_TI_BLD" "inside building temp computed with its equation evolution" "K" +state real TEB_TI_BLD_EQ ij misc 1 - r "TEB_TI_BLD_EQ" "inside building temperature" "K" +state real TEB_TI_BLDWFR ij misc 1 - r "TEB_TI_BLDWFR" "inside building temperature without heating" "K" +state real TEB_T_CANYON ij misc 1 - r "TEB_T_CANYON" "outdoor temperature without heating" "K" +state real TEB_Q_CANYON ij misc 1 - r "TEB_Q_CANYON" "outdoor air specific humidity" "kg kg-1" +state real TEB_WS_ROOF ij misc 1 - rh "TEB_WS_ROOF" "roof water content" "(kg/m2)" +state real TEB_WS_ROAD ij misc 1 - rh "TEB_WS_ROAD" "road water content" "(kg/m2)" +state real TEB_WSNOW_ROOF i{teb_snly}j misc 1 Z r "TEB_WSNOW_ROOF" "roof snow layers reservoir" "kg m{-2}" +state real TEB_TSNOW_ROOF i{teb_snly}j misc 1 Z r "TEB_TSNOW_ROOF" "roof snow layers temperature" "K" +state real TEB_RSNOW_ROOF i{teb_snly}j misc 1 Z r "TEB_RSNOW_ROOF" "roof snow layers density" "kg m{-3}" +state real TEB_ASNOW_ROOF ij misc 1 - r "TEB_ASNOW_ROOF" "roof snow albedo" "dimensionless" +state real TEB_ESNOW_ROOF ij misc 1 - r "TEB_ESNOW_ROOF" "roof snow emissivity" "dimensionless" +state real TEB_TSSNOW_ROOF ij misc 1 - r "TEB_TSSNOW_ROOF" "roof snow surface temperature" "K" +state real TEB_WSNOW_ROAD i{teb_snly}j misc 1 Z r "TEB_WSNOW_ROAD" "road snow layers reservoir" "kg m{-2}" +state real TEB_TSNOW_ROAD i{teb_snly}j misc 1 Z r "TEB_TSNOW_ROAD" "road snow layers temperature" "K" +state real TEB_RSNOW_ROAD i{teb_snly}j misc 1 Z r "TEB_RSNOW_ROAD" "road snow layers density" "kg m{-3}" +state real TEB_ASNOW_ROAD ij misc 1 - r "TEB_ASNOW_ROAD" "road snow albedo" "dimensionless" +state real TEB_ESNOW_ROAD ij misc 1 - r "TEB_ESNOW_ROAD" "road snow emissivity" "dimensionless" +state real TEB_TSSNOW_ROAD ij misc 1 - r "TEB_TSSNOW_ROAD" "road snow surface temperature" "K" +state real TEB_T_WIN1 ij misc 1 - r "TEB_T_WIN1" "outdoor window temperature" "K" +state real TEB_T_WIN2 ij misc 1 - r "TEB_T_WIN2" "indoor window temperature" "K" +state real TEB_AUX_MAX ij misc 1 - r "TEB_AUX_MAX" "auxiliar variable for autosize calcs" "W m-2" +state real TEB_QI_BLD ij misc 1 - r "TEB_QI_BLD" "indoor air specific humidity" "kg kg-1" +state real TEB_THER_PRODC_DAY ij misc 1 - r "TEB_THER_PRODC_DAY" "present day integrated thermal production of energy" "J/m2" +state real TEB_T_FLOOR i{teb_flly}j misc 1 Z r "TEB_T_FLOOR" "floor layers temperatures" "K" +state real TEB_T_MASS i{teb_flly}j misc 1 Z r "TEB_T_MASS" "internal mass layers temperatures" "K" +state real TEB_T_ROAD i{teb_rdly}j misc 1 Z r "TEB_T_ROAD" "road layers temperatures" "K" +state real TEB_T_ROOF i{teb_roly}j misc 1 Z r "TEB_T_ROOF" "roof layers temperatures" "K" +state real TEB_T_WALL_A i{teb_waly}j misc 1 Z r "TEB_T_WALL_A" "wall layers temperatures (wall A)" "K" +state real TEB_T_WALL_B i{teb_waly}j misc 1 Z r "TEB_T_WALL_B" "wall layers temperatures (wall B)" "K" +state real TEB_HVAC_COOL ij misc 1 - h "TEB_HVAC_COOL" "energy consumption of the cooling system" "W m-2" +state real TEB_HVAC_HEAT ij misc 1 - h "TEB_HVAC_HEAT" "energy consumption of the heating system" "W m-2" +state real TEB_THER_PROD_PANEL ij misc 1 - h "TEB_THER_PROD_PANEL" "thermal energy production of solar panel on roofs" "W m-2" +state real TEB_PHOT_PROD_PANEL ij misc 1 - h "TEB_PHOT_PROD_PANEL" "photovoltaic energy production of solar panel on roofs" "W m-2" +state real TEB_TSK_RURAL ij misc 1 - r "TEB_TSK_RURAL" "TSK for rural fraction (TEB only)" "K" + +# BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE +state real TEB_INPUT_TA ij misc 1 - h "TEB_INPUT_TA" " " " " +state real TEB_INPUT_PS ij misc 1 - h "TEB_INPUT_PS" " " " " +state real TEB_INPUT_QA ij misc 1 - h "TEB_INPUT_QA" " " " " +state real TEB_INPUT_WIND ij misc 1 - h "TEB_INPUT_WIND" " " " " +state real TEB_INPUT_DIR ij misc 1 - h "TEB_INPUT_DIR" " " " " +state real TEB_INPUT_DIR_SW ij misc 1 - h "TEB_INPUT_DIR_SW" " " " " +state real TEB_INPUT_SCA_SW ij misc 1 - h "TEB_INPUT_SCA_SW" " " " " +state real TEB_INPUT_LW ij misc 1 - h "TEB_INPUT_LW" " " " " +state real TEB_INPUT_RAIN ij misc 1 - h "TEB_INPUT_RAIN" " " " " +state real TEB_INPUT_SNOW ij misc 1 - h "TEB_INPUT_SNOW" " " " " +state real TEB_INPUT_DIR_CO2 ij misc 1 - h "TEB_INPUT_DIR_CO2" " " " " +# END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE # solar location variables from radiation driver state real COSZEN ij misc 1 - rh "COSZEN" "COS of SOLAR ZENITH ANGLE" "dimensionless" @@ -956,7 +1006,7 @@ state real Z0 ij misc 1 - r "Z0 state real QZ0 ij misc 1 - r "QZ0" "SPECIFIC HUMIDITY AT ZNT" "kg kg-1" state real UZ0 ij misc 1 - r "UZ0" "U WIND COMPONENT AT ZNT" "m s-1" state real VZ0 ij misc 1 - r "VZ0" "V WIND COMPONENT AT ZNT" "m s-1" -state real QSFC ij misc 1 - r "QSFC" "SPECIFIC HUMIDITY AT LOWER BOUNDARY" "kg kg-1" +state real QSFC ij misc 1 - rh "QSFC" "SPECIFIC HUMIDITY AT LOWER BOUNDARY" "kg kg-1" state real AKHS ij misc 1 - r "AKHS" "SFC EXCH COEFF FOR HEAT" "m s-1" state real AKMS ij misc 1 - r "AKMS" "SFC EXCH COEFF FOR MOMENTUM" "m s-1" state integer KPBL ij misc 1 - r "KPBL" "LEVEL OF PBL TOP" "" @@ -1542,7 +1592,7 @@ state real SWDOWNC ij misc 1 - - "S state real GSW ij misc 1 - rd "GSW" "NET SHORT WAVE FLUX AT GROUND SURFACE" "W m-2" state real GLW ij misc 1 - rhd "GLW" "DOWNWARD LONG WAVE FLUX AT GROUND SURFACE" "W m-2" state real SWNORM ij misc 1 - rhd "SWNORM" "NORMAL SHORT WAVE FLUX AT GROUND SURFACE (SLOPE-DEPENDENT)" "W m-2" -state real diffuse_frac ij misc 1 - rd "DIFFUSE_FRAC" "DIFFUSE FRACTION OF SURFACE SHORTWAVE IRRADIANCE" "" +state real diffuse_frac ij misc 1 - rdh "DIFFUSE_FRAC" "DIFFUSE FRACTION OF SURFACE SHORTWAVE IRRADIANCE" "" # WRF-Solar state real swddir ij misc 1 - rd "SWDDIR" "Shortwave surface downward direct irradiance" "W m-2" "" state real swddirc ij misc 1 - rd "SWDDIRC" "Clear-sky Shortwave surface downward direct irradiance" "W m-2" "" @@ -1683,8 +1733,8 @@ state integer STEPRA - misc 1 - r "S state real RUBLTEN ikj misc 1 - r "RUBLTEN" "X WIND TENDENCY DUE TO PBL PARAMETERIZATION" "m s-2" state real RVBLTEN ikj misc 1 - r "RVBLTEN" "Y WIND TENDENCY DUE TO PBL PARAMETERIZATION" "m s-2" -state real RTHBLTEN ikj misc 1 - r "RTHBLTEN" "THETA TENDENCY DUE TO PBL PARAMETERIZATION" "K s-1" -state real RQVBLTEN ikj misc 1 - r "RQVBLTEN" "Q_V TENDENCY DUE TO PBL PARAMETERIZATION" "kg kg-1 s-1" +state real RTHBLTEN ikj misc 1 - rh "RTHBLTEN" "THETA TENDENCY DUE TO PBL PARAMETERIZATION" "K s-1" +state real RQVBLTEN ikj misc 1 - rh "RQVBLTEN" "Q_V TENDENCY DUE TO PBL PARAMETERIZATION" "kg kg-1 s-1" state real RQCBLTEN ikj misc 1 - r "RQCBLTEN" "Q_C TENDENCY DUE TO PBL PARAMETERIZATION" "kg kg-1 s-1" state real RQIBLTEN ikj misc 1 - r "RQIBLTEN" "Q_I TENDENCY DUE TO PBL PARAMETERIZATION" "kg kg-1 s-1" state real RQNIBLTEN ikj misc 1 - r "RQNIBLTEN" "Q_NI TENDENCY DUE TO PBL PARAMETERIZATION" "# kg-1 s-1" @@ -2279,6 +2329,14 @@ rconfig integer maxpatch namelist,physics 1 10 rconfig integer num_snow_layers namelist,physics 1 3 irh "num_snow_layers" "" "" rconfig integer num_snso_layers namelist,physics 1 7 irh "num_snso_layers" "" "" +# TEB +rconfig integer teb_num_floor_layers namelist,physics 1 5 irh "teb_num_floor_layers" "" "" +rconfig integer teb_num_road_layers namelist,physics 1 5 irh "teb_num_road_layers" "" "" +rconfig integer teb_num_roof_layers namelist,physics 1 5 irh "teb_num_roof_layers" "" "" +rconfig integer teb_num_wall_layers namelist,physics 1 5 irh "teb_num_wall_layers" "" "" +rconfig integer teb_num_snow_layers namelist,physics 1 1 irh "teb_num_snow_layers" "" "" +rconfig integer teb_test_integration namelist,physics 1 0 irh "teb_test_integration" "" "" + rconfig integer num_urban_ndm derived 1 1 irh "num_urban_ndm" "maximum number of street dimensions (ndm in BEP or BEM header)" "" rconfig integer num_urban_ng derived 1 1 irh "num_urban_ng" "number of grid levels in the ground (ng_u in BEP or BEM header)" "" rconfig integer num_urban_nwr derived 1 1 irh "num_urban_nwr" "number of grid levels in the walls or roof (nwr_u in BEP or BEM header)" "" @@ -2856,6 +2914,10 @@ package sfclayscheme sf_sfclay_physics==91 - - package noahucmscheme sf_urban_physics==1 - state:trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,a_u_bep,a_v_bep,a_t_bep,a_q_bep,a_e_bep,b_u_bep,b_v_bep,b_t_bep,b_q_bep,b_e_bep,dlg_bep,dl_u_bep,sf_bep,vl_bep,mh_urb2d,stdh_urb2d,lf_urb2d,lp_urb2d,hgt_urb2d,lb_urb2d,tgr_urb2d,cmcr_urb2d,drelr_urb2d,drelb_urb2d,drelg_urb2d,flxhumr_urb2d,flxhumb_urb2d,flxhumg_urb2d,tgrl_urb3d,smr_urb3d,cmgr_sfcdif,chgr_sfcdif,trl_urb3d,tgl_urb3d,tbl_urb3d package bepscheme sf_urban_physics==2 - state:a_u_bep,a_v_bep,a_t_bep,a_q_bep,a_e_bep,b_u_bep,b_v_bep,b_t_bep,b_q_bep,b_e_bep,dlg_bep,dl_u_bep,sf_bep,vl_bep,trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,hi_urb2d,lp_urb2d,hgt_urb2d,lb_urb2d,trl_urb3d,tgl_urb3d,tbl_urb3d,tsk_rural package bep_bemscheme sf_urban_physics==3 - state:a_u_bep,a_v_bep,a_t_bep,a_q_bep,a_e_bep,b_u_bep,b_v_bep,b_t_bep,b_q_bep,b_e_bep,dlg_bep,dl_u_bep,sf_bep,vl_bep,trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d,tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d,tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d,cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d,sfwin1_urb3d,sfwin2_urb3d,sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d,hi_urb2d,lp_urb2d,hgt_urb2d,lb_urb2d,trl_urb3d,tgl_urb3d,tbl_urb3d,tsk_rural +package tebscheme sf_urban_physics==4 - state:lp_urb2d,hgt_urb2d,lb_urb2d,trl_urb3d,tgl_urb3d,tbl_urb3d,TEB_TI_BLD,TEB_T_CANYON,TEB_Q_CANYON,TEB_TI_BLD_EQ,TEB_TI_BLDWFR,TEB_WS_ROOF,TEB_WS_ROAD,TEB_WSNOW_ROOF,TEB_TSNOW_ROOF,TEB_RSNOW_ROOF,TEB_ASNOW_ROOF,TEB_TSSNOW_ROOF,TEB_ESNOW_ROOF,TEB_WSNOW_ROAD,TEB_TSNOW_ROAD,TEB_RSNOW_ROAD,TEB_ASNOW_ROAD,TEB_TSSNOW_ROAD,TEB_ESNOW_ROAD,TEB_T_WIN1,TEB_T_WIN2,TEB_AUX_MAX,TEB_QI_BLD,TEB_THER_PRODC_DAY,TEB_T_FLOOR,TEB_T_MASS,TEB_T_ROAD,TEB_T_ROOF,TEB_T_WALL_A,TEB_T_WALL_B,TEB_HVAC_COOL,TEB_HVAC_HEAT,TEB_TSK_RURAL + +package teb_notest teb_test_integration==0 - - +package teb_test teb_test_integration==1 - state:TEB_INPUT_TA,TEB_INPUT_PS,TEB_INPUT_QA,TEB_INPUT_WIND,TEB_INPUT_DIR,TEB_INPUT_DIR_SW,TEB_INPUT_SCA_SW,TEB_INPUT_LW,TEB_INPUT_RAIN,TEB_INPUT_SNOW,TEB_INPUT_DIR_CO2 package nolsmscheme sf_surface_physics==0 - - package slabscheme sf_surface_physics==1 - - diff --git a/Registry/registry.dimspec b/Registry/registry.dimspec index a680a6fce7..456b50a906 100644 --- a/Registry/registry.dimspec +++ b/Registry/registry.dimspec @@ -53,6 +53,12 @@ dimspec snly 2 namelist=num_snow_layers z snow_layers dimspec l 2 namelist=num_soil_layers z soil_layers dimspec snsl 2 namelist=num_snso_layers z snso_layers +dimspec teb_flly 2 namelist=teb_num_floor_layers z teb_floor_layers +dimspec teb_rdly 2 namelist=teb_num_road_layers z teb_road_layers +dimspec teb_roly 2 namelist=teb_num_roof_layers z teb_roof_layers +dimspec teb_waly 2 namelist=teb_num_wall_layers z teb_wall_layers +dimspec teb_snly 2 namelist=teb_num_snow_layers z teb_snow_layers + dimspec umap0 2 namelist=num_urban_ndm z num_urban_ndm dimspec umap1 2 namelist=urban_map_zrd z urban_map_zrd dimspec umap2 2 namelist=urban_map_zwd z urban_map_zwd diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index f2737e86b9..9952434805 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -690,6 +690,50 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,SF_BEP=grid%sf_bep,VL_BEP=grid%vl_bep & & ,A_E_BEP=grid%a_e_bep,B_E_BEP=grid%b_e_bep,DLG_BEP=grid%dlg_bep & & ,DL_U_BEP=grid%dl_u_bep & + ! + ! BEGIN: Optional urban for TEB + & ,teb_test_integration=grid%teb_test_integration & + & ,teb_num_snow_layers=grid%teb_num_snow_layers & + & ,teb_num_road_layers=grid%teb_num_road_layers & + & ,teb_num_roof_layers=grid%teb_num_roof_layers & + & ,teb_num_wall_layers=grid%teb_num_wall_layers & + & ,teb_num_floor_layers=grid%teb_num_floor_layers & + & ,TEB_TI_BLD=grid%TEB_TI_BLD & + & ,TEB_T_CANYON=grid%TEB_T_CANYON, TEB_Q_CANYON=grid%TEB_Q_CANYON & + & ,TEB_WS_ROOF=grid%TEB_WS_ROOF, TEB_WS_ROAD=grid%TEB_WS_ROAD & + & ,TEB_WSNOW_ROOF=grid%TEB_WSNOW_ROOF, TEB_TSNOW_ROOF=grid%TEB_TSNOW_ROOF & + & ,TEB_RSNOW_ROOF=grid%TEB_RSNOW_ROOF, TEB_ASNOW_ROOF=grid%TEB_ASNOW_ROOF & + & ,TEB_TSSNOW_ROOF=grid%TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF=grid%TEB_ESNOW_ROOF & + & ,TEB_WSNOW_ROAD=grid%TEB_WSNOW_ROAD, TEB_TSNOW_ROAD=grid%TEB_TSNOW_ROAD & + & ,TEB_RSNOW_ROAD=grid%TEB_RSNOW_ROAD, TEB_ASNOW_ROAD=grid%TEB_ASNOW_ROAD & + & ,TEB_TSSNOW_ROAD=grid%TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD=grid%TEB_ESNOW_ROAD & + & ,TEB_T_WIN1=grid%TEB_T_WIN1, TEB_T_WIN2=grid%TEB_T_WIN2 & + & ,TEB_AUX_MAX=grid%TEB_AUX_MAX, TEB_THER_PRODC_DAY=grid%TEB_THER_PRODC_DAY & + & ,TEB_QI_BLD=grid%TEB_QI_BLD & + & ,TEB_T_FLOOR=grid%TEB_T_FLOOR, TEB_T_MASS=grid%TEB_T_MASS & + & ,TEB_T_ROAD=grid%TEB_T_ROAD, TEB_T_ROOF=grid%TEB_T_ROOF & + & ,TEB_T_WALL_A=grid%TEB_T_WALL_A, TEB_T_WALL_B=grid%TEB_T_WALL_B & + & ,TEB_HVAC_COOL=grid%TEB_HVAC_COOL, TEB_HVAC_HEAT=grid%TEB_HVAC_HEAT & + & ,TEB_THER_PROD_PANEL=grid%TEB_THER_PROD_PANEL & + & ,TEB_PHOT_PROD_PANEL=grid%TEB_PHOT_PROD_PANEL & + & ,TEB_TSK_RURAL=grid%TEB_TSK_RURAL & + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + & ,TEB_INPUT_TA=grid%TEB_INPUT_TA & + & ,TEB_INPUT_PS=grid%TEB_INPUT_PS & + & ,TEB_INPUT_QA=grid%TEB_INPUT_QA & + & ,TEB_INPUT_WIND=grid%TEB_INPUT_WIND & + & ,TEB_INPUT_DIR=grid%TEB_INPUT_DIR & + & ,TEB_INPUT_DIR_SW=grid%TEB_INPUT_DIR_SW & + & ,TEB_INPUT_SCA_SW=grid%TEB_INPUT_SCA_SW & + & ,TEB_INPUT_LW=grid%TEB_INPUT_LW & + & ,TEB_INPUT_RAIN=grid%TEB_INPUT_RAIN & + & ,TEB_INPUT_SNOW=grid%TEB_INPUT_SNOW & + & ,TEB_INPUT_DIR_CO2=grid%TEB_INPUT_DIR_CO2 & + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + ! + ! END: Optional urban for TEB + ! & ,CMR_SFCDIF=grid%cmr_sfcdif, CHR_SFCDIF=grid%chr_sfcdif & !I/O urban & ,CMC_SFCDIF=grid%cmc_sfcdif, CHC_SFCDIF=grid%chc_sfcdif & !I/O urban & ,CMGR_SFCDIF=grid%cmgr_sfcdif, CHGR_SFCDIF=grid%chgr_sfcdif & !I/O urban diff --git a/dyn_em/module_initialize_real.F b/dyn_em/module_initialize_real.F index 279a789d51..cc16206381 100644 --- a/dyn_em/module_initialize_real.F +++ b/dyn_em/module_initialize_real.F @@ -3060,7 +3060,7 @@ SUBROUTINE init_domain_rk ( grid & ! Split NUDAPT Urban Parameters - IF ( ( config_flags%sf_urban_physics == 1 ) .OR. ( config_flags%sf_urban_physics == 2 ) .OR. ( config_flags%sf_urban_physics == 3 ) ) THEN + IF ( ( config_flags%sf_urban_physics == 1 ) .OR. ( config_flags%sf_urban_physics == 2 ) .OR. ( config_flags%sf_urban_physics == 3 ) .OR. ( config_flags%sf_urban_physics == 4 ) ) THEN DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) IF ( MMINLU == 'NLCD40' .OR. MMINLU == 'MODIFIED_IGBP_MODIS_NOAH') THEN diff --git a/dyn_em/start_em.F b/dyn_em/start_em.F index ae8447951c..0de31e129a 100644 --- a/dyn_em/start_em.F +++ b/dyn_em/start_em.F @@ -1097,6 +1097,24 @@ SUBROUTINE start_domain_em ( grid, allowed_to_read & grid%A_E_BEP,grid%B_U_BEP,grid%B_V_BEP,grid%B_T_BEP, & !multi-layer urban grid%B_Q_BEP,grid%B_E_BEP,grid%DLG_BEP, & !multi-layer urban grid%DL_U_BEP,grid%SF_BEP,grid%VL_BEP, & !multi-layer urban + ! optional urban for TEB + grid%teb_num_snow_layers, & + grid%teb_num_road_layers, & + grid%teb_num_roof_layers, & + grid%teb_num_wall_layers, & + grid%teb_num_floor_layers, & + grid%TEB_TI_BLD, & + grid%TEB_T_CANYON, grid%TEB_Q_CANYON, & + grid%TEB_WS_ROOF, grid%TEB_WS_ROAD, & + grid%TEB_WSNOW_ROOF, grid%TEB_TSNOW_ROOF, grid%TEB_RSNOW_ROOF, grid%TEB_ASNOW_ROOF, & + grid%TEB_TSSNOW_ROOF, grid%TEB_ESNOW_ROOF, & + grid%TEB_WSNOW_ROAD, grid%TEB_TSNOW_ROAD, grid%TEB_RSNOW_ROAD, grid%TEB_ASNOW_ROAD, & + grid%TEB_TSSNOW_ROAD, grid%TEB_ESNOW_ROAD, & + grid%TEB_T_WIN1, grid%TEB_T_WIN2, grid%TEB_AUX_MAX, grid%TEB_THER_PRODC_DAY, & + grid%TEB_QI_BLD, grid%TEB_T_FLOOR, grid%TEB_T_MASS, & + grid%TEB_T_ROAD, grid%TEB_T_ROOF, grid%TEB_T_WALL_A, grid%TEB_T_WALL_B, & + grid%TEB_TSK_RURAL, & + ! end for TEB grid%TML,grid%T0ML,grid%HML,grid%H0ML,grid%HUML,grid%HVML,grid%TMOML, & !Optional oml grid%lakedepth2d, grid%savedtke12d, grid%snowdp2d, grid%h2osno2d, & !lake grid%snl2d, grid%t_grnd2d, grid%t_lake3d, grid%lake_icefrac3d, & !lake diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index db727f1201..43e336eef8 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -155,6 +155,29 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & A_E_BEP,B_U_BEP,B_V_BEP, & !Optional multi-layer urban B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & !Optional multi-layer urban DL_U_BEP,SF_BEP,VL_BEP, & !Optional multi-layer urban + ! optional urban for TEB + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers, & + TEB_TI_BLD, & + TEB_T_CANYON,TEB_Q_CANYON, & + TEB_WS_ROOF, TEB_WS_ROAD, & + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & + TEB_T_WIN1, TEB_T_WIN2, & + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & + TEB_QI_BLD, & + TEB_T_FLOOR, TEB_T_MASS, & + TEB_T_ROAD, TEB_T_ROOF, & + TEB_T_WALL_A, TEB_T_WALL_B, & + TEB_TSK_RURAL, & + ! end for TEB TML,T0ML,HML,H0ML,HUML,HVML,TMOML, & !Optional oml lakedepth2d, savedtke12d, snowdp2d, h2osno2d, & !lake snl2d, t_grnd2d, t_lake3d, lake_icefrac3d, & !lake @@ -667,6 +690,42 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & REAL, OPTIONAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DLG_BEP REAL, OPTIONAL, DIMENSION(ims:ime, kms:kme,jms:jme), INTENT(INOUT) :: SF_BEP REAL, OPTIONAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DL_U_BEP + !for TEB + INTEGER, INTENT(IN) :: teb_num_snow_layers + INTEGER, INTENT(IN) :: teb_num_road_layers + INTEGER, INTENT(IN) :: teb_num_roof_layers + INTEGER, INTENT(IN) :: teb_num_wall_layers + INTEGER, INTENT(IN) :: teb_num_floor_layers + REAL, OPTIONAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TEB_TI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TEB_T_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TEB_Q_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TEB_WS_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TEB_WS_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN1 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN2 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_AUX_MAX + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_THER_PRODC_DAY + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_QI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_FLOOR + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_MASS + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_A + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_B + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSK_RURAL + !end TEB ! lake varibles: real, dimension(ims:ime,jms:jme ),intent(out) :: lakedepth2d, & savedtke12d @@ -1350,6 +1409,29 @@ SUBROUTINE phy_init ( id, config_flags, DT, restart, zfull, zhalf, & A_E_BEP,B_U_BEP,B_V_BEP, & !Optional multi-layer urban B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & !Optional multi-layer urban DL_U_BEP,SF_BEP,VL_BEP, & !Optional multi-layer urban + ! optional urban for TEB + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers, & + TEB_TI_BLD, & + TEB_T_CANYON,TEB_Q_CANYON, & + TEB_WS_ROOF, TEB_WS_ROAD, & + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & + TEB_T_WIN1, TEB_T_WIN2, & + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & + TEB_QI_BLD, & + TEB_T_FLOOR, TEB_T_MASS, & + TEB_T_ROAD, TEB_T_ROOF, & + TEB_T_WALL_A, TEB_T_WALL_B, & + TEB_TSK_RURAL, & + ! end for TEB ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & @@ -2341,6 +2423,29 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & A_E_BEP,B_U_BEP,B_V_BEP, & !Optional multi-layer urban B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & !Optional multi-layer urban DL_U_BEP,SF_BEP,VL_BEP, & !Optional multi-layer urban + ! optional urban for TEB + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers, & + TEB_TI_BLD, & + TEB_T_CANYON,TEB_Q_CANYON, & + TEB_WS_ROOF, TEB_WS_ROAD, & + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & + TEB_T_WIN1, TEB_T_WIN2, & + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & + TEB_QI_BLD, & + TEB_T_FLOOR, TEB_T_MASS, & + TEB_T_ROAD, TEB_T_ROOF, & + TEB_T_WALL_A, TEB_T_WALL_B, & + TEB_TSK_RURAL, & + ! end for TEB ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte, & @@ -2697,6 +2802,42 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & REAL, OPTIONAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DLG_BEP REAL, OPTIONAL, DIMENSION(ims:ime, kms:kme,jms:jme),INTENT(INOUT) :: SF_BEP REAL, OPTIONAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DL_U_BEP + !for TEB + INTEGER, INTENT(IN) :: teb_num_snow_layers + INTEGER, INTENT(IN) :: teb_num_road_layers + INTEGER, INTENT(IN) :: teb_num_roof_layers + INTEGER, INTENT(IN) :: teb_num_wall_layers + INTEGER, INTENT(IN) :: teb_num_floor_layers + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_Q_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN1 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN2 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_AUX_MAX + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_THER_PRODC_DAY + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_QI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_FLOOR + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_MASS + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_A + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_B + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSK_RURAL + !end TEB ! lake varibles: real, dimension(ims:ime,jms:jme ),intent(out) :: lakedepth2d, & savedtke12d @@ -2993,14 +3134,19 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & its,ite, jts,jte, kts,kte ) -!URBAN - IF ((SF_URBAN_PHYSICS.eq.1).OR.(SF_URBAN_PHYSICS.EQ.2).OR.(SF_URBAN_PHYSICS.EQ.3)) THEN +!URBAN. Add TEB with option 4. + IF ((SF_URBAN_PHYSICS.eq.1).OR.(SF_URBAN_PHYSICS.EQ.2).OR.(SF_URBAN_PHYSICS.EQ.3).OR.(SF_URBAN_PHYSICS.EQ.4)) THEN IF ( PRESENT( FRC_URB2D ) .AND. PRESENT( UTYPE_URB2D )) THEN CALL urban_param_init(DZR,DZB,DZG,num_soil_layers, & !urban - sf_urban_physics) -! num_roof_layers,num_wall_layers,road_soil_layers) !urban + sf_urban_physics, & +! num_roof_layers,num_wall_layers,road_soil_layers, !urban + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers) CALL urban_var_init(ISURBAN,TSK,TSLB,TMN,IVGTYP, & !urban @@ -3043,7 +3189,31 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & A_E_BEP,B_U_BEP,B_V_BEP, & !multi-layer urban B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & !multi-layer urban DL_U_BEP,SF_BEP,VL_BEP, & !multi-layer urban - FRC_URB2D, UTYPE_URB2D) !urban + FRC_URB2D, UTYPE_URB2D, & !urban + ! for TEB + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers, & + TEB_TI_BLD, & + TEB_T_CANYON, TEB_Q_CANYON, & + TEB_WS_ROOF, TEB_WS_ROAD, & + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & + TEB_T_WIN1, TEB_T_WIN2, & + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & + TEB_QI_BLD, & + TEB_T_FLOOR, TEB_T_MASS, & + TEB_T_ROAD, TEB_T_ROOF, & + TEB_T_WALL_A, TEB_T_WALL_B, & + TEB_TSK_RURAL & + ! end TEB + ) ELSE CALL wrf_error_fatal ( 'arguments not present for calling urban model' ) ENDIF @@ -3119,7 +3289,12 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & IF ((SF_URBAN_PHYSICS.eq.1).OR.(SF_URBAN_PHYSICS.EQ.2).OR.(SF_URBAN_PHYSICS.EQ.3)) THEN IF ( PRESENT( FRC_URB2D ) .AND. PRESENT( UTYPE_URB2D )) THEN CALL urban_param_init(DZR,DZB,DZG,num_soil_layers, & !urban - sf_urban_physics) + sf_urban_physics, & + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers) CALL urban_var_init(ISURBAN,TSK,TSLB,TMN,IVGTYP, & !urban ims,ime,jms,jme,kms,kme,num_soil_layers, & !urban LOW_DENSITY_RESIDENTIAL_TABLE, & @@ -3159,10 +3334,36 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & A_E_BEP,B_U_BEP,B_V_BEP, & !multi-layer urban B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & !multi-layer urban DL_U_BEP,SF_BEP,VL_BEP, & !multi-layer urban - FRC_URB2D, UTYPE_URB2D) !urban + FRC_URB2D, UTYPE_URB2D, & !urban + !for TEB + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers, & + TEB_TI_BLD, & + TEB_T_CANYON, TEB_Q_CANYON, & + TEB_WS_ROOF, TEB_WS_ROAD, & + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & + TEB_T_WIN1, TEB_T_WIN2, & + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & + TEB_QI_BLD, & + TEB_T_FLOOR, TEB_T_MASS, & + TEB_T_ROAD, TEB_T_ROOF, & + TEB_T_WALL_A, TEB_T_WALL_B, & + TEB_TSK_RURAL & + ! end TEB + ) ELSE CALL wrf_error_fatal ( 'arguments not present for calling urban model' ) ENDIF + ELSE + CALL wrf_error_fatal('TEB does not support NOAHMP') ENDIF CASE (RUCLSMSCHEME) @@ -3215,7 +3416,7 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & !-------------------------------------------------------------- ! CLM Init Coupling CASE (CLMSCHEME) - IF ((SF_URBAN_PHYSICS.eq.1).OR.(SF_URBAN_PHYSICS.EQ.2).OR.(SF_URBAN_PHYSICS.EQ.3)) THEN + IF ((SF_URBAN_PHYSICS.eq.1).OR.(SF_URBAN_PHYSICS.EQ.2).OR.(SF_URBAN_PHYSICS.EQ.3).OR.(SF_URBAN_PHYSICS.EQ.4)) THEN CALL wrf_error_fatal ( 'CLM DOES NOT WORK WITH URBAN SCHEME' ) ENDIF IF ( TRIM(mminlu) .EQ. 'NLCD40' ) THEN diff --git a/phys/module_sf_noahdrv.F b/phys/module_sf_noahdrv.F index 50c5809a99..2d3d58c0a0 100644 --- a/phys/module_sf_noahdrv.F +++ b/phys/module_sf_noahdrv.F @@ -16,6 +16,7 @@ MODULE module_sf_noahdrv USE module_sf_noahlsm_glacial_only, only: sflx_glacial USE module_sf_bep, only: bep USE module_sf_bep_bem, only: bep_bem + USE module_sf_teb, only: teb USE module_ra_gfdleta, only: cal_mon_day #if ( WRF_CHEM == 1 ) USE module_data_gocart_dust @@ -102,7 +103,53 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban a_e_bep,b_u_bep,b_v_bep, & !O multi-layer urban b_t_bep,b_q_bep,b_e_bep,dlg_bep, & !O multi-layer urban - dl_u_bep,sf_bep,vl_bep,sfcheadrt,INFXSRT, soldrain & !O multi-layer urban + dl_u_bep,sf_bep,vl_bep, & + ! + ! BEGIN: Optional urban for TEB + xtime, & !I + teb_test_integration, & !I + teb_num_snow_layers, & !I + teb_num_road_layers, & !I + teb_num_roof_layers, & !I + teb_num_wall_layers, & !I + teb_num_floor_layers, & !I + diffuse_frac, & !I + TEB_TI_BLD, & !H + TEB_T_CANYON, TEB_Q_CANYON, & !H + TEB_WS_ROOF, TEB_WS_ROAD, & !H + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & !H + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & !H + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & !H + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & !H + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & !H + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & !H + TEB_T_WIN1, TEB_T_WIN2, & !H + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & !H + TEB_QI_BLD, & !H + TEB_T_FLOOR, TEB_T_MASS, & !H + TEB_T_ROAD, TEB_T_ROOF, & !H + TEB_T_WALL_A, TEB_T_WALL_B, & !H + TEB_TSK_RURAL, & !H + TEB_HVAC_COOL, TEB_HVAC_HEAT, & !O + TEB_THER_PROD_PANEL, TEB_PHOT_PROD_PANEL, & !O + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + TEB_INPUT_TA, & + TEB_INPUT_PS, & + TEB_INPUT_QA, & + TEB_INPUT_WIND, & + TEB_INPUT_DIR, & + TEB_INPUT_DIR_SW, & + TEB_INPUT_SCA_SW, & + TEB_INPUT_LW, & + TEB_INPUT_RAIN, & + TEB_INPUT_SNOW, & + TEB_INPUT_DIR_CO2, & + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + ! + ! END: Optional urban for TEB + ! + sfcheadrt,INFXSRT, soldrain & !O multi-layer urban ,SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM, fasdas & !fasdas ,RC2,XLAI2 & ) @@ -632,9 +679,75 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep !Fraction air at the face of grid cell REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep !Length scale +!for TEB + INTEGER, INTENT(IN) :: teb_test_integration + INTEGER, INTENT(IN) :: teb_num_snow_layers + INTEGER, INTENT(IN) :: teb_num_road_layers + INTEGER, INTENT(IN) :: teb_num_roof_layers + INTEGER, INTENT(IN) :: teb_num_wall_layers + INTEGER, INTENT(IN) :: teb_num_floor_layers + REAL, OPTIONAL, INTENT(IN) :: XTIME + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: diffuse_frac + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_Q_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN1 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN2 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_AUX_MAX + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_THER_PRODC_DAY + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_QI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_FLOOR + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_MASS + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_A + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_B + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSK_RURAL + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_HVAC_COOL + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_HVAC_HEAT + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_THER_PROD_PANEL + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_PHOT_PROD_PANEL + ! + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_TA + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_PS + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_QA + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_WIND + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_DIR + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_DIR_SW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_SCA_SW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_LW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_RAIN + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_SNOW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_DIR_CO2 + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + ! + ! local variables for time calculations + INTEGER(KIND=8) :: ixhour + REAL(KIND=8) :: xhour + REAL :: xmin, gmtp + +!end TEB + + ! Local variables for multi-layer UCM (Martilli et al. 2002) REAL, DIMENSION( its:ite, jts:jte ) :: HFX_RURAL,LH_RURAL,GRDFLX_RURAL ! ,RN_RURAL - REAL, DIMENSION( its:ite, jts:jte ) :: QFX_RURAL ! ,QSFC_RURAL,UMOM_RURAL,VMOM_RURAL + REAL, DIMENSION( its:ite, jts:jte ) :: QFX_RURAL, QSFC_RURAL !,UMOM_RURAL,VMOM_RURAL REAL, DIMENSION( its:ite, jts:jte ) :: ALB_RURAL,EMISS_RURAL,TSK_RURAL ! ,UST_RURAL ! REAL, DIMENSION( ims:ime, jms:jme ) :: QSFC_URB REAL, DIMENSION( its:ite, jts:jte ) :: HFX_URB,UMOM_URB,VMOM_URB @@ -819,6 +932,8 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & LH_RURAL(I,J)=LH(I,J) EMISS_RURAL(I,J)=EMISS(I,J) GRDFLX_RURAL(I,J)=GRDFLX(I,J) + QSFC_RURAL(I,J)=QSFC(I,J) + ALB_RURAL(I,J)=ALBEDO(I,J) ELSE ! Land or sea-ice case @@ -907,7 +1022,7 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as ! the "NATURAL" category in the VEGPARM.TBL - IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN + IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 .OR. SF_URBAN_PHYSICS==4) THEN ! TEB opt 4 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LOW_DENSITY_RESIDENTIAL .or. & IVGTYP(I,J) == HIGH_DENSITY_RESIDENTIAL .or. IVGTYP(I,J) == HIGH_INTENSITY_INDUSTRIAL) THEN VEGTYP = NATURAL @@ -920,6 +1035,8 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J)) elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then T1=tsk_rural_bep(i,j) + elseif(sf_urban_physics.eq.4) then + T1 = TEB_TSK_RURAL(i,j) endif ELSE T1 = TSK(I,J) @@ -959,14 +1076,14 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & ENDIF ENDIF - IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN + IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3 .or. SF_URBAN_PHYSICS == 4) THEN IF(AOASIS > 1.0) THEN CALL wrf_error_fatal('Urban oasis option is for SF_URBAN_PHYSICS == 1 only') ENDIF IF(IRIOPTION == 1) THEN CALL wrf_error_fatal('Urban irrigation option is for SF_URBAN_PHYSICS == 1 only') ENDIF - ENDIF + ENDIF #if 0 IF(IPRINT) THEN @@ -1169,6 +1286,8 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & TSK_RURAL(I,J)=T1 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN TSK_RURAL_BEP(I,J)=T1 + ELSEIF (SF_URBAN_PHYSICS == 4) THEN + TEB_TSK_RURAL(I,J) = T1 END IF HFX(I,J)=SHEAT HFX_RURAL(I,J)=SHEAT @@ -1202,8 +1321,7 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & ! QSFC(I,J)=Q1 ! Convert QSFC back to mixing ratio QSFC(I,J)= Q1/(1.0-Q1) -! - ! QSFC_RURAL(I,J)= Q1/(1.0-Q1) + QSFC_RURAL(I,J)= Q1/(1.0-Q1) ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002) DO 80 NS=1,NSOIL @@ -1691,6 +1809,69 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & endif !Bep end +! for TEB +! + if (sf_urban_physics.eq.4) then + + CALL TEB( & + ims ,ime, jms, jme, kms, kme, & + its ,ite, jts, jte, kts, kte, & + teb_test_integration, & + teb_num_roof_layers, teb_num_wall_layers, teb_num_road_layers, & + teb_num_snow_layers, teb_num_floor_layers, & + XLAT, XLONG, & + julian, julyr, gmt, xtime, DT, & + U_PHY, V_PHY, & + DZ8W, P8W3D, T3D, QV3D, RHO, & + SWDOWN, diffuse_frac, GLW, & + RAINBL, SR, FRC_URB2D, UTYPE_URB2D, & + IVGTYP, ISURBAN, & + CQS2, CHS2, & + ! input from noah + ALB_RURAL, EMISS_RURAL, & + HFX_RURAL, QFX_RURAL, LH_RURAL, GRDFLX_RURAL, & + TEB_TSK_RURAL, QSFC_RURAL, & + ! input/output from TEB + TEB_TI_BLD, & + TEB_T_CANYON, TEB_Q_CANYON, & + TEB_WS_ROOF, TEB_WS_ROAD, & + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & + TEB_T_WIN1, TEB_T_WIN2, & + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & + TEB_QI_BLD, & + TEB_T_FLOOR, TEB_T_MASS, & + TEB_T_ROAD, TEB_T_ROOF, & + TEB_T_WALL_A, TEB_T_WALL_B, & + ! output of TEB + TEB_HVAC_COOL, TEB_HVAC_HEAT, & + TEB_THER_PROD_PANEL, TEB_PHOT_PROD_PANEL, & + RN_URB2D, U10_URB2D, V10_URB2D, & + ! output Noah + TEB + ALBEDO, HFX, QFX, LH, GRDFLX, TSK, EMISS, QSFC, UST, & + TH2_URB2D, Q2_URB2D, & + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + TEB_INPUT_TA, & + TEB_INPUT_PS, & + TEB_INPUT_QA, & + TEB_INPUT_WIND, & + TEB_INPUT_DIR, & + TEB_INPUT_DIR_SW, & + TEB_INPUT_SCA_SW, & + TEB_INPUT_LW, & + TEB_INPUT_RAIN, & + TEB_INPUT_SNOW, & + TEB_INPUT_DIR_CO2) + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + + endif +! end TEB + !------------------------------------------------------ END SUBROUTINE lsm !------------------------------------------------------ @@ -4158,6 +4339,8 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J)) elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then T1=tsk_rural_bep(i,j) + elseif(sf_urban_physics.eq.4)then + CALL wrf_error_fatal('NoahMP is not supported for TEB') endif ELSE T1 = TSK(I,J) diff --git a/phys/module_sf_psychrolib.F b/phys/module_sf_psychrolib.F new file mode 100644 index 0000000000..23906122f5 --- /dev/null +++ b/phys/module_sf_psychrolib.F @@ -0,0 +1,1430 @@ +! PsychroLib version 2.1.1 (https://github.com/psychrometrics/psychrolib) +! Copyright (c) 2018 D. Thevenard and D. Meyer for the current library implementation +! Copyright (c) 2017 ASHRAE Handbook — Fundamentals for ASHRAE equations and coefficients +! Licensed under the MIT License. + +! Please cite paper as https://doi.org/10.21105/joss.01137. +! For version-specific citation see also https://doi.org/10.5281/zenodo.2537945. + +module module_sf_psychrolib + !+ Module overview + !+ Contains functions for calculating thermodynamic properties of gas-vapor mixtures + !+ and standard atmosphere suitable for most engineering, physical, and meteorological + !+ applications. + !+ + !+ Most of the functions are an implementation of the formulae found in the + !+ 2017 ASHRAE Handbook - Fundamentals, in both International System (SI), + !+ and Imperial (IP) units. Please refer to the information included in + !+ each function for their respective reference. + !+ + !+ Example + !+ use psychrolib, only: GetTDewPointFromRelHum, SetUnitSystem, SI + !+ ! Set the unit system, for example to SI (can be either 'SI' or 'IP') + !+ call SetUnitSystem(SI) + !+ ! Calculate the dew point temperature for a dry bulb temperature of 25 C and a relative humidity of 80% + !+ print *, GetTDewPointFromRelHum(25.0, 0.80) + !+ 21.3094 + !+ + !+ Copyright + !+ - For the current library implementation + !+ Copyright (c) 2018 D. Thevenard and D. Meyer. + !+ - For equations and coefficients published ASHRAE Handbook — Fundamentals, Chapter 1 + !+ Copyright (c) 2017 ASHRAE Handbook — Fundamentals (https://www.ashrae.org) + !+ + !+ License + !+ MIT (https://github.com/psychrometrics/psychrolib/LICENSE.txt) + !+ + !+ Note from the Authors + !+ We have made every effort to ensure that the code is adequate, however, we make no + !+ representation with respect to its accuracy. Use at your own risk. Should you notice + !+ an error, or if you have a suggestion, please notify us through GitHub at + !+ https://github.com/psychrometrics/psychrolib/issues. + + + implicit none + + private + public :: IP + public :: SI + public :: SetUnitSystem + public :: GetUnitSystem + public :: isIP + public :: GetTRankineFromTFahrenheit + public :: GetTKelvinFromTCelsius + public :: GetTWetBulbFromTDewPoint + public :: GetTWetBulbFromRelHum + public :: GetRelHumFromTDewPoint + public :: GetRelHumFromTWetBulb + public :: GetTDewPointFromRelHum + public :: GetTDewPointFromTWetBulb + public :: GetVapPresFromRelHum + public :: GetRelHumFromVapPres + public :: GetTDewPointFromVapPres + public :: GetVapPresFromTDewPoint + public :: GetTWetBulbFromHumRatio + public :: GetHumRatioFromTWetBulb + public :: GetHumRatioFromRelHum + public :: GetRelHumFromHumRatio + public :: GetHumRatioFromTDewPoint + public :: GetTDewPointFromHumRatio + public :: GetHumRatioFromVapPres + public :: GetVapPresFromHumRatio + public :: GetDryAirEnthalpy + public :: GetDryAirDensity + public :: GetDryAirVolume + public :: GetTDryBulbFromEnthalpyAndHumRatio + public :: GetHumRatioFromEnthalpyAndTDryBulb + public :: GetSatVapPres + public :: GetSatHumRatio + public :: GetSatAirEnthalpy + public :: GetVaporPressureDeficit + public :: GetDegreeOfSaturation + public :: GetMoistAirEnthalpy + public :: GetMoistAirVolume + public :: GetMoistAirDensity + public :: GetStandardAtmPressure + public :: GetStandardAtmTemperature + public :: GetSeaLevelPressure + public :: GetStationPressure + public :: GetSpecificHumFromHumRatio + public :: GetHumRatioFromSpecificHum + public :: CalcPsychrometricsFromTWetBulb + public :: CalcPsychrometricsFromTDewPoint + public :: CalcPsychrometricsFromRelHum + public :: dLnPws_ + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Global constants + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + real, parameter :: R_DA_IP = 53.350 + !+ Universal gas constant for dry air (IP version) in ft lb_Force lb_DryAir⁻¹ R⁻¹. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1. + + real, parameter :: R_DA_SI = 287.042 + !+ Universal gas constant for dry air (SI version) in J kg_DryAir⁻¹ K⁻¹. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1. + + integer, parameter :: IP = 1 + integer, parameter :: SI = 2 + + integer :: PSYCHROLIB_UNITS = 0 ! 0 = undefined. + !+ Unit system to use. + + real :: PSYCHROLIB_TOLERANCE = 1.0 + !+ Tolerance of temperature calculations. + + integer, parameter :: MAX_ITER_COUNT = 100 + !+ Maximum number of iterations before exiting while loops. + + real, parameter :: MIN_HUM_RATIO = 1e-7 + !+ Minimum acceptable humidity ratio used/returned by any functions. + !+ Any value above 0 or below the MIN_HUM_RATIO will be reset to this value. + + real, parameter :: FREEZING_POINT_WATER_IP = 32.0 + !+ float: Freezing point of water in Fahrenheit. + + real, parameter :: FREEZING_POINT_WATER_SI = 0.0 + !+ float: Freezing point of water in Celsius. + + real, parameter :: TRIPLE_POINT_WATER_IP = 32.018 + !+ float: Triple point of water in Fahrenheit. + + real, parameter :: TRIPLE_POINT_WATER_SI = 0.01 + !+ float: Triple point of water in Celsius. + + + contains + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Helper functions + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine SetUnitSystem(UnitSystem) + !+ Set the system of units to use (SI or IP). + !+ Notes: this function *HAS TO BE CALLED* before the library can be used + + integer, intent(in) :: UnitSystem + !+ Units: string indicating the system of units chosen (SI or IP) + + if (.not. (UnitSystem == SI .or. UnitSystem == IP)) then + error stop "The system of units has to be either SI or IP." + end if + + PSYCHROLIB_UNITS = UnitSystem + + ! Define tolerance on temperature calculations + ! The tolerance is the same in IP and SI + if (UnitSystem == IP) then + PSYCHROLIB_TOLERANCE = 0.001 * 9.0 / 5.0 + else + PSYCHROLIB_TOLERANCE = 0.001 + end if + end subroutine SetUnitSystem + + function GetUnitSystem() result(UnitSystem) + !+ Return the system of units in use. + integer :: UnitSystem + UnitSystem = PSYCHROLIB_UNITS + end function GetUnitSystem + + function isIP() + !+ Check whether the system in use is IP or SI + logical :: isIP + if (PSYCHROLIB_UNITS == IP) then + isIP = .true. + else if (PSYCHROLIB_UNITS == SI) then + isIP = .false. + else + error stop "The system of units has not been defined." + end if + end function isIP + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Conversion between temperature units + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetTRankineFromTFahrenheit(TFahrenheit) result(TRankine) + real, intent(in) :: TFahrenheit + !+ Temperature in degree Fahrenheit + real :: TRankine + !+ Temperature in degree Rankine + real :: ZERO_FAHRENHEIT_AS_RANKINE + !+ Zero degree Fahrenheit (°F) expressed as degree Rankine (°R) + + ZERO_FAHRENHEIT_AS_RANKINE = 459.67 + + TRankine = TFahrenheit + ZERO_FAHRENHEIT_AS_RANKINE + end function GetTRankineFromTFahrenheit + + function GetTKelvinFromTCelsius(TCelsius) result(TKelvin) + real, intent(in) :: TCelsius + !+ Temperature in degree Celsius + real :: TKelvin + !+ Tempearatyre in Kelvin + real :: ZERO_CELSIUS_AS_KELVIN + ! Zero degree Fahrenheit (°F) expressed as Kelvin (K) + + ZERO_CELSIUS_AS_KELVIN = 273.15 + + TKelvin = TCelsius + ZERO_CELSIUS_AS_KELVIN + end function GetTKelvinFromTCelsius + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Conversions between dew point, wet bulb, and relative humidity + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetTWetBulbFromTDewPoint(TDryBulb, TDewPoint, Pressure) result(TWetBulb) + !+ Return wet-bulb temperature given dry-bulb temperature, dew-point temperature, and pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + if (TDewPoint > TDryBulb) then + error stop "Error: dew point temperature is above dry bulb temperature" + end if + + HumRatio = GetHumRatioFromTDewPoint(TDewPoint, Pressure) + TWetBulb = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) + end function GetTWetBulbFromTDewPoint + + function GetTWetBulbFromRelHum(TDryBulb, RelHum, Pressure) result(TWetBulb) + !+ Return wet-bulb temperature given dry-bulb temperature, relative humidity, and pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: RelHum + !+ Relative humidity in range [0, 1] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + if (RelHum < 0.0 .or. RelHum > 1.0) then + error stop "Error: relative humidity is outside range [0,1]" + end if + + HumRatio = GetHumRatioFromRelHum(TDryBulb, RelHum, Pressure) + TWetBulb = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) + end function GetTWetBulbFromRelHum + + function GetRelHumFromTDewPoint(TDryBulb, TDewPoint) result(RelHum) + !+ Return relative humidity given dry-bulb temperature and dew-point temperature. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 22 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real :: RelHum + !+ Relative humidity in range [0, 1] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real :: SatVapPres + !+ Vapor pressure of saturated air in Psi [IP] or Pa [SI] + + if (TDewPoint > TDryBulb) then + error stop "Error: dew point temperature is above dry bulb temperature" + end if + + VapPres = GetSatVapPres(TDewPoint) + SatVapPres = GetSatVapPres(TDryBulb) + RelHum = VapPres / SatVapPres + end function GetRelHumFromTDewPoint + + function GetRelHumFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(RelHum) + !+ Return relative humidity given dry-bulb temperature, wet bulb temperature and pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: RelHum + !+ Relative humidity in range [0, 1] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + if (TWetBulb > TDryBulb) then + error stop "Error: wet bulb temperature is above dry bulb temperature" + end if + + HumRatio = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) + RelHum = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure) + end function GetRelHumFromTWetBulb + + function GetTDewPointFromRelHum(TDryBulb, RelHum) result(TDewPoint) + !+ Return dew-point temperature given dry-bulb temperature and relative humidity. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: RelHum + !+ Relative humidity in range [0, 1] + real :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + + if (RelHum < 0.0 .or. RelHum > 1.0) then + error stop "Error: relative humidity is outside range [0,1]" + end if + + VapPres = GetVapPresFromRelHum(TDryBulb, RelHum) + TDewPoint = GetTDewPointFromVapPres(TDryBulb, VapPres) + end function GetTDewPointFromRelHum + + function GetTDewPointFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(TDewPoint) + !+ Return dew-point temperature given dry-bulb temperature, wet-bulb temperature, and pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + if (TWetBulb > TDryBulb) then + error stop "Error: wet bulb temperature is above dry bulb temperature" + end if + + HumRatio = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) + end function GetTDewPointFromTWetBulb + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Conversions between dew point, or relative humidity and vapor pressure + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetVapPresFromRelHum(TDryBulb, RelHum) result(VapPres) + !+ Return partial pressure of water vapor as a function of relative humidity and temperature. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 12, 22 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: RelHum + !+ Relative humidity in range [0, 1] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + + if (RelHum < 0.0 .or. RelHum > 1.0) then + error stop "Error: relative humidity is outside range [0,1]" + end if + + VapPres = RelHum * GetSatVapPres(TDryBulb) + end function GetVapPresFromRelHum + + function GetRelHumFromVapPres(TDryBulb, VapPres) result(RelHum) + !+ Return relative humidity given dry-bulb temperature and vapor pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 12, 22 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real :: RelHum + !+ Relative humidity in range [0, 1] + + if (VapPres < 0.0) then + error stop "Error: partial pressure of water vapor in moist air cannot be negative" + end if + + RelHum = VapPres / GetSatVapPres(TDryBulb) + end function GetRelHumFromVapPres + + function dLnPws_(TDryBulb) result(dLnPws) + !+ Helper function returning the derivative of the natural log of the saturation vapor pressure + !+ as a function of dry-bulb temperature. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: dLnPws + !+ Derivative of natural log of vapor pressure of saturated air in Psi [IP] or Pa [SI] + real :: T + !+ Dry bulb temperature in R [IP] or K [SI] + + if (isIP()) then + + T = GetTRankineFromTFahrenheit(TDryBulb) + + if (TDryBulb <= TRIPLE_POINT_WATER_IP) then + dLnPws = 1.0214165E+04 / T**2 - 5.3765794E-03 + 2 * 1.9202377E-07 * T & + + 2 * 3.5575832E-10 * T**2 - 4 * 9.0344688E-14 * T**3 + 4.1635019 / T + else + dLnPws = 1.0440397E+04 / T**2 - 2.7022355E-02 + 2 * 1.2890360E-05 * T & + - 3 * 2.4780681E-09 * T**2 + 6.5459673 / T + end if + + else + + T = GetTKelvinFromTCelsius(TDryBulb) + + if (TDryBulb <= TRIPLE_POINT_WATER_SI) then + dLnPws = 5.6745359E+03 / T**2 - 9.677843E-03 + 2 * 6.2215701E-07 * T & + + 3 * 2.0747825E-09 * T**2 - 4 * 9.484024E-13 * T**3 + 4.1635019 / T + else + dLnPws = 5.8002206E+03 / T**2 - 4.8640239E-02 + 2 * 4.1764768E-05 * T & + - 3 * 1.4452093E-08 * T**2 + 6.5459673 / T + end if + end if + end function dLnPws_ + + function GetTDewPointFromVapPres(TDryBulb, VapPres) result(TDewPoint) + !+ Return dew-point temperature given dry-bulb temperature and vapor pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn. 5 and 6 + !+ Notes: + !+ The dew point temperature is solved by inverting the equation giving water vapor pressure + !+ at saturation from temperature rather than using the regressions provided + !+ by ASHRAE (eqn. 37 and 38) which are much less accurate and have a + !+ narrower range of validity. + !+ The Newton-Raphson (NR) method is used on the logarithm of water vapour + !+ pressure as a function of temperature, which is a very smooth function + !+ Convergence is usually achieved in 3 to 5 iterations. + !+ TDryBulb is not really needed here, just used for convenience. + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real :: lnVP + !+ Natural logarithm of partial pressure of water vapor pressure in moist air + real :: d_lnVP + !+ Derivative of function, calculated numerically + real :: lnVP_iter + !+ Value of log of vapor water pressure used in NR calculation + real :: TDewPoint_iter + !+ Value of TDewPoint used in NR calculation + real, dimension(2) :: BOUNDS + !+ Valid temperature range in °F [IP] or °C [SI] + integer :: index + !+ Index used in the calculation + + ! Bounds and step size as a function of the system of units + if (isIP()) then + BOUNDS(1) = -148.0 + BOUNDS(2) = 392.0 + else + BOUNDS(1) = -100.0 + BOUNDS(2) = 200.0 + end if + + ! Validity check -- bounds outside which a solution cannot be found + if (VapPres < GetSatVapPres(BOUNDS(1)) .or. VapPres > GetSatVapPres(BOUNDS(2))) then + error stop "Error: partial pressure of water vapor is outside range of validity of equations" + end if + + ! We use NR to approximate the solution. + TDewPoint = TDryBulb + lnVP = log(VapPres) + index = 1 + + do while (.true.) + TDewPoint_iter = TDewPoint ! TDewPoint_iter used in NR calculation + lnVP_iter = log(GetSatVapPres(TDewPoint_iter)) + + ! Derivative of function, calculated analytically + d_lnVP = dLnPws_(TDewPoint_iter) + + ! New estimate, bounded by the search domain defined above + TDewPoint = TDewPoint_iter - (lnVP_iter - lnVP) / d_lnVP + TDewPoint = max(TDewPoint, BOUNDS(1)) + TDewPoint = min(TDewPoint, BOUNDS(2)) + + if (abs(TDewPoint - TDewPoint_iter) <= PSYCHROLIB_TOLERANCE) then + exit + end if + + if (index > MAX_ITER_COUNT) then + error stop "Convergence not reached in GetTDewPointFromVapPres. Stopping." + end if + + index = index + 1 + end do + + TDewPoint = min(TDewPoint, TDryBulb) + end function GetTDewPointFromVapPres + + function GetVapPresFromTDewPoint(TDewPoint) result(VapPres) + !+ Return vapor pressure given dew point temperature. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 36 + + real, intent(in) :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + + VapPres = GetSatVapPres(TDewPoint) + end function GetVapPresFromTDewPoint + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Conversions from wet-bulb temperature, dew-point temperature, or relative humidity to humidity ratio + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) result(TWetBulb) + !+ Return wet-bulb temperature given dry-bulb temperature, humidity ratio, and pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 33 and 35 solved for Tstar + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real :: TDewPoint + !+ TDewPoint : Dew-point temperature in °F [IP] or °C [SI] + real :: TWetBulbSup + !+ Upper value of wet bulb temperature in bissection method (initial guess is from dry bulb temperature) in °F [IP] or °C [SI] + real :: TWetBulbInf + !+ Lower value of wet bulb temperature in bissection method (initial guess is from dew point temperature) in °F [IP] or °C [SI] + real :: Wstar + !+ Humidity ratio at temperature Tstar in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + integer :: index + !+ index used in iteration + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio cannot be negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, BoundedHumRatio, Pressure) + + ! Initial guesses + TWetBulbSup = TDryBulb + TWetBulbInf = TDewPoint + TWetBulb = (TWetBulbInf + TWetBulbSup) / 2.0 + + index = 1 + ! Bisection loop + do while ((TWetBulbSup - TWetBulbInf) > PSYCHROLIB_TOLERANCE) + + ! Compute humidity ratio at temperature Tstar + Wstar = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) + + ! Get new bounds + if (Wstar > BoundedHumRatio) then + TWetBulbSup = TWetBulb + else + TWetBulbInf = TWetBulb + end if + + ! New guess of wet bulb temperature + TWetBulb = (TWetBulbSup + TWetBulbInf) / 2.0 + + if (index > MAX_ITER_COUNT) then + error stop "Convergence not reached in GetTWetBulbFromHumRatio. Stopping." + end if + + index = index + 1 + end do + end function GetTWetBulbFromHumRatio + + function GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) result(HumRatio) + !+ Return humidity ratio given dry-bulb temperature, wet-bulb temperature, and pressure. + !+ References: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 33 and 35 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: Wsstar + !+ Humidity ratio at temperature Tstar in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + if (TWetBulb > TDryBulb) then + error stop "Error: wet bulb temperature is above dry bulb temperature" + end if + + Wsstar = GetSatHumRatio(TWetBulb, Pressure) + + if (isIP()) then + if (TWetBulb >= FREEZING_POINT_WATER_IP) then + HumRatio = ((1093.0 - 0.556 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb)) & + / (1093.0 + 0.444 * TDryBulb - TWetBulb) + else + HumRatio = ((1220.0 - 0.04 * TWetBulb) * Wsstar - 0.240 * (TDryBulb - TWetBulb)) & + / (1220.0 + 0.444 * TDryBulb - 0.48 * TWetBulb) + end if + else + if (TWetBulb >= FREEZING_POINT_WATER_SI) then + HumRatio = ((2501.0 - 2.326 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) & + / (2501.0 + 1.86 * TDryBulb - 4.186 * TWetBulb) + else + HumRatio = ((2830.0 - 0.24 * TWetBulb) * Wsstar - 1.006 * (TDryBulb - TWetBulb)) & + / (2830.0 + 1.86 * TDryBulb - 2.1 * TWetBulb) + end if + end if + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) + end function GetHumRatioFromTWetBulb + + function GetHumRatioFromRelHum(TDryBulb, RelHum, Pressure) result(HumRatio) + !+ Return humidity ratio given dry-bulb temperature, relative humidity, and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: RelHum + !+ Relative humidity in range [0, 1] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + + if (RelHum < 0.0 .or. RelHum > 1.0) then + error stop "Error: relative humidity is outside range [0,1]" + end if + + VapPres = GetVapPresFromRelHum(TDryBulb, RelHum) + HumRatio = GetHumRatioFromVapPres(VapPres, Pressure) + end function GetHumRatioFromRelHum + + function GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure) result(RelHum) + !+ Return relative humidity given dry-bulb temperature, humidity ratio, and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: RelHum + !+ Relative humidity in range [0, 1] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio cannot be negative" + end if + + VapPres = GetVapPresFromHumRatio(HumRatio, Pressure) + RelHum = GetRelHumFromVapPres(TDryBulb, VapPres) + end function GetRelHumFromHumRatio + + function GetHumRatioFromTDewPoint(TDewPoint, Pressure) result(HumRatio) + !+ Return humidity ratio given dew-point temperature and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + + VapPres = GetSatVapPres(TDewPoint) + HumRatio = GetHumRatioFromVapPres(VapPres, Pressure) + end function GetHumRatioFromTDewPoint + + function GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) result(TDewPoint) + !+ Return dew-point temperature given dry-bulb temperature, humidity ratio, and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio cannot be negative" + end if + + VapPres = GetVapPresFromHumRatio(HumRatio, Pressure) + TDewPoint = GetTDewPointFromVapPres(TDryBulb, VapPres) + end function GetTDewPointFromHumRatio + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Conversions between humidity ratio and vapor pressure + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetHumRatioFromVapPres(VapPres, Pressure) result(HumRatio) + !+ Return humidity ratio given water vapor pressure and atmospheric pressure. + !+ Reference: + !+ ASHRAE Fundamentals (2005) ch. 6 eqn. 22; + !+ ASHRAE Fundamentals (2009) ch. 1 eqn. 22. + + real, intent(in) :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + if (VapPres < 0.0) then + error stop "Error: partial pressure of water vapor in moist air cannot be negative" + end if + + HumRatio = 0.621945 * VapPres / (Pressure-VapPres) + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) + end function GetHumRatioFromVapPres + + function GetVapPresFromHumRatio(HumRatio, Pressure) result(VapPres) + !+ Return vapor pressure given humidity ratio and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 20 solved for pw + + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio is negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + VapPres = Pressure * BoundedHumRatio / (0.621945 + BoundedHumRatio) + end function GetVapPresFromHumRatio + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Conversions between humidity ratio and specific humidity + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetSpecificHumFromHumRatio(HumRatio) result(SpecificHum) + !+ Return the specific humidity from humidity ratio (aka mixing ratio). + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 9b + + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Dry_Air⁻¹ [IP] or kg_H₂O kg_Dry_Air⁻¹ [SI] + real :: SpecificHum + !+ Specific humidity in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + + if (HumRatio < 0.0) then + print *, "Error: humidity ratio cannot be negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + SpecificHum = BoundedHumRatio / (1.0 + BoundedHumRatio) + end function GetSpecificHumFromHumRatio + + function GetHumRatioFromSpecificHum(SpecificHum) result(HumRatio) + !+ Return the humidity ratio (aka mixing ratio) from specific humidity. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 9b (solved for humidity ratio) + + real, intent(in) :: SpecificHum + !+ Specific humidity in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Dry_Air⁻¹ [IP] or kg_H₂O kg_Dry_Air⁻¹ [SI] + + if (SpecificHum < 0.0 .or. SpecificHum >= 1.0) then + print *, "Error: specific humidity is outside range [0, 1[" + end if + + HumRatio = SpecificHum / (1.0 - SpecificHum) + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) + end function GetHumRatioFromSpecificHum + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Dry Air Calculations + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetDryAirEnthalpy(TDryBulb) result(DryAirEnthalpy) + !+ Return dry-air enthalpy given dry-bulb temperature. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 28 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: DryAirEnthalpy + !+ Dry air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI] + + if (isIP()) then + DryAirEnthalpy = 0.240 * TDryBulb + else + DryAirEnthalpy = 1006 * TDryBulb + end if + end function GetDryAirEnthalpy + + function GetDryAirDensity(TDryBulb, Pressure) result(DryAirDensity) + !+ Return dry-air density given dry-bulb temperature and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + !+ Notes: + !+ Eqn 14 for the perfect gas relationship for dry air. + !+ Eqn 1 for the universal gas constant. + !+ The factor 144 in IP is for the conversion of Psi = lb in⁻² to lb ft⁻². + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: DryAirDensity + !+ Dry air density in lb ft⁻³ [IP] or kg m⁻³ [SI] + + if (isIP()) then + DryAirDensity = (144 * Pressure) / R_DA_IP / GetTRankineFromTFahrenheit(TDryBulb) + else + DryAirDensity = Pressure / R_DA_SI / GetTKelvinFromTCelsius(TDryBulb) + end if + end function GetDryAirDensity + + function GetDryAirVolume(TDryBulb, Pressure) result(DryAirVolume) + !+ Return dry-air volume given dry-bulb temperature and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + !+ Notes: + !+ Eqn 14 for the perfect gas relationship for dry air. + !+ Eqn 1 for the universal gas constant. + !+ The factor 144 in IP is for the conversion of Psi = lb in⁻² to lb ft⁻². + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: DryAirVolume + !+ Dry air volume in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI] + + if (isIP()) then + DryAirVolume = GetTRankineFromTFahrenheit(TDryBulb) * R_DA_IP / (144 * Pressure) + else + DryAirVolume = GetTKelvinFromTCelsius(TDryBulb) * R_DA_SI / Pressure + end if + end function GetDryAirVolume + + function GetTDryBulbFromEnthalpyAndHumRatio(MoistAirEnthalpy, HumRatio) result(TDryBulb) + !+ Return dry bulb temperature from enthalpy and humidity ratio. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 30 + !+ Notes: + !+ Based on the `GetMoistAirEnthalpy` function, rearranged for temperature. + + real, intent(in) :: MoistAirEnthalpy + !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio is negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + if (isIP()) then + TDryBulb = (MoistAirEnthalpy - 1061.0 * BoundedHumRatio) / (0.240 + 0.444 * BoundedHumRatio) + else + TDryBulb = (MoistAirEnthalpy / 1000.0 - 2501.0 * BoundedHumRatio) / (1.006 + 1.86 * BoundedHumRatio) + end if + end function GetTDryBulbFromEnthalpyAndHumRatio + + function GetHumRatioFromEnthalpyAndTDryBulb(MoistAirEnthalpy, TDryBulb) result(HumRatio) + !+ Return humidity ratio from enthalpy and dry-bulb temperature. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 30 + !+ Notes: + !+ Based on the `GetMoistAirEnthalpy` function, rearranged for humidity ratio. + + real, intent(in) :: MoistAirEnthalpy + !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + if (isIP()) then + HumRatio = (MoistAirEnthalpy - 0.240 * TDryBulb) / (1061.0 + 0.444 * TDryBulb) + else + HumRatio = (MoistAirEnthalpy / 1000.0 - 1.006 * TDryBulb) / (2501.0 + 1.86 * TDryBulb) + end if + + ! Validity check. + HumRatio = max(HumRatio, MIN_HUM_RATIO) + end function GetHumRatioFromEnthalpyAndTDryBulb + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Saturated Air Calculations + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetSatVapPres(TDryBulb) result(SatVapPres) + !+ Return saturation vapor pressure given dry-bulb temperature. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 5 + !+ Important note: the ASHRAE formulae are defined above and below the freezing point but have + !+ a discontinuity at the freezing point. This is a small inaccuracy on ASHRAE's part: the formulae + !+ should be defined above and below the triple point of water (not the feezing point) in which case + !+ the discontinuity vanishes. It is essential to use the triple point of water otherwise function + !+ GetTDewPointFromVapPres, which inverts the present function, does not converge properly around + !+ the freezing point. + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: SatVapPres + !+ Vapor pressure of saturated air in Psi [IP] or Pa [SI] + real :: LnPws + !+ Log of Vapor Pressure of saturated air (dimensionless) + real :: T + !+ Dry bulb temperature in R [IP] or K [SI] + + if (isIP()) then + if (TDryBulb < -148.0 .or. TDryBulb > 392.0) then + error stop "Error: dry bulb temperature must be in range [-148, 392]°F" + end if + + T = GetTRankineFromTFahrenheit(TDryBulb) + + if (TDryBulb <= TRIPLE_POINT_WATER_IP) then + LnPws = (-1.0214165E+04 / T - 4.8932428 - 5.3765794E-03 * T + 1.9202377E-07 * T**2 & + + 3.5575832E-10 * T**3 - 9.0344688E-14 * T**4 + 4.1635019 * log(T)) + else + LnPws = -1.0440397E+04 / T - 1.1294650E+01 - 2.7022355E-02* T + 1.2890360E-05 * T**2 & + - 2.4780681E-09 * T**3 + 6.5459673 * log(T) + end if + + else + if (TDryBulb < -100.0 .or. TDryBulb > 200.0) then + error stop "Error: dry bulb temperature must be in range [-100, 200]°C" + end if + + T = GetTKelvinFromTCelsius(TDryBulb) + + if (TDryBulb <= TRIPLE_POINT_WATER_SI) then + LnPws = -5.6745359E+03 / T + 6.3925247 - 9.677843E-03 * T + 6.2215701E-07 * T**2 & + + 2.0747825E-09 * T**3 - 9.484024E-13 * T**4 + 4.1635019 * log(T) + else + LnPws = -5.8002206E+03 / T + 1.3914993 - 4.8640239E-02 * T + 4.1764768E-05 * T**2 & + - 1.4452093E-08 * T**3 + 6.5459673 * log(T) + end if + end if + + SatVapPres = exp(LnPws) + end function GetSatVapPres + + function GetSatHumRatio(TDryBulb, Pressure) result(SatHumRatio) + !+ Return humidity ratio of saturated air given dry-bulb temperature and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 36, solved for W + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: SatHumRatio + !+ Humidity ratio of saturated air in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: SatVaporPres + !+ Vapor pressure of saturated air in in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + + SatVaporPres = GetSatVapPres(TDryBulb) + SatHumRatio = 0.621945 * SatVaporPres / (Pressure-SatVaporPres) + + ! Validity check. + SatHumRatio = max(SatHumRatio, MIN_HUM_RATIO) + end function GetSatHumRatio + + function GetSatAirEnthalpy(TDryBulb, Pressure) result(SatAirEnthalpy) + !+ Return saturated air enthalpy given dry-bulb temperature and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: SatAirEnthalpy + !+ Saturated air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI] + + SatAirEnthalpy = GetMoistAirEnthalpy(TDryBulb, GetSatHumRatio(TDryBulb, Pressure)) + end function GetSatAirEnthalpy + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Moist Air Calculations + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetVaporPressureDeficit(TDryBulb, HumRatio, Pressure) result(VaporPressureDeficit) + !+ Return Vapor pressure deficit given dry-bulb temperature, humidity ratio, and pressure. + !+ Reference: + !+ Oke (1987) eqn 2.13a + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: VaporPressureDeficit + !+ Vapor pressure deficit in Psi [IP] or Pa [SI] + real :: RelHum + !+ Relative humidity in range [0, 1] + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio is negative" + end if + + RelHum = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure) + VaporPressureDeficit = GetSatVapPres(TDryBulb) * (1.0 - RelHum) + end function GetVaporPressureDeficit + + function GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure) result(DegreeOfSaturation) + !+ Return the degree of saturation (i.e humidity ratio of the air / humidity ratio of the air at saturation + !+ at the same temperature and pressure) given dry-bulb temperature, humidity ratio, and atmospheric pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2009) ch. 1 eqn 12 + !+ Notes: + !+ This definition is absent from the 2017 Handbook. Using 2009 version instead. + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: DegreeOfSaturation + !+ Degree of saturation in arbitrary unit + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio is negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + DegreeOfSaturation = BoundedHumRatio / GetSatHumRatio(TDryBulb, Pressure) + end function GetDegreeOfSaturation + + function GetMoistAirEnthalpy(TDryBulb, HumRatio) result(MoistAirEnthalpy) + !+ Return moist air enthalpy given dry-bulb temperature and humidity ratio. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 30 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real :: MoistAirEnthalpy + !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio is negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + if (isIP()) then + MoistAirEnthalpy = 0.240 * TDryBulb + BoundedHumRatio * (1061.0 + 0.444 * TDryBulb) + else + MoistAirEnthalpy = (1.006 * TDryBulb + BoundedHumRatio * (2501.0 + 1.86 * TDryBulb)) * 1000.0 + end if + end function GetMoistAirEnthalpy + + function GetMoistAirVolume(TDryBulb, HumRatio, Pressure) result(MoistAirVolume) + !+ Return moist air specific volume given dry-bulb temperature, humidity ratio, and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 26 + !+ Notes: + !+ In IP units, R_DA_IP / 144 equals 0.370486 which is the coefficient appearing in eqn 26 + !+ The factor 144 is for the conversion of Psi = lb in⁻² to lb ft⁻². + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: MoistAirVolume + !+ Specific volume of moist air in ft³ lb⁻¹ of dry air [IP] or in m³ kg⁻¹ of dry air [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio is negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + if (isIP()) then + MoistAirVolume = R_DA_IP * GetTRankineFromTFahrenheit(TDryBulb) * (1.0 + 1.607858 * BoundedHumRatio) / (144.0 * Pressure) + else + MoistAirVolume = R_DA_SI * GetTKelvinFromTCelsius(TDryBulb) * (1.0 + 1.607858 * BoundedHumRatio) / Pressure + end if + end function GetMoistAirVolume + + function GetMoistAirDensity(TDryBulb, HumRatio, Pressure) result(MoistAirDensity) + !+ Return moist air density given humidity ratio, dry bulb temperature, and pressure. + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 11 + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real :: MoistAirDensity + !+ Moist air density in lb ft⁻³ [IP] or kg m⁻³ [SI] + real :: BoundedHumRatio + !+ Humidity ratio bounded to MIN_HUM_RATIO + + if (HumRatio < 0.0) then + error stop "Error: humidity ratio is negative" + end if + BoundedHumRatio = max(HumRatio, MIN_HUM_RATIO) + + MoistAirDensity = (1.0 + BoundedHumRatio) / GetMoistAirVolume(TDryBulb, BoundedHumRatio, Pressure) + end function GetMoistAirDensity + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Standard atmosphere + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + function GetStandardAtmPressure(Altitude) result(StandardAtmPressure) + !+ Return standard atmosphere barometric pressure, given the elevation (altitude). + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 3 + + real, intent(in) :: Altitude + !+ Altitude in ft [IP] or m [SI] + real :: StandardAtmPressure + !+ Standard atmosphere barometric pressure in Psi [IP] or Pa [SI] + + if (isIP()) then + StandardAtmPressure = 14.696 * (1.0 - 6.8754e-06 * Altitude)**5.2559 + else + StandardAtmPressure = 101325 * (1 - 2.25577e-05 * Altitude)**5.2559 + end if + end function GetStandardAtmPressure + + function GetStandardAtmTemperature(Altitude) result(StandardAtmTemperature) + !+ Return standard atmosphere temperature, given the elevation (altitude). + !+ Reference: + !+ ASHRAE Handbook - Fundamentals (2017) ch. 1 eqn 4 + + real, intent(in) :: Altitude + !+ Altitude in ft [IP] or m [SI] + real :: StandardAtmTemperature + !+ Standard atmosphere dry-bulb temperature in °F [IP] or °C [SI] + + if (isIP()) then + StandardAtmTemperature = 59.0 - 0.00356620 * Altitude + else + StandardAtmTemperature = 15.0 - 0.0065 * Altitude + end if + end function GetStandardAtmTemperature + + function GetSeaLevelPressure(StnPressure, Altitude, TDryBulb) result(SeaLevelPressure) + !+ Return sea level pressure given dry-bulb temperature, altitude above sea level and pressure. + !+ Reference: + !+ Hess SL, Introduction to theoretical meteorology, Holt Rinehart and Winston, NY 1959, + !+ ch. 6.5; Stull RB, Meteorology for scientists and engineers, 2nd edition, + !+ Brooks/Cole 2000, ch. 1. + !+ Notes: + !+ The standard procedure for the US is to use for TDryBulb the average + !+ of the current station temperature and the station temperature from 12 hours ago. + + real, intent(in) :: StnPressure + !+ Observed station pressure in Psi [IP] or Pa [SI] + real, intent(in) :: Altitude + !+ Altitude in ft [IP] or m [SI] + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: SeaLevelPressure + !+ Sea level barometric pressure in Psi [IP] or Pa [SI] + real :: TColumn + !+ Average temperature in column of air in R [IP] or K [SI] + real :: H + !+ scale height (dimensionless) + + if (isIP()) then + ! Calculate average temperature in column of air, assuming a lapse rate + ! of 3.6 °F/1000ft + TColumn = TDryBulb + 0.0036 * Altitude / 2.0 + + ! Determine the scale height + H = 53.351 * GetTRankineFromTFahrenheit(TColumn) + else + ! Calculate average temperature in column of air, assuming a lapse rate + ! of 6.5 °C/km + TColumn = TDryBulb + 0.0065 * Altitude / 2.0 + + ! Determine the scale height + H = 287.055 * GetTKelvinFromTCelsius(TColumn) / 9.807 + end if + + ! Calculate the sea level pressure + SeaLevelPressure = StnPressure * exp(Altitude / H) + end function GetSeaLevelPressure + + function GetStationPressure(SeaLevelPressure, Altitude, TDryBulb) result(StationPressure) + !+ Return station pressure from sea level pressure. + !+ Reference: + !+ See 'GetSeaLevelPressure' + !+ Notes: + !+ This function is just the inverse of 'GetSeaLevelPressure'. + + real, intent(in) :: SeaLevelPressure + !+ Sea level barometric pressure in Psi [IP] or Pa [SI] + real, intent(in) :: Altitude + !+ Altitude in ft [IP] or m [SI] + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real :: StationPressure + !+ Station pressure in Psi [IP] or Pa [SI] + + StationPressure = SeaLevelPressure / GetSeaLevelPressure(1.0, Altitude, TDryBulb) + end function GetStationPressure + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Functions to set all psychrometric values + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + subroutine CalcPsychrometricsFromTWetBulb(TDryBulb, & + TWetBulb, & + Pressure, & + HumRatio, & + TDewPoint, & + RelHum, & + VapPres, & + MoistAirEnthalpy, & + MoistAirVolume, & + DegreeOfSaturation) + + !+ Utility function to calculate humidity ratio, dew-point temperature, relative humidity, + !+ vapour pressure, moist air enthalpy, moist air volume, and degree of saturation of air given + !+ dry-bulb temperature, wet-bulb temperature, and pressure. + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real, intent(out) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(out) :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real, intent(out) :: RelHum + !+ Relative humidity in range [0, 1] + real, intent(out) :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real, intent(out) :: MoistAirEnthalpy + !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI] + real, intent(out) :: MoistAirVolume + !+ Specific volume of moist air in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI] + real, intent(out) :: DegreeOfSaturation + !+ Degree of saturation [unitless] + + HumRatio = GetHumRatioFromTWetBulb(TDryBulb, TWetBulb, Pressure) + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) + RelHum = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure) + VapPres = GetVapPresFromHumRatio(HumRatio, Pressure) + MoistAirEnthalpy = GetMoistAirEnthalpy(TDryBulb, HumRatio) + MoistAirVolume = GetMoistAirVolume(TDryBulb, HumRatio, Pressure) + DegreeOfSaturation = GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure) + end subroutine CalcPsychrometricsFromTWetBulb + + subroutine CalcPsychrometricsFromTDewPoint(TDryBulb, & + TDewPoint, & + Pressure, & + HumRatio, & + TWetBulb, & + RelHum, & + VapPres, & + MoistAirEnthalpy, & + MoistAirVolume, & + DegreeOfSaturation) + + !+ Utility function to calculate humidity ratio, wet-bulb temperature, relative humidity, + !+ vapour pressure, moist air enthalpy, moist air volume, and degree of saturation of air given + !+ dry-bulb temperature, dew-point temperature, and pressure. + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real, intent(out) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(out) :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real, intent(out) :: RelHum + !+ Relative humidity in range [0, 1] + real, intent(out) :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real, intent(out) :: MoistAirEnthalpy + !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI] + real, intent(out) :: MoistAirVolume + !+ Specific volume of moist air in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI] + real, intent(out) :: DegreeOfSaturation + !+ Degree of saturation [unitless] + + HumRatio = GetHumRatioFromTDewPoint(TDewPoint, Pressure) + TWetBulb = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) + RelHum = GetRelHumFromHumRatio(TDryBulb, HumRatio, Pressure) + VapPres = GetVapPresFromHumRatio(HumRatio, Pressure) + MoistAirEnthalpy = GetMoistAirEnthalpy(TDryBulb, HumRatio) + MoistAirVolume = GetMoistAirVolume(TDryBulb, HumRatio, Pressure) + DegreeOfSaturation = GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure) + end subroutine CalcPsychrometricsFromTDewPoint + + subroutine CalcPsychrometricsFromRelHum(TDryBulb, & + RelHum, & + Pressure, & + HumRatio, & + TWetBulb, & + TDewPoint, & + VapPres, & + MoistAirEnthalpy, & + MoistAirVolume, & + DegreeOfSaturation) + + !+ Utility function to calculate humidity ratio, wet-bulb temperature, dew-point temperature, + !+ vapour pressure, moist air enthalpy, moist air volume, and degree of saturation of air given + !+ dry-bulb temperature, relative humidity and pressure. + + real, intent(in) :: TDryBulb + !+ Dry-bulb temperature in °F [IP] or °C [SI] + real, intent(in) :: RelHum + !+ Relative humidity in range [0, 1] + real, intent(in) :: Pressure + !+ Atmospheric pressure in Psi [IP] or Pa [SI] + real, intent(out) :: HumRatio + !+ Humidity ratio in lb_H₂O lb_Air⁻¹ [IP] or kg_H₂O kg_Air⁻¹ [SI] + real, intent(out) :: TWetBulb + !+ Wet-bulb temperature in °F [IP] or °C [SI] + real, intent(out) :: TDewPoint + !+ Dew-point temperature in °F [IP] or °C [SI] + real, intent(out) :: VapPres + !+ Partial pressure of water vapor in moist air in Psi [IP] or Pa [SI] + real, intent(out) :: MoistAirEnthalpy + !+ Moist air enthalpy in Btu lb⁻¹ [IP] or J kg⁻¹ [SI] + real, intent(out) :: MoistAirVolume + !+ Specific volume of moist air in ft³ lb⁻¹ [IP] or in m³ kg⁻¹ [SI] + real, intent(out) :: DegreeOfSaturation + !+ Degree of saturation [unitless] + + HumRatio = GetHumRatioFromRelHum(TDryBulb, RelHum, Pressure) + TWetBulb = GetTWetBulbFromHumRatio(TDryBulb, HumRatio, Pressure) + TDewPoint = GetTDewPointFromHumRatio(TDryBulb, HumRatio, Pressure) + VapPres = GetVapPresFromHumRatio(HumRatio, Pressure) + MoistAirEnthalpy = GetMoistAirEnthalpy(TDryBulb, HumRatio) + MoistAirVolume = GetMoistAirVolume(TDryBulb, HumRatio, Pressure) + DegreeOfSaturation = GetDegreeOfSaturation(TDryBulb, HumRatio, Pressure) + end subroutine CalcPsychrometricsFromRelHum + + +end module module_sf_psychrolib \ No newline at end of file diff --git a/phys/module_sf_teb.F b/phys/module_sf_teb.F new file mode 100644 index 0000000000..8ad9ea1aa5 --- /dev/null +++ b/phys/module_sf_teb.F @@ -0,0 +1,639 @@ +! WRF-TEB (https://github.com/WRF-CMake/wrf-teb). +! Copyright 2018 D. Meyer. Licensed under the MIT License. + +MODULE MODULE_SF_TEB + + IMPLICIT NONE + + PRIVATE + PUBLIC :: TEB + + CONTAINS + + SUBROUTINE TEB ( & + ims ,ime, jms, jme, kms, kme, & + its ,ite, jts, jte, kts, kte, & + TEB_TEST_INTEGRATION, & + TEB_NUM_ROOF_LAYERS, TEB_NUM_WALL_LAYERS, TEB_NUM_ROAD_LAYERS, & + TEB_NUM_SNOW_LAYERS, TEB_NUM_FLOOR_LAYERS, & + XLAT, XLONG, & + JULDAY, JULYR, GMT, XTIME, XTSTEP_SURF, & + U_PHY, V_PHY, & + DZ8W, P8W3D, T3D, QV3D, RHO, & + SWDOWN, SW_DIFFUSE_FRAC, GLW, & + RAINBL, SR, FRAC_URB2D, UTYPE_URB2D, & + IVGTYP, ISURBAN, & + CQS2, CHS2, & + ! Input from noah + ALB_RURAL, EMISSI_RURAL, HFX_RURAL, QFX_RURAL, & + LH_RURAL, GRDFLX_RURAL, & + TEB_TSK_RURAL, QSFC_RURAL, & + ! End input from noah + TEB_TI_BLD, & + TEB_T_CANYON, TEB_Q_CANYON, & + TEB_WS_ROOF, TEB_WS_ROAD, & + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & + TEB_T_WIN1, TEB_T_WIN2, & + TEB_AUX_MAX, TEB_THER_PRODC_DAY, & + TEB_QI_BLD, & + TEB_T_FLOOR, TEB_T_MASS, & + TEB_T_ROAD, TEB_T_ROOF, & + TEB_T_WALL_A, TEB_T_WALL_B, & + TEB_HVAC_COOL, TEB_HVAC_HEAT, & + TEB_THER_PROD_PANEL, TEB_PHOT_PROD_PANEL, & + RN_URB2D, U10_URB2D, V10_URB2D, & + ! Output Noah + TEB + ALBEDO, HFX, QFX, LH, GRDFLX, TSK, EMISS, QSFC, UST, & + TH2_URB2D, Q2_URB2D, & + ! End output Noah + TEB + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + TEB_INPUT_TA, & + TEB_INPUT_PS, & + TEB_INPUT_QA, & + TEB_INPUT_WIND, & + TEB_INPUT_DIR, & + TEB_INPUT_DIR_SW, & + TEB_INPUT_SCA_SW, & + TEB_INPUT_LW, & + TEB_INPUT_RAIN, & + TEB_INPUT_SNOW, & + TEB_INPUT_DIR_CO2) + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + + USE module_sf_noahlsm, only : LOW_DENSITY_RESIDENTIAL, HIGH_DENSITY_RESIDENTIAL, HIGH_INTENSITY_INDUSTRIAL + + USE module_sf_urban, only: TEB_SOLAR_PANEL_TBL, & + TEB_USE_NOAH, TEB_USE_ZGARDEN_FROM_TBL, & + TEB_LGARDEN_TBL, TEB_GREENROOF_TBL, TEB_HROAD_DIR_TBL, TEB_ZROAD_DIR_TBL, TEB_WALL_OPT_TBL, & + TEB_COOL_COIL_TBL, TEB_HEAT_COIL_TBL, TEB_F_WATER_COND_TBL, TEB_HNATVENT_TBL, TEB_ZNATVENT_TBL, & + TEB_F_WASTE_CAN_TBL, TEB_QIN_TBL, TEB_QIN_FRAD_TBL, TEB_QIN_FLAT_TBL, TEB_GR_TBL, TEB_EFF_HEAT_TBL, & + TEB_INF_TBL, TEB_TCOOL_TARGET_TBL, TEB_THEAT_TARGET_TBL, TEB_HR_TARGET_TBL, TEB_V_VENT_TBL, & + TEB_CAP_SYS_HEAT_TBL, TEB_CAP_SYS_RAT_TBL, TEB_T_ADP_TBL, TEB_M_SYS_RAT_TBL, TEB_COP_RAT_TBL, & + TEB_SHGC_TBL, TEB_SHGC_SH_TBL, TEB_U_WIN_TBL, TEB_LSHADE_TBL, TEB_ZSHADE_TBL, TEB_FLOOR_HEIGHT_TBL, & + TEB_CH_BEM_TBL, TEB_ROUGH_ROOF_TBL, TEB_ROUGH_WALL_TBL, TEB_PAR_RD_IRRIG_TBL, TEB_RD_START_MONTH_TBL, & + TEB_RD_END_MONTH_TBL, TEB_RD_START_HOUR_TBL, TEB_RD_END_HOUR_TBL, TEB_RD_24H_IRRIG_TBL, TEB_EMIS_PANEL_TBL, & + TEB_ALB_PANEL_TBL, TEB_EFF_PANEL_TBL, TEB_FRAC_PANEL_TBL, TEB_RESIDENTIAL_TBL, TEB_DT_RES_TBL, TEB_DT_OFF_TBL, & + TEB_CBEM_TBL, TEB_Z0H_TBL, TEB_OPT_ZOH_TBL, TEB_SNOW_ROOF_TBL, TEB_SNOW_SCH_ROOF_TBL, TEB_SNOW_ROAD_TBL, & + TEB_SNOW_SCH_ROAD_TBL, TEB_ZGARDEN_TBL, TEB_FRAC_GR_TBL, TEB_Z0_TOWN_TBL, TEB_BLD_HEIGHT_TBL, TEB_FRAC_BLD_TBL, & + TEB_WALL_O_TBL, TEB_H_TRAFFIC_TBL, TEB_LE_TRAFFIC_TBL, TEB_H_INDUSTRY_TBL, TEB_LE_INDUSTRY_TBL, & + TEB_ALB_ROOF_TBL, TEB_EMIS_ROOF_TBL, TEB_ALB_ROAD_TBL, TEB_EMIS_ROAD_TBL, & + TEB_ALB_WALL_TBL, TEB_EMIS_WALL_TBL, & + TEB_HC_ROOF_TBL, TEB_TC_ROOF_TBL, TEB_D_ROOF_TBL, & + TEB_HC_ROAD_TBL, TEB_TC_ROAD_TBL, TEB_D_ROAD_TBL, & + TEB_HC_WALL_TBL, TEB_TC_WALL_TBL, TEB_D_WALL_TBL, & + TEB_HC_FLOOR_TBL, TEB_TC_FLOOR_TBL, TEB_D_FLOOR_TBL + + USE module_ra_gfdleta, only: cal_mon_day + + USE modd_wrf_teb_driver, only: TEB_DRIVER + USE module_sf_psychrolib, only: GetSpecificHumFromHumRatio, GetHumRatioFromSpecificHum + + ! TEB imports + USE MODD_CSTS, ONLY : XCPD, & ! Cpd (dry air) + XRV, & ! gaz constant for vapor + XRD, & ! Gaz constant for dry air + XG ! Gravity constant + USE MODD_SURF_PAR, ONLY : XUNDEF ! "no-data" value + USE MODI_ADD_FORECAST_TO_DATE_SURF + + IMPLICIT NONE + + INTEGER, INTENT(IN) :: TEB_TEST_INTEGRATION + + ! array dimension sizes + INTEGER, INTENT(IN) :: ims, ime, jms, jme, kms, kme + INTEGER, INTENT(IN) :: its, ite, jts, jte, kts, kte + INTEGER, INTENT(IN) :: TEB_NUM_ROOF_LAYERS + INTEGER, INTENT(IN) :: TEB_NUM_WALL_LAYERS + INTEGER, INTENT(IN) :: TEB_NUM_ROAD_LAYERS + INTEGER, INTENT(IN) :: TEB_NUM_SNOW_LAYERS + INTEGER, INTENT(IN) :: TEB_NUM_FLOOR_LAYERS + + ! spatial coordinates + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: XLAT + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: XLONG + + ! temporal coordinates + INTEGER, INTENT(IN) :: JULDAY, JULYR ! date at simulation start + REAL, INTENT(IN) :: GMT ! hours since midnight at simulation start + REAL, INTENT(IN) :: XTIME ! minutes since simulation start + REAL, INTENT(IN) :: XTSTEP_SURF ! duration of time step (seconds) + + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: U_PHY ! zonal speed + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: V_PHY ! meridional speed + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: DZ8W ! thickness of each atmospheric model layer + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: P8W3D ! 3D pressure (Pa) This is the pressure at full levels (Pa) + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: T3D ! 3D temperature (K) + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: QV3D ! 3D mixing ratio + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN) :: RHO ! density of the lowest level (kg/m^3) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: SWDOWN ! downward short wave flux at ground surface (W/m^2) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: SW_DIFFUSE_FRAC ! diffuse fraction of surface shortwave irradiance + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: GLW ! downward long wave flux at ground surface (W/m^2) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: RAINBL ! in mm (Accumulation between PBL calls) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: SR ! fraction of frozen precipitation + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: FRAC_URB2D ! urban fraction + INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: UTYPE_URB2D ! urban type + INTEGER, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: IVGTYP ! dominant vegetation category + INTEGER, INTENT(IN) :: ISURBAN ! vegetation category for urban + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: CQS2 ! 2m surface exchange coefficient for moisture + REAL, DIMENSION(ims:ime, jms:jme), INTENT(IN) :: CHS2 ! 2m surface exchange coefficient for heat + ! input from Noah + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: ALB_RURAL + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: EMISSI_RURAL + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: HFX_RURAL + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: QFX_RURAL + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: LH_RURAL + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: GRDFLX_RURAL + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: TEB_TSK_RURAL + REAL, DIMENSION(its:ite, jts:jte), INTENT(IN) :: QSFC_RURAL ! mixing ratio from Noah + ! end input from Noah + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TI_BLD ! indoor air temperature + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_CANYON ! air canyon temperature + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_Q_CANYON ! canyon air humidity ratio [kg kg-1] + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROOF ! roof water content (kg/m2) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROAD ! road water content (kg/m2) + REAL, DIMENSION(ims:ime, TEB_NUM_SNOW_LAYERS, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROOF ! snow layers reservoir + REAL, DIMENSION(ims:ime, TEB_NUM_SNOW_LAYERS, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROOF ! snow layers temperature + REAL, DIMENSION(ims:ime, TEB_NUM_SNOW_LAYERS, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROOF ! snow layers density + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROOF ! snow albedo + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROOF ! snow emissivity + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROOF ! snow surface temperature + REAL, DIMENSION(ims:ime, TEB_NUM_SNOW_LAYERS, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROAD ! snow layers reservoir + REAL, DIMENSION(ims:ime, TEB_NUM_SNOW_LAYERS, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROAD ! snow layers temperature + REAL, DIMENSION(ims:ime, TEB_NUM_SNOW_LAYERS, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROAD ! snow layers density + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROAD ! snow albedo + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROAD ! snow emissivity + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROAD ! snow surface temperature + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN1 ! outdoor window temperature [K] + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN2 ! Indoor window temperature [K] + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_AUX_MAX ! Auxiliar variable for autosize calcs + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_THER_PRODC_DAY ! Present day integrated thermal production of energy (J/m2 panel) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_QI_BLD ! Indoor air specific humidity [kg kg-1] + REAL, DIMENSION(ims:ime, TEB_NUM_FLOOR_LAYERS, jms:jme), INTENT(INOUT) :: TEB_T_FLOOR ! Floor layers temperatures [K] + REAL, DIMENSION(ims:ime, TEB_NUM_FLOOR_LAYERS, jms:jme), INTENT(INOUT) :: TEB_T_MASS ! Internal mass layers temperatures [K] + REAL, DIMENSION(ims:ime, TEB_NUM_ROAD_LAYERS, jms:jme), INTENT(INOUT) :: TEB_T_ROAD ! road layers temperatures + REAL, DIMENSION(ims:ime, TEB_NUM_ROOF_LAYERS, jms:jme), INTENT(INOUT) :: TEB_T_ROOF ! roof layers temperatures + REAL, DIMENSION(ims:ime, TEB_NUM_WALL_LAYERS, jms:jme), INTENT(INOUT) :: TEB_T_WALL_A ! wall layers temperatures (wall 'A') + REAL, DIMENSION(ims:ime, TEB_NUM_WALL_LAYERS, jms:jme), INTENT(INOUT) :: TEB_T_WALL_B ! wall layers temperatures (wall 'B') + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_HVAC_COOL ! Energy consumption of the cooling system [W m-2(bld)] + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_HVAC_HEAT ! Energy consumption of the heating system [W m-2(bld)] + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_THER_PROD_PANEL ! Thermal energy production of solar panel on roofs [W m-2(panel)] + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_PHOT_PROD_PANEL ! Photovoltaic energy production of solar panel on roofs [W m-2(panel)] + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: RN_URB2D + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: U10_URB2D + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: V10_URB2D + ! output Noah + TEB + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: ALBEDO + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: HFX ! upward heat flux at the surface (W m-2) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: QFX ! upward moisture flux at the surface (kg m-2 s-1) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: LH ! latent heat flux at the surface (W m-2) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: GRDFLX ! ground heat flux (W m-2) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TSK ! surface skin temperature (K) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: EMISS ! surface emissivity + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: QSFC ! mixing ratio + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: UST ! U* in similarity theory (m s-1) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TH2_URB2D ! 2 m temperature (not potential) + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: Q2_URB2D ! 2 m mixing ratio + ! end output Noah + TEB + + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_TA + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_PS + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_QA + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_WIND + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_DIR + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_DIR_SW + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_SCA_SW + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_LW + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_RAIN + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_SNOW + REAL, DIMENSION(ims:ime, jms:jme), INTENT(OUT) :: TEB_INPUT_DIR_CO2 + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! LOCAL + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + INTEGER :: IYEAR, IMONTH, IDAY ! date at current time step + REAL :: ZTIME_START ! time at beginning of time step, since midnight (UTC, s) + INTEGER :: I, J ,L + ! Index in the vertical dimension for quantities with half levels. + ! 1 is at the centre of the grid cell + INTEGER, PARAMETER :: KHL = 1 + ! Index in the vertical dimension for quantities with full levels. + ! 1 is at the surface + INTEGER, PARAMETER :: KFL = 1 + REAL :: FRAC_NOAH ! Noah fraction to use during averaging + INTEGER :: URBTYPE ! urban type + REAL :: QA ! specific humidity at the lowest level + REAL :: Q2_RURAL ! 2m specific humidity for vegetation (Noah) + REAL :: T2_RURAL ! 2m temperature for vegetation (Noah) + REAL :: RHO_TSK + REAL :: EMISS_TEMP ! emissivity (Noah + TEB) + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! LOCAL for TEB 1D + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! IN + REAL, DIMENSION(1) :: LAT ! latitude + REAL, DIMENSION(1) :: LON ! longitude + REAL, DIMENSION(1) :: ZGARDEN ! fraction of GARDEN areas + REAL, DIMENSION(1) :: LW ! atmospheric infrared radiation + REAL, DIMENSION(1) :: TA ! air temperature forcing (K) + REAL, DIMENSION(1) :: WIND ! wind speed + REAL, DIMENSION(1) :: DIR ! wind direction + REAL, DIMENSION(1) :: EVAP_TOWN ! evaporation flux (kg/m2/s) + REAL, DIMENSION(1) :: RAIN ! rain rate (kg/m2/s) + REAL, DIMENSION(1) :: SNOW ! snow rate (kg/s/m2 of snow) + REAL, DIMENSION(1) :: PA ! pressure at the first atmospheric level + REAL, DIMENSION(1) :: PS ! pressure at the surface + REAL, DIMENSION(1) :: RHOA ! air density at the lowest level + REAL, DIMENSION(1,1) :: DIR_SW ! incoming direct solar radiation on an horizontal surface + REAL, DIMENSION(1,1) :: SCA_SW ! scattered incoming solar rad. + REAL, DIMENSION(1) :: ZREF ! reference height of the first atmospheric level (temperature) + REAL ,DIMENSION(1) :: QA_KGKG ! air humidity at forcing level (kg/kg) + REAL ,DIMENSION(1) :: BLD ! Horizontal building area density + ! TODO: make into actual input (ims:ime, jms:jme) from file + ! This is not critical for now as we are not considering CO2 in the evaluation. + REAL, DIMENSION(1) :: CO2 = 0. ! CO2 concentration in the air (kg/m3) + + ! OUT + REAL, DIMENSION(1) :: EMIS_TOWN ! town equivalent emissivity + REAL, DIMENSION(1) :: H_TOWN ! sensible heat flux over town + REAL, DIMENSION(1) :: LE_TOWN ! latent heat flux over town + REAL, DIMENSION(1) :: U_CANYON ! canyon hor. wind + REAL, DIMENSION(1) :: ALB_TOWN ! Town eqivalent albedo + REAL, DIMENSION(1) :: DIR_CANYON ! Canyon wind direction + REAL, DIMENSION(1) :: Q_TOWN ! Town eqivalent specific humidity + REAL, DIMENSION(1) :: USTAR_TOWN ! friction velocity over town + REAL, DIMENSION(1) :: GFLUX_TOWN ! flux through the ground + REAL, DIMENSION(1) :: RN_TOWN ! net radiation over town + REAL, DIMENSION(1) :: TS_TOWN ! town surface temperature + + ! Date at simulation start. + CALL cal_mon_day(JULDAY, JULYR, IMONTH, IDAY) + + ! TEB needs year, month, day, GMT time for beginning of current timestep. + IYEAR = JULYR + ZTIME_START = GMT*3600 + XTIME*60 ! seconds since midnight from simulation start + CALL ADD_FORECAST_TO_DATE_SURF(IYEAR, IMONTH, IDAY, ZTIME_START) + ! ZTIME_START is now seconds since midnight at current timestep. + ! IYEAR, IMONTH, IDAY have been adjusted as well to current timestep. + + ! 1. Grid loop + ! ============ + ! Loop over all grid cells and identify whether they are urban. + ! If urban call TEB, otherwise pass through values as given by the LSM or + ! zero in the case of variables that apply to TEB only. + DO J = jts, jte + DO I = its, ite + + ! 2. Run TEB + ! ========== + ! If the cell is of urban type prepare parameters and run TEB. + ! + ! FIXME: does this make sense? I think there should be an and statement to avoid running TEB if + ! it is not defined in the namelist as TEB. If we do that the pass through conditions would be used. + ! We need to check if those pass through conditions would represent something that is representative of + ! urban areas of if they are classified as rural... + IF ( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LOW_DENSITY_RESIDENTIAL .or. & + IVGTYP(I,J) == HIGH_DENSITY_RESIDENTIAL .or. IVGTYP(I,J) == HIGH_INTENSITY_INDUSTRIAL ) THEN + + URBTYPE = UTYPE_URB2D(I,J) + + ! Initialise TEB inputs + ! ========================================= + ! + LAT(1) = XLAT(I,J) + LON(1) = XLONG(I,J) + LW(1) = GLW(I,J) + TA(1) = T3D(I,KHL,J) + + if (TEB_USE_NOAH) then + ! In TEB, `ROAD` is used to account for other impervious surfaces other than buildings. + ! This is defined as: `ROAD = 1 - BLD - GARDEN` + ! When using the Noah LSM to compute the fraction of gardens, we want to run TEB without any gardens + ! as gardens will be accounted by the the Noah LSM. + ! Since `FRAC_URB2D` represents the total impervious surface in a grid box, i.e. buildings + roads in TEB + ! We can obrain the fraction of buildings over the total impervious surface (FRAC_URB2D) as: + BLD(1) = TEB_FRAC_BLD_TBL(URBTYPE) * FRAC_URB2D(I,J) + ZGARDEN(1) = 0. + FRAC_NOAH = 1 - FRAC_URB2D(I,J) + else + FRAC_NOAH = 0 + if (TEB_USE_ZGARDEN_FROM_TBL) then + BLD(1) = TEB_FRAC_BLD_TBL(URBTYPE) + ZGARDEN(1) = TEB_ZGARDEN_TBL(URBTYPE) + else + ! See why this in comment above when TEB_USE_NOAH is ON. + BLD(1) = TEB_FRAC_BLD_TBL(URBTYPE) * FRAC_URB2D(I,J) + ZGARDEN(1) = 1 - FRAC_URB2D(I,J) + endif + endif + + ! + ! Reference height of the first atmospheric level. + ! Mass quantities and horizontal windspeed are at the same height. + ZREF(1) = (0.5 * DZ8W(I,KHL,J)) + IF (ZREF(1) < (TEB_Z0_TOWN_TBL(URBTYPE) * 2)) THEN + CALL WRF_ERROR_FATAL('The first vertical grid spacing must be set to at least ' // & + 'twice the value of roughness lenght for momentum (TEB_Z0_TOWN) ' // & + 'in the urban parameter table (URBPARM.TBL).') + END IF + ! + ! Atmospheric pressure at the surface. + ! NOTE: WRF surface pressure is same as pressure at the roof level in TEB. + ! the pressure is considered to be the same throughout the canyon. + PS(1) = P8W3D(I,KFL,J) + ! + ! Atmospheric pressure at forcing level. + PA(1) = 0.5 * (P8W3D(I,KFL,J) + P8W3D(I,KFL+1,J)) + ! + ! Wind Speed + WIND(1) = CALC_WIND_SPEED_FROM_U_V(U_PHY(I,KHL,J), V_PHY(I,KHL,J)) + ! + ! Wind direction + DIR(1) = CALC_WIND_DIR_FROM_U_V(U_PHY(I,KHL,J), V_PHY(I,KHL,J)) + ! + ! WRF uses humidity in terms of mixing ratio (kg/kg) but TEB require the + ! humidity in terms of specific humidity (kg/kg). + ! Convert from mixing ratio to specific humidity. + QA = GetSpecificHumFromHumRatio(QV3D(I,KHL,J)) + QA_KGKG(1) = GetSpecificHumFromHumRatio(QV3D(I,KHL,J)) + ! + RHOA(1) = RHO(I,KHL,J) + ! + ! RAINBL contains liquid and solid precipitation. + RAIN(1) = 0. + SNOW(1) = 0. + IF (RAINBL(I,J) > 0.) THEN + ! Snow is defined when fraction of frozen precipitation SR(I,J) > 0.5. + IF (SR(I,J) > 0.5 .AND. T3D(I,1,J) <= 273.15) THEN + SNOW(1) = RAINBL(I,J)/XTSTEP_SURF + ELSE + RAIN(1) = RAINBL(I,J)/XTSTEP_SURF + END IF + END IF + ! + ! Radiation + SCA_SW(1,1) = SWDOWN(I,J) * SW_DIFFUSE_FRAC(I,J) + DIR_SW(1,1) = SWDOWN(I,J) - SCA_SW(1,1) + + + !=======================USED IN TECHNICAL EVALUATION=========================================== + if (TEB_TEST_INTEGRATION == 1) then + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + TEB_INPUT_TA(I,J) = T3D(I,KHL,J) + TEB_INPUT_PS(I,J) = PS(1) + TEB_INPUT_QA(I,J) = QA_KGKG(1) + TEB_INPUT_WIND(I,J) = WIND(1) + TEB_INPUT_DIR(I,J) = DIR(1) + TEB_INPUT_DIR_SW(I,J) = DIR_SW(1,1) + TEB_INPUT_SCA_SW(I,J) = SCA_SW(1,1) + TEB_INPUT_LW(I,J) = GLW(I,J) + TEB_INPUT_RAIN(I,J) = RAIN(1) + TEB_INPUT_SNOW(I,J) = SNOW(1) + TEB_INPUT_DIR_CO2(I,J) = CO2(1) + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + ! + ZREF(1) = 33. + RHOA(1) = PS(1) / (T3D(I,KHL,J) * XRD * (1. + ((XRV/XRD) - 1.) * QA) + ZREF(1) * XG) + ! Humidity in kg/m3 + QA = QA * RHOA(1) ! Taken from: src\driver\ol_time_interp_atm.F90 + ! Forcing level pressure from hydrostatism + PA(1) = PS(1) - RHOA(1) * XG * ZREF(1) ! Taken from: src\driver\ol_time_interp_atm.F90 + ! specific humidity (conversion from kg/m3 to kg/kg) + QA_KGKG(1) = QA / RHOA(1) + endif + !=========================USED IN TECHNICAL EVALUATION======================================== + + CALL TEB_DRIVER ( & + ! inputs + TEB_NUM_ROOF_LAYERS, TEB_NUM_WALL_LAYERS, TEB_NUM_ROAD_LAYERS, & + TEB_NUM_SNOW_LAYERS, TEB_NUM_FLOOR_LAYERS, & + LON, LAT, & + IYEAR, IMONTH, IDAY, ZTIME_START, XTSTEP_SURF, & + PS, PA, & + TA, QA_KGKG, RHOA, CO2, & + LW, & + RAIN, SNOW, & + ZREF, & + DIR_SW, SCA_SW, WIND, DIR, & + ! inputs / outputs + TEB_T_CANYON(I,J), TEB_Q_CANYON(I,J), & + TEB_TI_BLD(I,J), & + TEB_T_ROOF(I,:,J), TEB_T_ROAD(I,:,J), & + TEB_T_WALL_A(I,:,J), TEB_T_WALL_B(I,:,J), & + TEB_WS_ROOF(I,J), TEB_WS_ROAD(I,J), & + TEB_WSNOW_ROOF(I,:,J), TEB_TSNOW_ROOF(I,:,J), & + TEB_RSNOW_ROOF(I,:,J), TEB_ASNOW_ROOF(I,J), & + TEB_TSSNOW_ROOF(I,J), TEB_ESNOW_ROOF(I,J), & + TEB_WSNOW_ROAD(I,:,J), TEB_TSNOW_ROAD(I,:,J), & + TEB_RSNOW_ROAD(I,:,J), TEB_ASNOW_ROAD(I,J), & + TEB_TSSNOW_ROAD(I,J), TEB_ESNOW_ROAD(I,J), & + TEB_AUX_MAX(I,J), TEB_T_FLOOR(I,:,J), & + TEB_T_MASS(I,:,J), & + TEB_T_WIN1(I,J), TEB_T_WIN2(I,J), TEB_QI_BLD(I,J), & + TEB_THER_PRODC_DAY(I,J), & + ! parameters + TEB_LGARDEN_TBL(URBTYPE), TEB_GREENROOF_TBL(URBTYPE), TEB_SOLAR_PANEL_TBL(URBTYPE), & + TEB_Z0H_TBL(TEB_OPT_ZOH_TBL(URBTYPE)), TEB_HROAD_DIR_TBL(URBTYPE), TEB_WALL_OPT_TBL(URBTYPE), & + TEB_SNOW_ROAD_TBL(TEB_SNOW_SCH_ROAD_TBL(URBTYPE)), & + TEB_SNOW_ROOF_TBL(TEB_SNOW_SCH_ROOF_TBL(URBTYPE)), & + TEB_H_TRAFFIC_TBL(URBTYPE), TEB_LE_TRAFFIC_TBL(URBTYPE), & + TEB_H_INDUSTRY_TBL(URBTYPE), TEB_LE_INDUSTRY_TBL(URBTYPE), & + TEB_Z0_TOWN_TBL(URBTYPE), & + BLD, ZGARDEN, & + TEB_ZROAD_DIR_TBL(URBTYPE), TEB_FRAC_GR_TBL(URBTYPE), & + TEB_BLD_HEIGHT_TBL(URBTYPE), TEB_WALL_O_TBL(URBTYPE), & + TEB_ALB_ROOF_TBL(URBTYPE), TEB_EMIS_ROOF_TBL(URBTYPE), & + TEB_HC_ROOF_TBL(:, URBTYPE), TEB_TC_ROOF_TBL(:, URBTYPE), TEB_D_ROOF_TBL(:, URBTYPE), & + TEB_ALB_ROAD_TBL(URBTYPE), TEB_EMIS_ROAD_TBL(URBTYPE), & + TEB_HC_ROAD_TBL(:, URBTYPE), TEB_TC_ROAD_TBL(:, URBTYPE), TEB_D_ROAD_TBL(:, URBTYPE), & + TEB_ALB_WALL_TBL(URBTYPE), TEB_EMIS_WALL_TBL(URBTYPE), & + TEB_HC_WALL_TBL(:, URBTYPE), TEB_TC_WALL_TBL(:, URBTYPE), TEB_D_WALL_TBL(:, URBTYPE), & + TEB_COOL_COIL_TBL(URBTYPE), TEB_F_WATER_COND_TBL(URBTYPE), TEB_HEAT_COIL_TBL(URBTYPE), & + TEB_HNATVENT_TBL(URBTYPE), TEB_ZNATVENT_TBL(URBTYPE), & + TEB_F_WASTE_CAN_TBL(URBTYPE), & + TEB_QIN_TBL(URBTYPE), TEB_QIN_FRAD_TBL(URBTYPE), & + TEB_QIN_FLAT_TBL(URBTYPE), TEB_GR_TBL(URBTYPE), TEB_EFF_HEAT_TBL(URBTYPE), TEB_INF_TBL(URBTYPE), & + TEB_TCOOL_TARGET_TBL(URBTYPE), TEB_THEAT_TARGET_TBL(URBTYPE), TEB_HR_TARGET_TBL(URBTYPE), & + TEB_V_VENT_TBL(URBTYPE), TEB_CAP_SYS_HEAT_TBL(URBTYPE), TEB_CAP_SYS_RAT_TBL(URBTYPE), TEB_T_ADP_TBL(URBTYPE), & + TEB_M_SYS_RAT_TBL(URBTYPE), TEB_COP_RAT_TBL(URBTYPE), & + TEB_HC_FLOOR_TBL(:, URBTYPE), & + TEB_TC_FLOOR_TBL(:, URBTYPE), & + TEB_D_FLOOR_TBL(:, URBTYPE), & + TEB_SHGC_TBL(URBTYPE), & + TEB_SHGC_SH_TBL(URBTYPE), & + TEB_LSHADE_TBL(URBTYPE), TEB_ZSHADE_TBL(URBTYPE), & + TEB_CBEM_TBL(URBTYPE), & + TEB_CH_BEM_TBL(URBTYPE), TEB_ROUGH_ROOF_TBL(URBTYPE), TEB_ROUGH_WALL_TBL(URBTYPE), & + TEB_PAR_RD_IRRIG_TBL(URBTYPE), TEB_RD_START_MONTH_TBL(URBTYPE), TEB_RD_END_MONTH_TBL(URBTYPE), & + TEB_RD_START_HOUR_TBL(URBTYPE), TEB_RD_END_HOUR_TBL(URBTYPE), TEB_RD_24H_IRRIG_TBL(URBTYPE), & + TEB_EMIS_PANEL_TBL(URBTYPE), TEB_ALB_PANEL_TBL(URBTYPE), TEB_EFF_PANEL_TBL(URBTYPE), TEB_FRAC_PANEL_TBL(URBTYPE), & + TEB_RESIDENTIAL_TBL(URBTYPE), & + TEB_DT_RES_TBL(URBTYPE), TEB_DT_OFF_TBL(URBTYPE), & + TEB_FLOOR_HEIGHT_TBL(URBTYPE), TEB_U_WIN_TBL(URBTYPE), & + ! outputs + TEB_HVAC_COOL(I,J), TEB_HVAC_HEAT(I,J), & + TEB_THER_PROD_PANEL(I,J), TEB_PHOT_PROD_PANEL(I,J), & + U_CANYON, & + RN_TOWN, H_TOWN, LE_TOWN, GFLUX_TOWN, EVAP_TOWN, & + USTAR_TOWN, & + TS_TOWN, EMIS_TOWN, & + ALB_TOWN, DIR_CANYON, Q_TOWN & + ) + + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! Update of INOUT variables + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + + ! TEB does not provide 10 m wind diagnostics, instead windspeed and direction + ! are provided at half the height of the canyon. + ! TEB only accounts for a reduction in windspeed. Here we recalculate the U and V + ! components with the updated windspeed as returned by TEB. The wind direction remains + ! unchanged, i.e. DIR_CANYON == DIR. + ! Note: U10/V10 is set to U10_URB2D/V10_URB2D in module_surface_driver. + U10_URB2D(I,J) = CALC_U_FROM_WIND_SPEED_AND_DIR(U_CANYON(1), DIR_CANYON(1)) + V10_URB2D(I,J) = CALC_V_FROM_WIND_SPEED_AND_DIR(U_CANYON(1), DIR_CANYON(1)) + + ! Net all-wave radiation (RN_TOWN) does not exist as general WRF variable, + ! only as RN_URB2D, so use that for diagnostic purposes. + RN_URB2D(I,J) = RN_TOWN(1) + + ! Calculate 2 meter temperature and specific humidity for vegetation + IF (CQS2(I,J) < 1.E-5) THEN + Q2_RURAL = QSFC_RURAL(I,J) + ELSE + Q2_RURAL = QSFC_RURAL(I,J) - QFX_RURAL(I,J) / (RHO(I,KHL,J) * CQS2(I,J)) + ENDIF + IF (CHS2(I,J) < 1.E-5) THEN + T2_RURAL = TEB_TSK_RURAL(I,J) + ELSE + T2_RURAL = TEB_TSK_RURAL(I,J) - HFX_RURAL(I,J)/(RHO(I,KHL,J) * XCPD * CHS2(I,J)) + ENDIF + + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! OUT: Town + Vegetation (TEB + Noah) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + ! TEB does not provide 2 m air temperature and humidity diagnostics, + ! instead air temperature and humidity are provided at half the height of the canyon. + ! TH2_URB2D/Q2_URB2D are assigned to Q2/T2/TH2 in the surface driver. + ! Writing directly to Q2 etc. is not possible here as they get overridden. + TH2_URB2D(I,J) = (1.-FRAC_NOAH) * TEB_T_CANYON(I,J) + & + FRAC_NOAH * T2_RURAL + Q2_URB2D(I,J) = (1.-FRAC_NOAH) * GetHumRatioFromSpecificHum(TEB_Q_CANYON(I,J)) + & + FRAC_NOAH * GetHumRatioFromSpecificHum(Q2_RURAL) + EMISS_TEMP = (1.-FRAC_NOAH) * EMIS_TOWN(1) + & + FRAC_NOAH * EMISSI_RURAL(I,J) + EMISS(I,J) = MAX(EMISS_TEMP, 1.E-10) + ALBEDO(I,J) = (1.-FRAC_NOAH) * ALB_TOWN(1) + & + FRAC_NOAH * ALB_RURAL(I,J) + HFX(I,J) = (1.-FRAC_NOAH) * H_TOWN(1) + & + FRAC_NOAH * HFX_RURAL(I,J) + QFX(I,J) = (1.-FRAC_NOAH) * EVAP_TOWN(1) + & + FRAC_NOAH * QFX_RURAL(I,J) + LH(I,J) = (1.-FRAC_NOAH) * LE_TOWN(1) + & + FRAC_NOAH * LH_RURAL(I,J) + GRDFLX(I,J) = (1.-FRAC_NOAH) * (-1)*GFLUX_TOWN(1) + & + FRAC_NOAH * GRDFLX_RURAL(I,J) + TSK(I,J) = (((1.-FRAC_NOAH) * EMIS_TOWN(1) * TS_TOWN(1)**4 + & + FRAC_NOAH * EMISSI_RURAL(I,J) * TEB_TSK_RURAL(I,J)**4) & + / EMISS(I,J))**0.25 + QSFC(I,J) = (1.-FRAC_NOAH) * GetHumRatioFromSpecificHum(Q_TOWN(1)) + & + FRAC_NOAH * QSFC_RURAL(I,J) + UST(I,J) = (1.-FRAC_NOAH) * USTAR_TOWN(1) + & + FRAC_NOAH * UST(I,J) + + ! 3. Pass-through condition + ! ========================= + ! For non-urban cells, return the same conditions as given by the + ! land surface model. + ELSE + ! Calculate 2 meter temperature and specific humidity for vegetation + ! TODO: In the case of water cells, rural diagnostics should be undefined. + IF (CQS2(I,J) < 1.E-5) THEN + Q2_RURAL = QSFC_RURAL(I,J) + ELSE + Q2_RURAL = QSFC_RURAL(I,J) - QFX_RURAL(I,J) / (RHO(I,KHL,J) * CQS2(I,J)) + ENDIF + IF (CHS2(I,J) < 1.E-5) THEN + T2_RURAL = TEB_TSK_RURAL(I,J) + ELSE + T2_RURAL = TEB_TSK_RURAL(I,J) - HFX_RURAL(I,J)/(RHO(I,KHL,J) * XCPD * CHS2(I,J)) + ENDIF + + ! TH2_URB2D/Q2_URB2D are assigned to Q2/T2/TH2 in the surface driver. + ! Writing directly to Q2 etc. is not possible here as they get overridden. + Q2_URB2D(I,J) = Q2_RURAL + TH2_URB2D(I,J) = T2_RURAL + + ! No change + ! Note: UST does not have a _RURAL variant, so no copy needed. + ! Note: TSK is not updated for non-urban cells as done in BEP/BEM. + EMISS(I,J) = EMISSI_RURAL(I,J) + ALBEDO(I,J) = ALB_RURAL(I,J) + HFX(I,J) = HFX_RURAL(I,J) + QFX(I,J) = QFX_RURAL(I,J) + LH(I,J) = LH_RURAL(I,J) + GRDFLX(I,J) = GRDFLX_RURAL(I,J) + QSFC(I,J) = QSFC_RURAL(I,J) + ! Only valid for urban-type grid points + U10_URB2D(I,J) = XUNDEF + V10_URB2D(I,J) = XUNDEF + RN_URB2D(I,J) = XUNDEF + TEB_HVAC_COOL(I,J) = XUNDEF + TEB_HVAC_HEAT(I,J) = XUNDEF + TEB_THER_PROD_PANEL(I,J) = XUNDEF + TEB_PHOT_PROD_PANEL(I,J) = XUNDEF + ! + ENDIF + END DO + END DO + END SUBROUTINE TEB + + +FUNCTION CALC_WIND_SPEED_FROM_U_V(U, V) RESULT(WIND_SPEED) + real, intent(in) :: U, V + real :: WIND_SPEED + + WIND_SPEED = SQRT(U**2. + V**2.) +END FUNCTION CALC_WIND_SPEED_FROM_U_V + + +FUNCTION CALC_WIND_DIR_FROM_U_V(U, V) RESULT(WIND_DIR) + USE MODD_CSTS, ONLY: XPI + real, intent(in) :: U, V + real :: WIND_DIR + + WIND_DIR = 180. + (ATAN2(U, V) * 180./XPI) +END FUNCTION CALC_WIND_DIR_FROM_U_V + + +FUNCTION CALC_U_FROM_WIND_SPEED_AND_DIR(WIND_SPEED, WIND_DIR) RESULT(U) + USE MODD_CSTS, ONLY: XPI + real, intent(in) :: WIND_SPEED, WIND_DIR + real :: U + + U = -WIND_SPEED * SIN(WIND_DIR * XPI/180.) +END FUNCTION CALC_U_FROM_WIND_SPEED_AND_DIR + + +FUNCTION CALC_V_FROM_WIND_SPEED_AND_DIR(WIND_SPEED, WIND_DIR) RESULT(V) + USE MODD_CSTS, ONLY: XPI + real, intent(in) :: WIND_SPEED, WIND_DIR + real :: V + + V = -WIND_SPEED * COS(WIND_DIR * XPI/180.) +END FUNCTION CALC_V_FROM_WIND_SPEED_AND_DIR + + +END MODULE MODULE_SF_TEB \ No newline at end of file diff --git a/phys/module_sf_urban.F b/phys/module_sf_urban.F index e0bffadef4..ad98191f0d 100644 --- a/phys/module_sf_urban.F +++ b/phys/module_sf_urban.F @@ -6,6 +6,10 @@ MODULE module_sf_urban ! Last Update: 2006/08/24 by Fei Chen and Mukul Tewari (NCAR/RAL) !=============================================================================== + REAL, PARAMETER :: UNINIT_REAL = -9999. + INTEGER, PARAMETER :: UNINIT_INT = -9999 + CHARACTER, PARAMETER :: UNINIT_CHAR = '?' + CHARACTER(LEN=4) :: LU_DATA_TYPE INTEGER :: ICATE @@ -90,10 +94,219 @@ MODULE module_sf_urban ! INTEGER :: num_wall_layers ! INTEGER :: num_road_layers + +!for TEB + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ALB_ROOF_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ALB_ROAD_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ALB_WALL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_EMIS_ROOF_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_EMIS_ROAD_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_EMIS_WALL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_BLD_HEIGHT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_FRAC_BLD_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_WALL_O_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_H_TRAFFIC_TBL, TEB_LE_TRAFFIC_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_H_INDUSTRY_TBL, TEB_LE_INDUSTRY_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_Z0_TOWN_TBL + CHARACTER(LEN=6), DIMENSION(3) :: TEB_Z0H_TBL ! not in urbparm.tbl + CHARACTER(LEN=4), DIMENSION(3) :: TEB_SNOW_ROOF_TBL ! not in urbparm.tbl + CHARACTER(LEN=4), DIMENSION(3) :: TEB_SNOW_ROAD_TBL ! not in urbparm.tbl + INTEGER, ALLOCATABLE, DIMENSION(:) :: TEB_OPT_ZOH_TBL + INTEGER, ALLOCATABLE, DIMENSION(:) :: TEB_SNOW_SCH_ROOF_TBL + INTEGER, ALLOCATABLE, DIMENSION(:) :: TEB_SNOW_SCH_ROAD_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_FLOOR_TBL ! function of teb_num_floor_layers + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_MASS_TBL ! function of teb_num_floor_layers + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_TI_BLD_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_CANYON_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_Q_CANYON_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_WIN1_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_WIN2_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_QI_BLD_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_WS_ROOF_TBL ! roof water reservoir(kg/m2) + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_WS_ROAD_TBL ! road water reservoir(kg/m2) + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_WSNOW_ROOF_TBL ! snow layers reservoir + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_TSNOW_ROOF_TBL ! snow layers temperature + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RSNOW_ROOF_TBL ! snow layers density + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ASNOW_ROOF_TBL ! snow albedo + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ESNOW_ROOF_TBL ! snow emissivity + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_TSSNOW_ROOF_TBL ! snow surface temperature + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_WSNOW_ROAD_TBL ! snow layers reservoir + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_TSNOW_ROAD_TBL ! snow layers temperature + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RSNOW_ROAD_TBL ! snow layers density + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ASNOW_ROAD_TBL ! snow albedo + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ESNOW_ROAD_TBL ! snow emissivity + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_TSSNOW_ROAD_TBL ! snow surface temperature + LOGICAL, ALLOCATABLE, DIMENSION(:) :: TEB_SOLAR_PANEL_TBL + LOGICAL :: TEB_SOLAR_PANEL_INITED = .FALSE. + LOGICAL, ALLOCATABLE, DIMENSION(:) :: TEB_LGARDEN_TBL + LOGICAL :: TEB_LGARDEN_INITED = .FALSE. + LOGICAL, ALLOCATABLE, DIMENSION(:) :: TEB_GREENROOF_TBL + LOGICAL :: TEB_GREENROOF_INITED = .FALSE. + CHARACTER(LEN=4), ALLOCATABLE, DIMENSION(:) :: TEB_HROAD_DIR_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ZROAD_DIR_TBL + CHARACTER(LEN=4), ALLOCATABLE, DIMENSION(:) :: TEB_WALL_OPT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ZGARDEN_TBL + LOGICAL :: TEB_USE_ZGARDEN_FROM_TBL + LOGICAL :: TEB_USE_ZGARDEN_FROM_TBL_INITED = .FALSE. + LOGICAL :: TEB_USE_NOAH + LOGICAL :: TEB_USE_NOAH_INITED = .FALSE. + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_FRAC_GR_TBL + CHARACTER(LEN=12),ALLOCATABLE, DIMENSION(:) :: TEB_COOL_COIL_TBL + CHARACTER(LEN=6), ALLOCATABLE, DIMENSION(:) :: TEB_HEAT_COIL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_F_WATER_COND_TBL + CHARACTER(LEN=4), ALLOCATABLE, DIMENSION(:) :: TEB_HNATVENT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ZNATVENT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_F_WASTE_CAN_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_QIN_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_QIN_FRAD_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_QIN_FLAT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_GR_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_EFF_HEAT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_INF_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_TCOOL_TARGET_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_THEAT_TARGET_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_HR_TARGET_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_V_VENT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_CAP_SYS_HEAT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_CAP_SYS_RAT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_ADP_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_M_SYS_RAT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_COP_RAT_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_SHGC_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_SHGC_SH_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_U_WIN_TBL + LOGICAL, ALLOCATABLE, DIMENSION(:) :: TEB_LSHADE_TBL + LOGICAL :: TEB_LSHADE_INITED + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ZSHADE_TBL + CHARACTER(LEN=3), ALLOCATABLE, DIMENSION(:) :: TEB_CBEM_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_FLOOR_HEIGHT_TBL + CHARACTER(LEN=5), ALLOCATABLE, DIMENSION(:) :: TEB_CH_BEM_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ROUGH_ROOF_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ROUGH_WALL_TBL + LOGICAL, ALLOCATABLE, DIMENSION(:) :: TEB_PAR_RD_IRRIG_TBL + LOGICAL :: TEB_PAR_RD_IRRIG_INITED + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RD_START_MONTH_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RD_END_MONTH_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RD_START_HOUR_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RD_END_HOUR_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RD_24H_IRRIG_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_EMIS_PANEL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_ALB_PANEL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_EFF_PANEL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_FRAC_PANEL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_RESIDENTIAL_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_DT_RES_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_DT_OFF_TBL + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_ROAD_TBL ! function of teb_num_floor_layers + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_ROOF_TBL ! function of teb_num_floor_layers + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_WALL_A_TBL ! function of teb_num_floor_layers + REAL, ALLOCATABLE, DIMENSION(:) :: TEB_T_WALL_B_TBL ! function of teb_num_floor_layers + + ! ROOF + INTEGER, ALLOCATABLE, DIMENSION(:) :: teb_numroof_tbl + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_HC_ROOF_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_TC_ROOF_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_D_ROOF_TBL + + ! ROAD + INTEGER, ALLOCATABLE, DIMENSION(:) :: teb_numroad_tbl + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_HC_ROAD_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_TC_ROAD_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_D_ROAD_TBL + + ! WALL + INTEGER, ALLOCATABLE, DIMENSION(:) :: teb_numwall_tbl + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_HC_WALL_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_TC_WALL_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_D_WALL_TBL + + ! FLOOR + INTEGER, ALLOCATABLE, DIMENSION(:) :: teb_numfloor_tbl + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_HC_FLOOR_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_TC_FLOOR_TBL + REAL, ALLOCATABLE, DIMENSION(:,:) :: TEB_D_FLOOR_TBL + +!end TEB + + + CHARACTER (LEN=256) , PRIVATE :: mesg CONTAINS + + SUBROUTINE check_init_integer(x, name) + IMPLICIT NONE + INTEGER, INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: name + if (x == UNINIT_INT) then + CALL wrf_error_fatal('URBPARM.TBL: "'//name//'" uninitialized') + end if + END SUBROUTINE + + SUBROUTINE check_init_integer_array(x, name) + IMPLICIT NONE + INTEGER, DIMENSION(:), INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER :: i + do i = 1, size(x) + if (x(i) == UNINIT_INT) then + CALL wrf_error_fatal('URBPARM.TBL: "'//name//'" uninitialized') + end if + enddo + END SUBROUTINE + + SUBROUTINE check_init_real(x, name) + IMPLICIT NONE + REAL, INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: name + if (x == UNINIT_REAL) then + CALL wrf_error_fatal('URBPARM.TBL: "'//name//'" uninitialized') + end if + END SUBROUTINE + + SUBROUTINE check_init_real_array(x, name) + IMPLICIT NONE + REAL, DIMENSION(:), INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER :: i + do i = 1, size(x) + if (x(i) == UNINIT_REAL) then + CALL wrf_error_fatal('URBPARM.TBL: "'//name//'" uninitialized') + end if + enddo + END SUBROUTINE + + SUBROUTINE check_init_string(x, name) + IMPLICIT NONE + CHARACTER(LEN=*), INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: name + if (x == repeat(UNINIT_CHAR, len(x))) then + CALL wrf_error_fatal('URBPARM.TBL: "'//name//'" uninitialized') + end if + END SUBROUTINE + + SUBROUTINE check_init_string_array(x, name) + IMPLICIT NONE + CHARACTER(LEN=*), DIMENSION(:), INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: name + INTEGER :: i + do i = 1, size(x) + if (x(i) == repeat(UNINIT_CHAR, len(x(i)))) then + CALL wrf_error_fatal('URBPARM.TBL: "'//name//'" uninitialized') + end if + enddo + END SUBROUTINE + + SUBROUTINE check_init_logical(x, name) + IMPLICIT NONE + LOGICAL, INTENT(IN) :: x + CHARACTER(LEN=*), INTENT(IN) :: name + if (.not. x) then + CALL wrf_error_fatal('URBPARM.TBL: "'//name//'" uninitialized') + end if + END SUBROUTINE + !=============================================================================== ! ! Author: @@ -304,6 +517,11 @@ SUBROUTINE urban(LSOLAR, & ! L CMCR,TGR,TGRL,SMR,CMGR_URB,CHGR_URB,jmonth, & ! H DRELR,DRELB,DRELG,FLXHUMR,FLXHUMB,FLXHUMG) + ! TEB + USE MODI_INI_CSTS, ONLY: INI_CSTS + USE MODI_INIT_SURFCONSPHY, ONLY: INIT_SURFCONSPHY + ! end TEB + IMPLICIT NONE REAL, PARAMETER :: CP=0.24 ! heat capacity of dry air [cgs unit] @@ -1934,12 +2152,22 @@ END SUBROUTINE read_param ! !=============================================================================== SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & - sf_urban_physics) + sf_urban_physics, & + teb_num_snow_layers, & + teb_num_road_layers, & + teb_num_roof_layers, & + teb_num_wall_layers, & + teb_num_floor_layers) ! num_roof_layers,num_wall_layers,num_road_layers) IMPLICIT NONE INTEGER, INTENT(IN) :: num_soil_layers + INTEGER, INTENT(IN) :: teb_num_snow_layers + INTEGER, INTENT(IN) :: teb_num_road_layers + INTEGER, INTENT(IN) :: teb_num_roof_layers + INTEGER, INTENT(IN) :: teb_num_wall_layers + INTEGER, INTENT(IN) :: teb_num_floor_layers ! REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: DZR ! REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: DZB @@ -2086,6 +2314,345 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & ! if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROOF_WIDTH in urban_param_init') ! ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status ) ! if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROAD_WIDTH in urban_param_init') + !for TEB + ALLOCATE( TEB_ALB_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ALB_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_ALB_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ALB_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_ALB_WALL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ALB_WALL_TBL in urban_param_init') + ALLOCATE( TEB_EMIS_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_EMIS_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_EMIS_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_EMIS_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_EMIS_WALL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_EMIS_WALL_TBL in urban_param_init') + ALLOCATE( TEB_BLD_HEIGHT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_BLD_HEIGHT_TBL in urban_param_init') + ALLOCATE( TEB_FRAC_BLD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_FRAC_BLD_TBL in urban_param_init') + ALLOCATE( TEB_WALL_O_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_WALL_O_TBL in urban_param_init') + ALLOCATE( TEB_H_TRAFFIC_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_H_TRAFFIC_TBL in urban_param_init') + ALLOCATE( TEB_LE_TRAFFIC_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_LE_TRAFFIC_TBL in urban_param_init') + ALLOCATE( TEB_H_INDUSTRY_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_H_INDUSTRY_TBL in urban_param_init') + ALLOCATE( TEB_LE_INDUSTRY_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_LE_INDUSTRY_TBL in urban_param_init') + ALLOCATE( TEB_Z0_TOWN_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_Z0_TOWN_TBL in urban_param_init') + ALLOCATE( TEB_OPT_ZOH_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_OPT_ZOH_TBL in urban_param_init') + ALLOCATE( TEB_SNOW_SCH_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_SNOW_SCH_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_SNOW_SCH_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_SNOW_SCH_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_TI_BLD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TI_BLD_TBL in urban_param_init') + ALLOCATE( TEB_T_CANYON_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_CANYON_TBL in urban_param_init') + ALLOCATE( TEB_Q_CANYON_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_Q_CANYON_TBL in urban_param_init') + ALLOCATE( TEB_T_WIN1_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_WIN1_TBL in urban_param_init') + ALLOCATE( TEB_T_WIN2_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_WIN2_TBL in urban_param_init') + ALLOCATE( TEB_QI_BLD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_QI_BLD_TBL in urban_param_init') + ALLOCATE( TEB_WS_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_WS_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_WS_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_WS_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_WSNOW_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_WSNOW_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_TSNOW_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TSNOW_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_RSNOW_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RSNOW_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_ASNOW_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ASNOW_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_ESNOW_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ESNOW_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_TSSNOW_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TSSNOW_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_WSNOW_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_WSNOW_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_TSNOW_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TSNOW_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_RSNOW_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RSNOW_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_ASNOW_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ASNOW_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_ESNOW_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ESNOW_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_TSSNOW_ROAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TSSNOW_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_SOLAR_PANEL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_SOLAR_PANEL_TBL in urban_param_init') + ALLOCATE( TEB_LGARDEN_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_LGARDEN_TBL in urban_param_init') + ALLOCATE( TEB_GREENROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_GREENROOF_TBL in urban_param_init') + ALLOCATE( TEB_HROAD_DIR_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HROAD_DIR_TBL in urban_param_init') + ALLOCATE( TEB_ZROAD_DIR_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ZROAD_DIR_TBL in urban_param_init') + ALLOCATE( TEB_WALL_OPT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_WALL_OPT_TBL in urban_param_init') + ALLOCATE( TEB_ZGARDEN_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ZGARDEN_TBL in urban_param_init') + ALLOCATE( TEB_FRAC_GR_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_FRAC_GR_TBL in urban_param_init') + ALLOCATE( TEB_COOL_COIL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_COOL_COIL_TBL in urban_param_init') + ALLOCATE( TEB_HEAT_COIL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HEAT_COIL_TBL in urban_param_init') + ALLOCATE( TEB_F_WATER_COND_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_F_WATER_COND_TBL in urban_param_init') + ALLOCATE( TEB_HNATVENT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HNATVENT_TBL in urban_param_init') + ALLOCATE( TEB_ZNATVENT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ZNATVENT_TBL in urban_param_init') + ALLOCATE( TEB_F_WASTE_CAN_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_F_WASTE_CAN_TBL in urban_param_init') + ALLOCATE( TEB_QIN_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_QIN_TBL in urban_param_init') + ALLOCATE( TEB_QIN_FRAD_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_QIN_FRAD_TBL in urban_param_init') + ALLOCATE( TEB_QIN_FLAT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_QIN_FLAT_TBL in urban_param_init') + ALLOCATE( TEB_GR_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_GR_TBL in urban_param_init') + ALLOCATE( TEB_EFF_HEAT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_EFF_HEAT_TBL in urban_param_init') + ALLOCATE( TEB_INF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_INF_TBL in urban_param_init') + ALLOCATE( TEB_TCOOL_TARGET_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TCOOL_TARGET_TBL in urban_param_init') + ALLOCATE( TEB_THEAT_TARGET_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_THEAT_TARGET_TBL in urban_param_init') + ALLOCATE( TEB_HR_TARGET_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HR_TARGET_TBL in urban_param_init') + ALLOCATE( TEB_V_VENT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_V_VENT_TBL in urban_param_init') + ALLOCATE( TEB_CAP_SYS_HEAT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_CAP_SYS_HEAT_TBL in urban_param_init') + ALLOCATE( TEB_CAP_SYS_RAT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_CAP_SYS_RAT_TBL in urban_param_init') + ALLOCATE( TEB_T_ADP_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_ADP_TBL in urban_param_init') + ALLOCATE( TEB_M_SYS_RAT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_M_SYS_RAT_TBL in urban_param_init') + ALLOCATE( TEB_COP_RAT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_COP_RAT_TBL in urban_param_init') + ALLOCATE( TEB_SHGC_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_SHGC_TBL in urban_param_init') + ALLOCATE( TEB_SHGC_SH_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_SHGC_SH_TBL in urban_param_init') + ALLOCATE( TEB_U_WIN_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_U_WIN_TBL in urban_param_init') + ALLOCATE( TEB_LSHADE_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_LSHADE_TBL in urban_param_init') + ALLOCATE( TEB_ZSHADE_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ZSHADE_TBL in urban_param_init') + ALLOCATE( TEB_CBEM_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_CBEM_TBL in urban_param_init') + ALLOCATE( TEB_FLOOR_HEIGHT_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_FLOOR_HEIGHT_TBL in urban_param_init') + ALLOCATE( TEB_CH_BEM_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_CH_BEM_TBL in urban_param_init') + ALLOCATE( TEB_ROUGH_ROOF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ROUGH_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_ROUGH_WALL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ROUGH_WALL_TBL in urban_param_init') + ALLOCATE( TEB_PAR_RD_IRRIG_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_PAR_RD_IRRIG_TBL in urban_param_init') + ALLOCATE( TEB_RD_START_MONTH_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RD_START_MONTH_TBL in urban_param_init') + ALLOCATE( TEB_RD_END_MONTH_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RD_END_MONTH_TBL in urban_param_init') + ALLOCATE( TEB_RD_START_HOUR_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RD_START_HOUR_TBL in urban_param_init') + ALLOCATE( TEB_RD_END_HOUR_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RD_END_HOUR_TBL in urban_param_init') + ALLOCATE( TEB_RD_24H_IRRIG_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RD_24H_IRRIG_TBL in urban_param_init') + ALLOCATE( TEB_EMIS_PANEL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_EMIS_PANEL_TBL in urban_param_init') + ALLOCATE( TEB_ALB_PANEL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_ALB_PANEL_TBL in urban_param_init') + ALLOCATE( TEB_EFF_PANEL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_EFF_PANEL_TBL in urban_param_init') + ALLOCATE( TEB_FRAC_PANEL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_FRAC_PANEL_TBL in urban_param_init') + ALLOCATE( TEB_RESIDENTIAL_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_RESIDENTIAL_TBL in urban_param_init') + ALLOCATE( TEB_DT_RES_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_DT_RES_TBL in urban_param_init') + ALLOCATE( TEB_DT_OFF_TBL(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_DT_OFF_TBL in urban_param_init') + + ALLOCATE( TEB_T_FLOOR_TBL(teb_num_floor_layers), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_FLOOR in urban_param_init') + ALLOCATE( TEB_T_MASS_TBL(teb_num_floor_layers), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_MASS in urban_param_init') + + ALLOCATE( TEB_T_ROAD_TBL(teb_num_road_layers), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_T_ROOF_TBL(teb_num_roof_layers), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_T_WALL_A_TBL(teb_num_wall_layers), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_WALL_A_TBL in urban_param_init') + ALLOCATE( TEB_T_WALL_B_TBL(teb_num_wall_layers), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_T_WALL_B_TBL in urban_param_init') + ! ROOF + ALLOCATE( teb_numroof_tbl(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating teb_numroof_tbl in urban_param_init') + ALLOCATE( TEB_HC_ROOF_TBL(teb_num_roof_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HC_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_TC_ROOF_TBL(teb_num_roof_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TC_ROOF_TBL in urban_param_init') + ALLOCATE( TEB_D_ROOF_TBL(teb_num_roof_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_D_ROOF_TBL in urban_param_init') + ! ROAD + ALLOCATE( teb_numroad_tbl(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating teb_numroad_tbl in urban_param_init') + ALLOCATE( TEB_HC_ROAD_TBL(teb_num_road_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HC_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_TC_ROAD_TBL(teb_num_road_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TC_ROAD_TBL in urban_param_init') + ALLOCATE( TEB_D_ROAD_TBL(teb_num_road_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_D_ROAD_TBL in urban_param_init') + ! WALL + ALLOCATE( teb_numwall_tbl(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating teb_numwall_tbl in urban_param_init') + ALLOCATE( TEB_HC_WALL_TBL(teb_num_wall_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HC_WALL_TBL in urban_param_init') + ALLOCATE( TEB_TC_WALL_TBL(teb_num_wall_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TC_WALL_TBL in urban_param_init') + ALLOCATE( TEB_D_WALL_TBL(teb_num_wall_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_D_WALL_TBL in urban_param_init') + ! FLOOR + ALLOCATE( teb_numfloor_tbl(ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating teb_numfloor_tbl in urban_param_init') + ALLOCATE( TEB_HC_FLOOR_TBL(teb_num_floor_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_HC_FLOOR_TBL in urban_param_init') + ALLOCATE( TEB_TC_FLOOR_TBL(teb_num_floor_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_TC_FLOOR_TBL in urban_param_init') + ALLOCATE( TEB_D_FLOOR_TBL(teb_num_floor_layers , ICATE), stat=allocate_status ) + if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TEB_D_FLOOR_TBL in urban_param_init') + + TEB_EMIS_ROOF_TBL = UNINIT_REAL + TEB_EMIS_ROAD_TBL = UNINIT_REAL + TEB_EMIS_WALL_TBL = UNINIT_REAL + TEB_ALB_ROOF_TBL = UNINIT_REAL + TEB_ALB_ROAD_TBL = UNINIT_REAL + TEB_ALB_WALL_TBL = UNINIT_REAL + TEB_BLD_HEIGHT_TBL = UNINIT_REAL + TEB_FRAC_BLD_TBL = UNINIT_REAL + TEB_WALL_O_TBL = UNINIT_REAL + TEB_H_TRAFFIC_TBL = UNINIT_REAL + TEB_LE_TRAFFIC_TBL = UNINIT_REAL + TEB_H_INDUSTRY_TBL = UNINIT_REAL + TEB_LE_INDUSTRY_TBL = UNINIT_REAL + TEB_Z0_TOWN_TBL = UNINIT_REAL + TEB_OPT_ZOH_TBL = UNINIT_INT + TEB_SNOW_SCH_ROOF_TBL = UNINIT_INT + TEB_SNOW_SCH_ROAD_TBL = UNINIT_INT + TEB_TI_BLD_TBL = UNINIT_REAL + TEB_T_CANYON_TBL = UNINIT_REAL + TEB_Q_CANYON_TBL = UNINIT_REAL + TEB_T_WIN1_TBL = UNINIT_REAL + TEB_T_WIN2_TBL = UNINIT_REAL + TEB_QI_BLD_TBL = UNINIT_REAL + TEB_WS_ROOF_TBL = UNINIT_REAL + TEB_WS_ROAD_TBL = UNINIT_REAL + TEB_WSNOW_ROOF_TBL = UNINIT_REAL + TEB_TSNOW_ROOF_TBL = UNINIT_REAL + TEB_RSNOW_ROOF_TBL = UNINIT_REAL + TEB_ASNOW_ROOF_TBL = UNINIT_REAL + TEB_ESNOW_ROOF_TBL = UNINIT_REAL + TEB_TSSNOW_ROOF_TBL = UNINIT_REAL + TEB_WSNOW_ROAD_TBL = UNINIT_REAL + TEB_TSNOW_ROAD_TBL = UNINIT_REAL + TEB_RSNOW_ROAD_TBL = UNINIT_REAL + TEB_ASNOW_ROAD_TBL = UNINIT_REAL + TEB_ESNOW_ROAD_TBL = UNINIT_REAL + TEB_TSSNOW_ROAD_TBL = UNINIT_REAL + TEB_HROAD_DIR_TBL = repeat(UNINIT_CHAR, 4) + TEB_ZROAD_DIR_TBL = UNINIT_REAL + TEB_WALL_OPT_TBL = repeat(UNINIT_CHAR, 4) + TEB_ZGARDEN_TBL = UNINIT_REAL + TEB_FRAC_GR_TBL = UNINIT_REAL + TEB_COOL_COIL_TBL = repeat(UNINIT_CHAR, 12) + TEB_HEAT_COIL_TBL = repeat(UNINIT_CHAR, 6) + TEB_F_WATER_COND_TBL = UNINIT_REAL + TEB_HNATVENT_TBL = repeat(UNINIT_CHAR, 4) + TEB_ZNATVENT_TBL = UNINIT_REAL + TEB_F_WASTE_CAN_TBL = UNINIT_REAL + TEB_QIN_TBL = UNINIT_REAL + TEB_QIN_FRAD_TBL = UNINIT_REAL + TEB_QIN_FLAT_TBL = UNINIT_REAL + TEB_GR_TBL = UNINIT_REAL + TEB_EFF_HEAT_TBL = UNINIT_REAL + TEB_INF_TBL = UNINIT_REAL + TEB_TCOOL_TARGET_TBL = UNINIT_REAL + TEB_THEAT_TARGET_TBL = UNINIT_REAL + TEB_HR_TARGET_TBL = UNINIT_REAL + TEB_V_VENT_TBL = UNINIT_REAL + TEB_CAP_SYS_HEAT_TBL = UNINIT_REAL + TEB_CAP_SYS_RAT_TBL = UNINIT_REAL + TEB_T_ADP_TBL = UNINIT_REAL + TEB_M_SYS_RAT_TBL = UNINIT_REAL + TEB_COP_RAT_TBL = UNINIT_REAL + TEB_SHGC_TBL = UNINIT_REAL + TEB_SHGC_SH_TBL = UNINIT_REAL + TEB_U_WIN_TBL = UNINIT_REAL + TEB_ZSHADE_TBL = UNINIT_REAL + TEB_CBEM_TBL = repeat(UNINIT_CHAR, 3) + TEB_FLOOR_HEIGHT_TBL = UNINIT_REAL + TEB_CH_BEM_TBL = repeat(UNINIT_CHAR, 5) + TEB_ROUGH_ROOF_TBL = UNINIT_REAL + TEB_ROUGH_WALL_TBL = UNINIT_REAL + TEB_RD_START_MONTH_TBL = UNINIT_REAL + TEB_RD_END_MONTH_TBL = UNINIT_REAL + TEB_RD_START_HOUR_TBL = UNINIT_REAL + TEB_RD_END_HOUR_TBL = UNINIT_REAL + TEB_RD_24H_IRRIG_TBL = UNINIT_REAL + TEB_EMIS_PANEL_TBL = UNINIT_REAL + TEB_ALB_PANEL_TBL = UNINIT_REAL + TEB_EFF_PANEL_TBL = UNINIT_REAL + TEB_FRAC_PANEL_TBL = UNINIT_REAL + TEB_RESIDENTIAL_TBL = UNINIT_REAL + TEB_DT_RES_TBL = UNINIT_REAL + TEB_DT_OFF_TBL = UNINIT_REAL + + TEB_T_FLOOR_TBL = UNINIT_REAL + TEB_T_MASS_TBL = UNINIT_REAL + TEB_T_ROAD_TBL = UNINIT_REAL + TEB_T_ROOF_TBL = UNINIT_REAL + TEB_T_WALL_A_TBL = UNINIT_REAL + TEB_T_WALL_B_TBL = UNINIT_REAL + + ! ROOF + TEB_HC_ROOF_TBL = UNINIT_REAL + TEB_TC_ROOF_TBL = UNINIT_REAL + TEB_D_ROOF_TBL = UNINIT_REAL + ! ROAD + TEB_HC_ROAD_TBL = UNINIT_REAL + TEB_TC_ROAD_TBL = UNINIT_REAL + TEB_D_ROAD_TBL = UNINIT_REAL + ! WALL + TEB_HC_WALL_TBL = UNINIT_REAL + TEB_TC_WALL_TBL = UNINIT_REAL + TEB_D_WALL_TBL = UNINIT_REAL + ! FLOOR + TEB_HC_FLOOR_TBL = UNINIT_REAL + TEB_TC_FLOOR_TBL = UNINIT_REAL + TEB_D_FLOOR_TBL = UNINIT_REAL + !end TEB !for BEP ALLOCATE( NUMDIR_TBL(ICATE), stat=allocate_status ) if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating NUMDIR_TBL in urban_param_init') @@ -2137,6 +2704,13 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & numhgt_tbl = 0 height_bin_tbl = -1.E36 hpercent_bin_tbl = -1.E36 + + ! TEB + ! Note: This must be here as urban_param_init may be called more than once. + teb_numroof_tbl = 0 + teb_numroad_tbl = 0 + teb_numwall_tbl = 0 + teb_numfloor_tbl = 0 !end BEP else if (name == "ZR") then @@ -2322,6 +2896,268 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & else if (name == "HSEQUIP_SCALE_FACTOR") then read(string(indx+1:),*) hsesf_tbl(1:icate) !end BEP + +!for TEB + else if (name == "TEB_USE_ZGARDEN_FROM_TBL") then + read(string(indx+1:),*) TEB_USE_ZGARDEN_FROM_TBL + TEB_USE_ZGARDEN_FROM_TBL_INITED = .true. + else if (name == "TEB_USE_NOAH") then + read(string(indx+1:),*) TEB_USE_NOAH + TEB_USE_NOAH_INITED = .true. + else if (name == "TEB_ALB_ROOF") then + read(string(indx+1:),*) TEB_ALB_ROOF_TBL(1:icate) + else if (name == "TEB_ALB_ROAD") then + read(string(indx+1:),*) TEB_ALB_ROAD_TBL(1:icate) + else if (name == "TEB_ALB_WALL") then + read(string(indx+1:),*) TEB_ALB_WALL_TBL(1:icate) + else if (name == "TEB_EMIS_ROOF") then + read(string(indx+1:),*) TEB_EMIS_ROOF_TBL(1:icate) + else if (name == "TEB_EMIS_ROAD") then + read(string(indx+1:),*) TEB_EMIS_ROAD_TBL(1:icate) + else if (name == "TEB_EMIS_WALL") then + read(string(indx+1:),*) TEB_EMIS_WALL_TBL(1:icate) + else if (name == "TEB_BLD_HEIGHT") then + read(string(indx+1:),*) TEB_BLD_HEIGHT_TBL(1:icate) + else if (name == "TEB_FRAC_BLD") then + read(string(indx+1:),*) TEB_FRAC_BLD_TBL(1:icate) + else if (name == "TEB_WALL_O_HOR") then + read(string(indx+1:),*) TEB_WALL_O_TBL(1:icate) + else if (name == "TEB_H_TRAFFIC") then + read(string(indx+1:),*) TEB_H_TRAFFIC_TBL(1:icate) + else if (name == "TEB_LE_TRAFFIC") then + read(string(indx+1:),*) TEB_LE_TRAFFIC_TBL(1:icate) + else if (name == "TEB_H_INDUSTRY") then + read(string(indx+1:),*) TEB_H_INDUSTRY_TBL(1:icate) + else if (name == "TEB_LE_INDUSTRY") then + read(string(indx+1:),*) TEB_LE_INDUSTRY_TBL(1:icate) + else if (name == "TEB_Z0_TOWN") then + read(string(indx+1:),*) TEB_Z0_TOWN_TBL(1:icate) + else if (name == "TEB_SNOW_SCH_ROOF") then + read(string(indx+1:),*) TEB_SNOW_SCH_ROOF_TBL(1:icate) + else if (name == "TEB_SNOW_SCH_ROAD") then + read(string(indx+1:),*) TEB_SNOW_SCH_ROAD_TBL(1:icate) + else if (name == "TEB_OPT_ZOH") then + read(string(indx+1:),*)TEB_OPT_ZOH_TBL(1:icate) + else if (name == "TEB_SOLAR_PANEL") then + read(string(indx+1:),*)TEB_SOLAR_PANEL_TBL(1:icate) + TEB_SOLAR_PANEL_INITED = .TRUE. + else if (name == "TEB_LGARDEN") then + read(string(indx+1:),*)TEB_LGARDEN_TBL(1:icate) + TEB_LGARDEN_INITED = .TRUE. + else if (name == "TEB_GREENROOF") then + read(string(indx+1:),*)TEB_GREENROOF_TBL(1:icate) + TEB_GREENROOF_INITED = .TRUE. + else if (name == "TEB_HROAD_DIR") then + read(string(indx+1:),*)TEB_HROAD_DIR_TBL(1:icate) + else if (name == "TEB_WALL_OPT") then + read(string(indx+1:),*)TEB_WALL_OPT_TBL(1:icate) + else if (name == "TEB_ZGARDEN") then + read(string(indx+1:),*)TEB_ZGARDEN_TBL(1:icate) + else if (name == "TEB_ZROAD_DIR") then + read(string(indx+1:),*)TEB_ZROAD_DIR_TBL(1:icate) + else if (name == "TEB_FRAC_GR") then + read(string(indx+1:),*)TEB_FRAC_GR_TBL(1:icate) + else if (name == "TEB_COOL_COIL") then + read(string(indx+1:),*)TEB_COOL_COIL_TBL(1:icate) + else if (name == "TEB_HEAT_COIL") then + read(string(indx+1:),*)TEB_HEAT_COIL_TBL(1:icate) + else if (name == "TEB_F_WATER_COND") then + read(string(indx+1:),*)TEB_F_WATER_COND_TBL(1:icate) + else if (name == "TEB_HNATVENT") then + read(string(indx+1:),*)TEB_HNATVENT_TBL(1:icate) + else if (name == "TEB_ZNATVENT") then + read(string(indx+1:),*)TEB_ZNATVENT_TBL(1:icate) + else if (name == "TEB_F_WASTE_CAN") then + read(string(indx+1:),*)TEB_F_WASTE_CAN_TBL(1:icate) + else if (name == "TEB_QIN") then + read(string(indx+1:),*)TEB_QIN_TBL(1:icate) + else if (name == "TEB_QIN_FRAD") then + read(string(indx+1:),*)TEB_QIN_FRAD_TBL(1:icate) + else if (name == "TEB_QIN_FLAT") then + read(string(indx+1:),*)TEB_QIN_FLAT_TBL(1:icate) + else if (name == "TEB_GR") then + read(string(indx+1:),*)TEB_GR_TBL(1:icate) + else if (name == "TEB_EFF_HEAT") then + read(string(indx+1:),*)TEB_EFF_HEAT_TBL(1:icate) + else if (name == "TEB_INF") then + read(string(indx+1:),*)TEB_INF_TBL(1:icate) + else if (name == "TEB_TCOOL_TARGET") then + read(string(indx+1:),*)TEB_TCOOL_TARGET_TBL(1:icate) + else if (name == "TEB_THEAT_TARGET") then + read(string(indx+1:),*)TEB_THEAT_TARGET_TBL(1:icate) + else if (name == "TEB_HR_TARGET") then + read(string(indx+1:),*)TEB_HR_TARGET_TBL(1:icate) + else if (name == "TEB_V_VENT") then + read(string(indx+1:),*)TEB_V_VENT_TBL(1:icate) + else if (name == "TEB_CAP_SYS_HEAT") then + read(string(indx+1:),*)TEB_CAP_SYS_HEAT_TBL(1:icate) + else if (name == "TEB_CAP_SYS_RAT") then + read(string(indx+1:),*)TEB_CAP_SYS_RAT_TBL(1:icate) + else if (name == "TEB_T_ADP") then + read(string(indx+1:),*)TEB_T_ADP_TBL(1:icate) + else if (name == "TEB_M_SYS_RAT") then + read(string(indx+1:),*)TEB_M_SYS_RAT_TBL(1:icate) + else if (name == "TEB_COP_RAT") then + read(string(indx+1:),*)TEB_COP_RAT_TBL(1:icate) + else if (name == "TEB_T_FLOOR") then + read(string(indx+1:),*)TEB_T_FLOOR_TBL(1:teb_num_floor_layers) + else if (name == "TEB_T_MASS") then + read(string(indx+1:),*)TEB_T_MASS_TBL(1:teb_num_floor_layers) + else if (name == "TEB_TI_BLD") then + read(string(indx+1:),*)TEB_TI_BLD_TBL(1:icate) + else if (name == "TEB_T_CANYON") then + read(string(indx+1:),*)TEB_T_CANYON_TBL(1:icate) + else if (name == "TEB_Q_CANYON") then + read(string(indx+1:),*)TEB_Q_CANYON_TBL(1:icate) + else if (name == "TEB_T_WIN1") then + read(string(indx+1:),*)TEB_T_WIN1_TBL(1:icate) + else if (name == "TEB_T_WIN2") then + read(string(indx+1:),*)TEB_T_WIN2_TBL(1:icate) + else if (name == "TEB_QI_BLD") then + read(string(indx+1:),*)TEB_QI_BLD_TBL(1:icate) + else if (name == "TEB_WS_ROOF") then + read(string(indx+1:),*)TEB_WS_ROOF_TBL(1:icate) + else if (name == "TEB_WS_ROAD") then + read(string(indx+1:),*)TEB_WS_ROAD_TBL(1:icate) + else if (name == "TEB_WSNOW_ROOF") then + read(string(indx+1:),*)TEB_WSNOW_ROOF_TBL(1:icate) + else if (name == "TEB_TSNOW_ROOF") then + read(string(indx+1:),*)TEB_TSNOW_ROOF_TBL(1:icate) + else if (name == "TEB_RSNOW_ROOF") then + read(string(indx+1:),*)TEB_RSNOW_ROOF_TBL(1:icate) + else if (name == "TEB_ASNOW_ROOF") then + read(string(indx+1:),*)TEB_ASNOW_ROOF_TBL(1:icate) + else if (name == "TEB_ESNOW_ROOF") then + read(string(indx+1:),*)TEB_ESNOW_ROOF_TBL(1:icate) + else if (name == "TEB_TSSNOW_ROOF") then + read(string(indx+1:),*)TEB_TSSNOW_ROOF_TBL(1:icate) + else if (name == "TEB_WSNOW_ROAD") then + read(string(indx+1:),*)TEB_WSNOW_ROAD_TBL(1:icate) + else if (name == "TEB_TSNOW_ROAD") then + read(string(indx+1:),*)TEB_TSNOW_ROAD_TBL(1:icate) + else if (name == "TEB_RSNOW_ROAD") then + read(string(indx+1:),*)TEB_RSNOW_ROAD_TBL(1:icate) + else if (name == "TEB_ASNOW_ROAD") then + read(string(indx+1:),*)TEB_ASNOW_ROAD_TBL(1:icate) + else if (name == "TEB_ESNOW_ROAD") then + read(string(indx+1:),*)TEB_ESNOW_ROAD_TBL(1:icate) + else if (name == "TEB_TSSNOW_ROAD") then + read(string(indx+1:),*)TEB_TSSNOW_ROAD_TBL(1:icate) + else if (name == "TEB_SHGC") then + read(string(indx+1:),*)TEB_SHGC_TBL(1:icate) + else if (name == "TEB_SHGC_SH") then + read(string(indx+1:),*)TEB_SHGC_SH_TBL(1:icate) + else if (name == "TEB_U_WIN") then + read(string(indx+1:),*)TEB_U_WIN_TBL(1:icate) + else if (name == "TEB_LSHADE") then + read(string(indx+1:),*)TEB_LSHADE_TBL(1:icate) + TEB_LSHADE_INITED = .TRUE. + else if (name == "TEB_ZSHADE") then + read(string(indx+1:),*)TEB_ZSHADE_TBL(1:icate) + else if (name == "TEB_CBEM") then + read(string(indx+1:),*)TEB_CBEM_TBL(1:icate) + else if (name == "TEB_FLOOR_HEIGHT") then + read(string(indx+1:),*)TEB_FLOOR_HEIGHT_TBL(1:icate) + else if (name == "TEB_CH_BEM") then + read(string(indx+1:),*)TEB_CH_BEM_TBL(1:icate) + else if (name == "TEB_ROUGH_ROOF") then + read(string(indx+1:),*)TEB_ROUGH_ROOF_TBL(1:icate) + else if (name == "TEB_ROUGH_WALL") then + read(string(indx+1:),*)TEB_ROUGH_WALL_TBL(1:icate) + else if (name == "TEB_PAR_RD_IRRIG") then + read(string(indx+1:),*)TEB_PAR_RD_IRRIG_TBL(1:icate) + TEB_PAR_RD_IRRIG_INITED = .TRUE. + else if (name == "TEB_RD_START_MONTH") then + read(string(indx+1:),*)TEB_RD_START_MONTH_TBL(1:icate) + else if (name == "TEB_RD_END_MONTH") then + read(string(indx+1:),*)TEB_RD_END_MONTH_TBL(1:icate) + else if (name == "TEB_RD_START_HOUR") then + read(string(indx+1:),*)TEB_RD_START_HOUR_TBL(1:icate) + else if (name == "TEB_RD_END_HOUR") then + read(string(indx+1:),*)TEB_RD_END_HOUR_TBL(1:icate) + else if (name == "TEB_RD_24H_IRRIG") then + read(string(indx+1:),*)TEB_RD_24H_IRRIG_TBL(1:icate) + else if (name == "TEB_EMIS_PANEL") then + read(string(indx+1:),*)TEB_EMIS_PANEL_TBL(1:icate) + else if (name == "TEB_ALB_PANEL") then + read(string(indx+1:),*)TEB_ALB_PANEL_TBL(1:icate) + else if (name == "TEB_EFF_PANEL") then + read(string(indx+1:),*)TEB_EFF_PANEL_TBL(1:icate) + else if (name == "TEB_FRAC_PANEL") then + read(string(indx+1:),*)TEB_FRAC_PANEL_TBL(1:icate) + else if (name == "TEB_RESIDENTIAL") then + read(string(indx+1:),*)TEB_RESIDENTIAL_TBL(1:icate) + else if (name == "TEB_DT_RES") then + read(string(indx+1:),*)TEB_DT_RES_TBL(1:icate) + else if (name == "TEB_DT_OFF") then + read(string(indx+1:),*)TEB_DT_OFF_TBL(1:icate) + else if (name == "TEB_T_ROAD") then + read(string(indx+1:),*)TEB_T_ROAD_TBL(1:teb_num_floor_layers) + else if (name == "TEB_T_ROOF") then + read(string(indx+1:),*)TEB_T_ROOF_TBL(1:teb_num_floor_layers) + else if (name == "TEB_T_WALL_A") then + read(string(indx+1:),*)TEB_T_WALL_A_TBL(1:teb_num_floor_layers) + else if (name == "TEB_T_WALL_B") then + read(string(indx+1:),*)TEB_T_WALL_B_TBL(1:teb_num_floor_layers) + ! ROOF + else if (name == "TEB ROOF PARAMETERS") then + + TEBROOFLOOP : do + read(11,'(A512)', iostat=iostatus) string + if (string(1:1) == "#") cycle TEBROOFLOOP + if (trim(string) == "") cycle TEBROOFLOOP + if (string == "END TEB ROOF PARAMETERS") exit TEBROOFLOOP + read(string, *) k + teb_numroof_tbl(k) = teb_numroof_tbl(k) + 1 + read(string, *) k, TEB_HC_ROOF_TBL(teb_numroof_tbl(k),k), & + TEB_TC_ROOF_TBL(teb_numroof_tbl(k),k), & + TEB_D_ROOF_TBL(teb_numroof_tbl(k),k) + enddo TEBROOFLOOP + ! ROAD + else if (name == "TEB ROAD PARAMETERS") then + + TEBROADLOOP : do + read(11,'(A512)', iostat=iostatus) string + if (string(1:1) == "#") cycle TEBROADLOOP + if (trim(string) == "") cycle TEBROADLOOP + if (string == "END TEB ROAD PARAMETERS") exit TEBROADLOOP + read(string, *) k + teb_numroad_tbl(k) = teb_numroad_tbl(k) + 1 + read(string, *) k, TEB_HC_ROAD_TBL(teb_numroad_tbl(k),k), & + TEB_TC_ROAD_TBL(teb_numroad_tbl(k),k), & + TEB_D_ROAD_TBL(teb_numroad_tbl(k),k) + enddo TEBROADLOOP + ! WALL + else if (name == "TEB WALL PARAMETERS") then + + TEBWALLLOOP : do + read(11,'(A512)', iostat=iostatus) string + if (string(1:1) == "#") cycle TEBWALLLOOP + if (trim(string) == "") cycle TEBWALLLOOP + if (string == "END TEB WALL PARAMETERS") exit TEBWALLLOOP + read(string, *) k + teb_numwall_tbl(k) = teb_numwall_tbl(k) + 1 + read(string, *) k, TEB_HC_WALL_TBL(teb_numwall_tbl(k),k), & + TEB_TC_WALL_TBL(teb_numwall_tbl(k),k), & + TEB_D_WALL_TBL(teb_numwall_tbl(k),k) + enddo TEBWALLLOOP + ! FLOOR + else if (name == "TEB FLOOR PARAMETERS") then + + TEBFLOORLOOP : do + read(11,'(A512)', iostat=iostatus) string + if (string(1:1) == "#") cycle TEBFLOORLOOP + if (trim(string) == "") cycle TEBFLOORLOOP + if (string == "END TEB FLOOR PARAMETERS") exit TEBFLOORLOOP + read(string, *) k + teb_numfloor_tbl(k) = teb_numfloor_tbl(k) + 1 + read(string, *) k, TEB_HC_FLOOR_TBL(teb_numfloor_tbl(k),k), & + TEB_TC_FLOOR_TBL(teb_numfloor_tbl(k),k), & + TEB_D_FLOOR_TBL(teb_numfloor_tbl(k),k) + enddo TEBFLOORLOOP + +! END: TEB + else CALL wrf_error_fatal('URBPARM.TBL: Unrecognized NAME = "'//trim(name)//'" in Subr URBAN_PARAM_INIT') endif @@ -2329,6 +3165,114 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & CLOSE(11) + ! TEB: check whether all required parameters are initialized + if (SF_URBAN_PHYSICS == 4 ) then + ! TODO some parameters are not required if certain options are disabled + ! -> disable check in that case, otherwise user is forced to define fake values + call check_init_real_array(TEB_ALB_ROOF_TBL, "TEB_ALB_ROOF") + call check_init_real_array(TEB_ALB_ROAD_TBL, "TEB_ALB_ROAD") + call check_init_real_array(TEB_ALB_WALL_TBL, "TEB_ALB_WALL") + call check_init_real_array(TEB_EMIS_ROOF_TBL, "TEB_EMIS_ROOF") + call check_init_real_array(TEB_EMIS_ROAD_TBL, "TEB_EMIS_ROAD") + call check_init_real_array(TEB_EMIS_WALL_TBL, "TEB_EMIS_WALL") + call check_init_real_array(TEB_BLD_HEIGHT_TBL, "TEB_BLD_HEIGHT") + call check_init_real_array(TEB_FRAC_BLD_TBL, "TEB_FRAC_BLD") + call check_init_real_array(TEB_WALL_O_TBL, "TEB_WALL_O_HOR") + call check_init_real_array(TEB_H_TRAFFIC_TBL, "TEB_H_TRAFFIC") + call check_init_real_array(TEB_LE_TRAFFIC_TBL, "TEB_LE_TRAFFIC") + call check_init_real_array(TEB_H_INDUSTRY_TBL, "TEB_H_INDUSTRY") + call check_init_real_array(TEB_LE_INDUSTRY_TBL, "TEB_LE_INDUSTRY") + call check_init_real_array(TEB_Z0_TOWN_TBL, "TEB_Z0_TOWN") + call check_init_integer_array(TEB_OPT_ZOH_TBL, "TEB_OPT_ZOH") + call check_init_integer_array(TEB_SNOW_SCH_ROOF_TBL, "TEB_SNOW_SCH_ROOF") + call check_init_integer_array(TEB_SNOW_SCH_ROAD_TBL, "TEB_SNOW_SCH_ROAD") + call check_init_real_array(TEB_T_FLOOR_TBL, "TEB_T_FLOOR") + call check_init_real_array(TEB_T_MASS_TBL, "TEB_T_MASS") + call check_init_real_array(TEB_TI_BLD_TBL, "TEB_TI_BLD") + call check_init_real_array(TEB_T_CANYON_TBL, "TEB_T_CANYON") + call check_init_real_array(TEB_Q_CANYON_TBL, "TEB_Q_CANYON") + call check_init_real_array(TEB_T_WIN1_TBL, "TEB_T_WIN1") + call check_init_real_array(TEB_T_WIN2_TBL, "TEB_T_WIN2") + call check_init_real_array(TEB_QI_BLD_TBL, "TEB_QI_BLD") + call check_init_real_array(TEB_WS_ROOF_TBL, "TEB_WS_ROOF") + call check_init_real_array(TEB_WS_ROAD_TBL, "TEB_WS_ROAD") + call check_init_real_array(TEB_WSNOW_ROOF_TBL, "TEB_WSNOW_ROOF") + call check_init_real_array(TEB_TSNOW_ROOF_TBL, "TEB_TSNOW_ROOF") + call check_init_real_array(TEB_RSNOW_ROOF_TBL, "TEB_RSNOW_ROOF") + call check_init_real_array(TEB_ASNOW_ROOF_TBL, "TEB_ASNOW_ROOF") + call check_init_real_array(TEB_ESNOW_ROOF_TBL, "TEB_ESNOW_ROOF") + call check_init_real_array(TEB_TSSNOW_ROOF_TBL, "TEB_TSSNOW_ROOD") + call check_init_real_array(TEB_WSNOW_ROAD_TBL, "TEB_WSNOW_ROAD") + call check_init_real_array(TEB_TSNOW_ROAD_TBL, "TEB_TSNOW_ROAD") + call check_init_real_array(TEB_RSNOW_ROAD_TBL, "TEB_RSNOW_ROAD") + call check_init_real_array(TEB_ASNOW_ROAD_TBL, "TEB_ASNOW_ROAD") + call check_init_real_array(TEB_ESNOW_ROAD_TBL, "TEB_ESNOW_ROAD") + call check_init_real_array(TEB_TSSNOW_ROAD_TBL, "TEB_TSSNOW_ROAD") + call check_init_logical(TEB_SOLAR_PANEL_INITED, "TEB_SOLAR_PANEL") + call check_init_logical(TEB_LGARDEN_INITED, "TEB_LGARDEN") + call check_init_logical(TEB_GREENROOF_INITED, "TEB_GREENROOF") + call check_init_logical(TEB_USE_NOAH_INITED, "TEB_USE_NOAH") + call check_init_logical(TEB_USE_ZGARDEN_FROM_TBL_INITED, "TEB_USE_ZGARDEN_FROM_TBL") + call check_init_string_array(TEB_HROAD_DIR_TBL, "TEB_HROAD_DIR") + call check_init_real_array(TEB_ZROAD_DIR_TBL, "TEB_ZROAD_DIR") + call check_init_string_array(TEB_WALL_OPT_TBL, "TEB_WALL_OPT") + call check_init_real_array(TEB_ZGARDEN_TBL, "TEB_ZGARDEN") + call check_init_real_array(TEB_FRAC_GR_TBL, "TEB_FRAC_GR") + call check_init_string_array(TEB_COOL_COIL_TBL, "TEB_COOL_COIL") + call check_init_string_array(TEB_HEAT_COIL_TBL, "TEB_HEAT_COIL") + call check_init_real_array(TEB_F_WATER_COND_TBL, "TEB_F_WATER_COND") + call check_init_string_array(TEB_HNATVENT_TBL, "TEB_HNATVENT") + call check_init_real_array(TEB_ZNATVENT_TBL, "TEB_ZNATVENT") + call check_init_real_array(TEB_F_WASTE_CAN_TBL, "TEB_F_WASTE_CAN") + call check_init_real_array(TEB_QIN_TBL, "TEB_QIN") + call check_init_real_array(TEB_QIN_FRAD_TBL, "TEB_QIN_FRAD") + call check_init_real_array(TEB_QIN_FLAT_TBL, "TEB_QIN_FLAT") + call check_init_real_array(TEB_GR_TBL, "TEB_GR") + call check_init_real_array(TEB_EFF_HEAT_TBL, "TEB_EFF_HEAT") + call check_init_real_array(TEB_INF_TBL, "TEB_INF") + call check_init_real_array(TEB_TCOOL_TARGET_TBL, "TEB_TCOOL_TARGET") + call check_init_real_array(TEB_THEAT_TARGET_TBL, "TEB_THEAT_TARGET") + call check_init_real_array(TEB_HR_TARGET_TBL, "TEB_HR_TARGET") + call check_init_real_array(TEB_V_VENT_TBL, "TEB_V_VENT") + call check_init_real_array(TEB_CAP_SYS_HEAT_TBL, "TEB_CAP_SYS_HEAT") + call check_init_real_array(TEB_CAP_SYS_RAT_TBL, "TEB_CAP_SYS_RAT") + call check_init_real_array(TEB_T_ADP_TBL, "TEB_T_ADP") + call check_init_real_array(TEB_M_SYS_RAT_TBL, "TEB_M_SYS_RAT") + call check_init_real_array(TEB_COP_RAT_TBL, "TEB_COP_RAT") + call check_init_real_array(TEB_SHGC_TBL, "TEB_SHGC") + call check_init_real_array(TEB_SHGC_SH_TBL, "TEB_SHGC_SH") + call check_init_real_array(TEB_U_WIN_TBL, "TEB_U_WIN") + call check_init_logical(TEB_LSHADE_INITED, "TEB_LSHADE") + call check_init_real_array(TEB_ZSHADE_TBL, "TEB_ZSHADE") + call check_init_string_array(TEB_CBEM_TBL, "TEB_CBEM") + call check_init_real_array(TEB_FLOOR_HEIGHT_TBL, "TEB_FLOOR_HEIGHT") + call check_init_string_array(TEB_CH_BEM_TBL, "TEB_CH_BEM") + call check_init_real_array(TEB_ROUGH_ROOF_TBL, "TEB_ROUGH_ROOF") + call check_init_real_array(TEB_ROUGH_WALL_TBL, "TEB_ROUGH_WALL") + call check_init_logical(TEB_PAR_RD_IRRIG_INITED, "TEB_PAR_RD_IRRIG") + call check_init_real_array(TEB_RD_START_MONTH_TBL, "TEB_RD_START_MONTH") + call check_init_real_array(TEB_RD_END_MONTH_TBL, "TEB_RD_END_MONTH") + call check_init_real_array(TEB_RD_START_HOUR_TBL, "TEB_RD_START_HOUR") + call check_init_real_array(TEB_RD_END_HOUR_TBL, "TEB_RD_END_HOUR") + call check_init_real_array(TEB_RD_24H_IRRIG_TBL, "TEB_RD_24H_IRRIG") + call check_init_real_array(TEB_EMIS_PANEL_TBL, "TEB_EMIS_PANEL") + call check_init_real_array(TEB_ALB_PANEL_TBL, "TEB_ALB_PANEL") + call check_init_real_array(TEB_EFF_PANEL_TBL, "TEB_EFF_PANEL") + call check_init_real_array(TEB_FRAC_PANEL_TBL, "TEB_FRAC_PANEL") + call check_init_real_array(TEB_RESIDENTIAL_TBL, "TEB_RESIDENTIAL") + call check_init_real_array(TEB_DT_RES_TBL, "TEB_DT_RES") + call check_init_real_array(TEB_DT_OFF_TBL, "TEB_DT_OFF") + call check_init_real_array(TEB_T_ROAD_TBL, "TEB_T_ROAD") + call check_init_real_array(TEB_T_ROOF_TBL, "TEB_T_ROOF") + call check_init_real_array(TEB_T_WALL_A_TBL, "TEB_T_WALL_A") + call check_init_real_array(TEB_T_WALL_B_TBL, "TEB_T_WALL_B") + + ! initialize constants for TEB + CALL INI_CSTS + CALL INIT_SURFCONSPHY + end if + + IF (SF_URBAN_PHYSICS /= 4) THEN ! Assign a few table values that do not need to come from URBPARM.TBL Z0HB_TBL = 0.1 * Z0B_TBL @@ -2404,6 +3348,22 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & deallocate(roof_width) deallocate(road_width) + END IF + + !for TEB + IF (SF_URBAN_PHYSICS == 4) THEN + TEB_Z0H_TBL(1) = 'MASC95' ! Mascart et al 1995 + TEB_Z0H_TBL(2) = 'BRUT82' ! Brustaert 1982 + TEB_Z0H_TBL(3) = 'KAND07' ! Kanda 2007 + TEB_SNOW_ROOF_TBL(1) = 'NONE' + TEB_SNOW_ROOF_TBL(2) = 'D95 ' + TEB_SNOW_ROOF_TBL(3) = '1-L ' + TEB_SNOW_ROAD_TBL(1) = 'NONE' + TEB_SNOW_ROAD_TBL(2) = 'D95 ' + TEB_SNOW_ROAD_TBL(3) = '1-L ' + ENDIF +!end TEB + END SUBROUTINE urban_param_init !=========================================================================== ! @@ -2452,7 +3412,33 @@ SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, A_E_BEP,B_U_BEP,B_V_BEP, & ! inout multi-layer urban B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP, & ! inout multi-layer urban DL_U_BEP,SF_BEP,VL_BEP, & ! inout multi-layer urban - FRC_URB2D, UTYPE_URB2D) ! inout + FRC_URB2D, UTYPE_URB2D, & ! inout + +!for TEB + teb_num_snow_layers, & ! in + teb_num_road_layers, & ! in + teb_num_roof_layers, & ! in + teb_num_wall_layers, & ! in + teb_num_floor_layers, & ! in + TEB_TI_BLD, & ! inout + TEB_T_CANYON, TEB_Q_CANYON, & ! inout + TEB_WS_ROOF, TEB_WS_ROAD, & ! inout + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & ! inout + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & ! inout + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & ! inout + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & ! inout + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & ! inout + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & ! inout + TEB_T_WIN1, TEB_T_WIN2, & ! inout + TEB_AUX_MAX, TEB_THER_PRODC_DAY, TEB_QI_BLD, & ! inout + TEB_T_FLOOR, TEB_T_MASS, & ! inout + TEB_T_ROAD, TEB_T_ROOF, & ! inout + TEB_T_WALL_A, TEB_T_WALL_B, & ! inout + TEB_TSK_RURAL & ! inout + ) +! end TEB + + IMPLICIT NONE INTEGER, INTENT(IN) :: ISURBAN, sf_urban_physics @@ -2555,6 +3541,42 @@ SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, REAL, DIMENSION(ims:ime, kms:kme,jms:jme),INTENT(INOUT) :: SF_BEP REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DL_U_BEP ! +!for TEB + INTEGER, INTENT(IN) :: teb_num_snow_layers + INTEGER, INTENT(IN) :: teb_num_road_layers + INTEGER, INTENT(IN) :: teb_num_roof_layers + INTEGER, INTENT(IN) :: teb_num_wall_layers + INTEGER, INTENT(IN) :: teb_num_floor_layers + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TI_BLD + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_CANYON + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_Q_CANYON + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROOF + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROAD + REAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROOF + REAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROOF + REAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROOF + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROOF + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROOF + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROOF + REAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROAD + REAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROAD + REAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROAD + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROAD + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROAD + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROAD + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN1 + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN2 + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_AUX_MAX + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_THER_PRODC_DAY + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_QI_BLD + REAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_FLOOR + REAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_MASS + REAL, DIMENSION(ims:ime, teb_num_road_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROAD + REAL, DIMENSION(ims:ime, teb_num_roof_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROOF + REAL, DIMENSION(ims:ime, teb_num_wall_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_A + REAL, DIMENSION(ims:ime, teb_num_wall_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_B + REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSK_RURAL +!end TEB REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UTYPE_URB2D INTEGER :: UTYPE_URB @@ -2880,6 +3902,54 @@ SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, DL_U_BEP(I,K,J)=0. END DO ENDIF !sf_urban_physics=2 + +!for TEB + IF( sf_urban_physics .EQ. 4 )THEN + + TEB_TSK_RURAL(I,J) = TSURFACE0_URB(I,J) + + IF (UTYPE_URB2D(I,J) > 0) THEN + TEB_TI_BLD(I,J) = TEB_TI_BLD_TBL(UTYPE_URB2D(I,J)) ! indoor air temperature + TEB_T_CANYON(I,J) = TEB_T_CANYON_TBL(UTYPE_URB2D(I,J)) ! indoor air temperature + TEB_Q_CANYON(I,J) = TEB_Q_CANYON_TBL(UTYPE_URB2D(I,J)) ! indoor air temperature + + TEB_WS_ROOF(I,J) = TEB_WS_ROOF_TBL(UTYPE_URB2D(I,J)) ! roof water reservoir(kg/m2) + TEB_WS_ROAD(I,J) = TEB_WS_ROAD_TBL(UTYPE_URB2D(I,J)) ! road water reservoir(kg/m2) + TEB_ASNOW_ROOF(I,J) = TEB_ASNOW_ROOF_TBL(UTYPE_URB2D(I,J)) ! snow albedo + TEB_ESNOW_ROOF(I,J) = TEB_ESNOW_ROOF_TBL(UTYPE_URB2D(I,J)) ! snow emissivity + TEB_TSSNOW_ROOF(I,J) = TEB_TSSNOW_ROOF_TBL(UTYPE_URB2D(I,J)) ! snow surface temperature + TEB_ASNOW_ROAD(I,J) = TEB_ASNOW_ROAD_TBL(UTYPE_URB2D(I,J)) ! snow albedo + TEB_ESNOW_ROAD(I,J) = TEB_ESNOW_ROAD_TBL(UTYPE_URB2D(I,J)) ! snow emissivity + TEB_TSSNOW_ROAD(I,J) = TEB_TSSNOW_ROAD_TBL(UTYPE_URB2D(I,J)) ! snow surface temperature + ! TODO allow to initialize each layer differently + DO K=1, teb_num_snow_layers + TEB_WSNOW_ROOF(I,K,J) = TEB_WSNOW_ROOF_TBL(UTYPE_URB2D(I,J)) ! snow layers reservoir + TEB_TSNOW_ROOF(I,K,J) = TEB_TSNOW_ROOF_TBL(UTYPE_URB2D(I,J)) ! snow layers temperature + TEB_RSNOW_ROOF(I,K,J) = TEB_RSNOW_ROOF_TBL(UTYPE_URB2D(I,J)) ! snow layers density + TEB_WSNOW_ROAD(I,K,J) = TEB_WSNOW_ROAD_TBL(UTYPE_URB2D(I,J)) ! snow layers reservoir + TEB_TSNOW_ROAD(I,K,J) = TEB_TSNOW_ROAD_TBL(UTYPE_URB2D(I,J)) ! snow layers temperature + TEB_RSNOW_ROAD(I,K,J) = TEB_RSNOW_ROAD_TBL(UTYPE_URB2D(I,J)) ! snow layers density + END DO + TEB_T_WIN1(I,J) = TEB_T_WIN1_TBL(UTYPE_URB2D(I,J)) ! External window temperature + TEB_T_WIN2(I,J) = TEB_T_WIN2_TBL(UTYPE_URB2D(I,J)) ! Internal window temperature + TEB_QI_BLD(I,J) = TEB_QI_BLD_TBL(UTYPE_URB2D(I,J)) ! Indoor air specific humidity [kg kg-1] + + TEB_AUX_MAX(I,J) = 5. ! Auxiliar variable for autosize calcs (NOT USED IN TEB!) + TEB_THER_PRODC_DAY(I,J) = 0. ! Present day integrated thermal production + ! of energy (J/m2 panel). zero value at start + + ! TODO allow to initialize each urban type differently + TEB_T_FLOOR(I,:,J) = TEB_T_FLOOR_TBL ! building floor temperature + TEB_T_MASS(I,:,J) = TEB_T_MASS_TBL ! building mass temperature + TEB_T_ROAD(I,:,J) = TEB_T_ROAD_TBL ! road temperature + TEB_T_ROOF(I,:,J) = TEB_T_ROOF_TBL ! roof temperature + TEB_T_WALL_A(I,:,J) = TEB_T_WALL_A_TBL ! wall temperature + TEB_T_WALL_B(I,:,J) = TEB_T_WALL_B_TBL ! wall temperature + ENDIF + + ENDIF +!end TEB + ENDIF !restart diff --git a/phys/module_surface_driver.F b/phys/module_surface_driver.F index 9481299df3..d05726f3d9 100644 --- a/phys/module_surface_driver.F +++ b/phys/module_surface_driver.F @@ -268,7 +268,49 @@ SUBROUTINE surface_driver( & & ,b_u_bep,b_v_bep,b_t_bep,b_q_bep & & ,sf_bep,vl_bep & & ,a_e_bep,b_e_bep,dlg_bep & - & ,dl_u_bep & + & ,dl_u_bep & + ! Optional urban Bep end + ! BEGIN: Optional urban for TEB + & ,teb_test_integration & + & ,teb_num_snow_layers & + & ,teb_num_road_layers & + & ,teb_num_roof_layers & + & ,teb_num_wall_layers & + & ,teb_num_floor_layers & + & ,TEB_TI_BLD & + & ,TEB_T_CANYON, TEB_Q_CANYON & + & ,TEB_WS_ROOF, TEB_WS_ROAD & + & ,TEB_WSNOW_ROOF, TEB_TSNOW_ROOF & + & ,TEB_RSNOW_ROOF, TEB_ASNOW_ROOF & + & ,TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF & + & ,TEB_WSNOW_ROAD, TEB_TSNOW_ROAD & + & ,TEB_RSNOW_ROAD, TEB_ASNOW_ROAD & + & ,TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD & + & ,TEB_T_WIN1, TEB_T_WIN2 & + & ,TEB_AUX_MAX, TEB_THER_PRODC_DAY, TEB_QI_BLD & + & ,TEB_T_FLOOR, TEB_T_MASS & + & ,TEB_T_ROAD, TEB_T_ROOF & + & ,TEB_T_WALL_A, TEB_T_WALL_B & + & ,TEB_HVAC_COOL, TEB_HVAC_HEAT & + & ,TEB_THER_PROD_PANEL, TEB_PHOT_PROD_PANEL & + & ,TEB_TSK_RURAL & + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + & ,TEB_INPUT_TA & + & ,TEB_INPUT_PS & + & ,TEB_INPUT_QA & + & ,TEB_INPUT_WIND & + & ,TEB_INPUT_DIR & + & ,TEB_INPUT_DIR_SW & + & ,TEB_INPUT_SCA_SW & + & ,TEB_INPUT_LW & + & ,TEB_INPUT_RAIN & + & ,TEB_INPUT_SNOW & + & ,TEB_INPUT_DIR_CO2 & + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + ! + ! END: Optional urban for TEB + ! & ,tsk_save & !for fractional seaice & ,cldfra & !ssib & ,sf_surface_mosaic,mosaic_cat,mosaic_cat_index & !danli mosaic @@ -1251,6 +1293,63 @@ SUBROUTINE surface_driver( & REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D !urban INTEGER, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UTYPE_URB2D !urban + ! BEGIN: Optional urban for TEB + INTEGER , INTENT(IN) :: teb_test_integration + INTEGER , INTENT(IN) :: teb_num_snow_layers + INTEGER , INTENT(IN) :: teb_num_road_layers + INTEGER , INTENT(IN) :: teb_num_roof_layers + INTEGER , INTENT(IN) :: teb_num_wall_layers + INTEGER , INTENT(IN) :: teb_num_floor_layers + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_Q_CANYON + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_WS_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_WSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_TSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_snow_layers, jms:jme), INTENT(INOUT) :: TEB_RSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ASNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_ESNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_TSSNOW_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN1 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_T_WIN2 + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_AUX_MAX + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_THER_PRODC_DAY + REAL, OPTIONAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT) :: TEB_QI_BLD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_FLOOR + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_MASS + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROAD + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_ROOF + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_A + REAL, OPTIONAL, DIMENSION(ims:ime, teb_num_floor_layers, jms:jme), INTENT(INOUT) :: TEB_T_WALL_B + REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: TEB_HVAC_COOL + REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: TEB_HVAC_HEAT + REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: TEB_THER_PROD_PANEL + REAL, DIMENSION(ims:ime,jms:jme), INTENT(OUT) :: TEB_PHOT_PROD_PANEL + REAL, DIMENSION(ims:ime,jms:jme), INTENT(INOUT) :: TEB_TSK_RURAL + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_TA + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_PS + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_QA + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_WIND + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_DIR + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_DIR_SW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_SCA_SW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_LW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_RAIN + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_SNOW + REAL, OPTIONAL, DIMENSION (ims:ime,jms:jme), INTENT(OUT) :: TEB_INPUT_DIR_CO2 + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + ! + ! END: Optional urban for TEB + REAL, DIMENSION( ims:ime, jms:jme ) :: PSIM_URB2D !urban local var REAL, DIMENSION( ims:ime, jms:jme ) :: PSIH_URB2D !urban local var REAL, DIMENSION( ims:ime, jms:jme ) :: GZ1OZ0_URB2D !urban local var @@ -2788,8 +2887,50 @@ SUBROUTINE surface_driver( & a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban a_e_bep,b_u_bep,b_v_bep, & !O multi-layer urban b_t_bep,b_q_bep,b_e_bep,dlg_bep, & !O multi-layer urban - dl_u_bep,sf_bep,vl_bep & !O multi-layer urban - ,sfcheadrt,INFXSRT, soldrain & + dl_u_bep,sf_bep,vl_bep, & !O multi-layer urban + ! BEGIN: Optional urban for TEB + XTIME, & !I + teb_test_integration, & !I + teb_num_snow_layers, & !I + teb_num_road_layers, & !I + teb_num_roof_layers, & !I + teb_num_wall_layers, & !I + teb_num_floor_layers, & !I + diffuse_frac, & !I + TEB_TI_BLD, & !H + TEB_T_CANYON,TEB_Q_CANYON, & !H + TEB_WS_ROOF, TEB_WS_ROAD, & !H + TEB_WSNOW_ROOF, TEB_TSNOW_ROOF, & !H + TEB_RSNOW_ROOF, TEB_ASNOW_ROOF, & !H + TEB_TSSNOW_ROOF, TEB_ESNOW_ROOF, & !H + TEB_WSNOW_ROAD, TEB_TSNOW_ROAD, & !H + TEB_RSNOW_ROAD, TEB_ASNOW_ROAD, & !H + TEB_TSSNOW_ROAD, TEB_ESNOW_ROAD, & !H + TEB_T_WIN1, TEB_T_WIN2, & !H + TEB_AUX_MAX, TEB_THER_PRODC_DAY, TEB_QI_BLD, & !H + TEB_T_FLOOR, TEB_T_MASS, & !H + TEB_T_ROAD, TEB_T_ROOF, & !H + TEB_T_WALL_A, TEB_T_WALL_B, & !H + TEB_TSK_RURAL, & !H + TEB_HVAC_COOL, TEB_HVAC_HEAT, & !O + TEB_THER_PROD_PANEL, TEB_PHOT_PROD_PANEL, & !O + ! + ! BEGIN: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + TEB_INPUT_TA, & + TEB_INPUT_PS, & + TEB_INPUT_QA, & + TEB_INPUT_WIND, & + TEB_INPUT_DIR, & + TEB_INPUT_DIR_SW, & + TEB_INPUT_SCA_SW, & + TEB_INPUT_LW, & + TEB_INPUT_RAIN, & + TEB_INPUT_SNOW, & + TEB_INPUT_DIR_CO2, & + ! END: INPUT quantities used for evaluating WRF-TEB with TEB OFFLINE + ! + ! END: Optional urban for TEB + sfcheadrt,INFXSRT, soldrain & ! ! FASDAS ! @@ -2920,6 +3061,25 @@ SUBROUTINE surface_driver( & ENDDO !urban ENDIF + ! TEB + IF(SF_URBAN_PHYSICS.eq.4) THEN + DO j=j_start(ij),j_end(ij) !urban + DO i=i_start(ij),i_end(ij) !urban + IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LOW_DENSITY_RESIDENTIAL .or. & + IVGTYP(I,J) == HIGH_DENSITY_RESIDENTIAL .or. IVGTYP(I,J) == HIGH_INTENSITY_INDUSTRIAL ) THEN + ! TH2_URB2D(I,J): temperature urb (not potential temp) + T2(I,J) = TH2_URB2D(i,j) + TH2(I,J) = T2(I,J)*((1.E5/PSFC(I,J))**RCP) + Q2(I,J) = Q2_URB2D(i,j) + U10(I,J) = U10_URB2D(I,J) !urban + V10(I,J) = V10_URB2D(I,J) !urban + END IF + ENDDO + ENDDO + END IF + !end TEB +! + !------------------------------------------------------------------ ELSE @@ -2927,6 +3087,11 @@ SUBROUTINE surface_driver( & ENDIF CASE (NOAHMPSCHEME) + + IF (SF_URBAN_PHYSICS .eq. 4) THEN + call wrf_error_fatal('TEB does not support NOAHMP') + ENDIF + IF (PRESENT(qv_curr) .AND. PRESENT(rainbl) .AND. & ! PRESENT(emiss) .AND. PRESENT(t2) .AND. & ! PRESENT(declin) .AND. PRESENT(coszen) .AND. & diff --git a/run/URBPARM_TEB.TBL b/run/URBPARM_TEB.TBL new file mode 100644 index 0000000000..3a4c5397de --- /dev/null +++ b/run/URBPARM_TEB.TBL @@ -0,0 +1,886 @@ +# The parameters in this table may vary greatly from city to city. +# The default values are probably not appropriate for any given city. +# Users should adapt these values based on the city they are working +# with. + +# Urban Parameters depending on Urban type +# USGS + +Number of urban categories: 3 + +# Main TEB options +############################################################################## + +# +# TEB_CBEM: Building Energy model 'DEF' or 'BEM' +# (sf_urban_physics=4) +# + +TEB_CBEM: BEM BEM BEM + +# +# TEB_HROAD_DIR: road direction option ('UNIF' = uniform roads, 'ORIE' = specified road orientation) +# (sf_urban_physics=4) +# + +TEB_HROAD_DIR: UNIF UNIF UNIF + +# +# TEB_WALL_OPT: Wall option ('UNIF' = uniform walls, 'TWO' = 2 opposite walls) +# (sf_urban_physics=4) +# + +TEB_WALL_OPT: UNIF UNIF UNIF + +# +# TEB_LGARDEN: Flag to use a vegetation scheme in gardens +# (sf_urban_physics=4) +# + +TEB_LGARDEN: .TRUE. .TRUE. .TRUE. + +# +# TEB_GREENROOF: Flag to use a green roofs scheme +# (sf_urban_physics=4) +# + +TEB_GREENROOF: .FALSE. .FALSE. .FALSE. + +# +# TEB_FRAC_GR: fraction of greenroofs on roofs +# (sf_urban_physics=4) +# + +TEB_FRAC_GR: 0.0 0.0 0.0 + +# +# TEB_SOLAR_PANEL: Flag to use a solar panels on roofs +# (sf_urban_physics=4) +# + +TEB_SOLAR_PANEL: .FALSE. .FALSE. .FALSE. + +# +# TEB_FRAC_PANEL: fraction of solar panels on roofs +# (sf_urban_physics=4) +# + +TEB_FRAC_PANEL: 0.0 0.0 0.0 + +# +# TEB_PAR_RD_IRRIG: Flag for road watering +# (sf_urban_physics=4) +# + +TEB_PAR_RD_IRRIG: .FALSE. .FALSE. .FALSE. + +# +# TEB_HNATVENT: flag to activate natural ventilation 'NONE', 'MANU', 'AUTO' +# (sf_urban_physics=4) +# + +TEB_HNATVENT: NONE NONE NONE + +# +# TEB_COOL_COIL: option for cooling device type ('IDEAL', 'DXCOIL', 'MinimalDX') +# (sf_urban_physics=4) +# + +TEB_COOL_COIL: MinimalDX MinimalDX MinimalDX + +# +# TEB_HEAT_COIL: option for heating device type ('IDEAL', 'FINCAP') +# (sf_urban_physics=4) +# + +TEB_HEAT_COIL: IDEAL IDEAL IDEAL + +# +# TEB_OPT_ZOH : TEB option for z0h roof & road +# [ 1: Mascart et al 1995, 2: Brustaert 1982, 3: Kanda 2007 ] +# (sf_urban_physics=4) + + TEB_OPT_ZOH: 3, 3, 3 + + +# Urban geometry +############################################################################## + +# +# TEB_Z0_TOWN : roughness length for momentum over town +# (sf_urban_physics=4) +# + + TEB_Z0_TOWN : 2.0, 2.0, 2.0 + +# +# TEB_FRAC_BLD : fraction of building (= roof fraction ). +# (sf_urban_physics=4) +# + + TEB_FRAC_BLD: 0.55, 0.55, 0.55 + +# +# TEB_USE_NOAH: whether to use Noah or TEB for vegetation +# (sf_urban_physics=4) +# + +TEB_USE_NOAH: .FALSE. + +# +# TEB_USE_ZGARDEN_FROM_TBL: whether to use TEB_ZGARDEN or the urban fraction from WRF +# (requires TEB_USE_NOAH=.FALSE.) +# (sf_urban_physics=4) +# + +TEB_USE_ZGARDEN_FROM_TBL: .TRUE. + +# +# TEB_ZGARDEN: fraction of GARDEN areas +# (requires TEB_USE_NOAH=.FALSE. and TEB_USE_ZGARDEN_FROM_TBL=.TRUE.) +# (sf_urban_physics=4) +# + +TEB_ZGARDEN: 0.10 0.10 0.10 + +# +# TEB_BLD_HEIGHT: Roof level (building height) [ m ] +# (sf_urban_physics=4) + +TEB_BLD_HEIGHT: 15.0, 15.0, 15.0 + +# +# TEB_WALL_O_HOR : wall surfuce. / hor. surface +# (sf_urban_physics=4) +# + + TEB_WALL_O_HOR: 1.4, 1.4, 1.4 + +# +# TEB_ZROAD_DIR: road direction (° from North, clockwise) +# (sf_urban_physics=4) +# + +TEB_ZROAD_DIR: 0. 0. 0. + +# Roof +############################################################################## + +# +# TEB_ALB_ROOF: Surface albedo of roof [ fraction ] +# (sf_urban_physics=4) +# + +TEB_ALB_ROOF: 0.15, 0.15, 0.15 + +# +# TEB_EMIS_ROOF: Surface emissivity of roof [ - ] +# (sf_urban_physics=4) +# + +TEB_EMIS_ROOF: 0.90, 0.90, 0.90 + + +TEB ROOF PARAMETERS: +# (sf_urban_physics=4) + +# urban heat thermal thickcness +# category capacity conductivity +# [index] [J m-3 K-1] [W/m K] [m] + 1 1580000. 1.15 0.001 + 1 1580000. 1.15 0.098 + 1 1580000. 1.15 0.132 + 1 1127845.62 0.095454545 0.098 + 1 52030. 0.03 0.001 + 2 1580000. 1.15 0.001 + 2 1580000. 1.15 0.098 + 2 1580000. 1.15 0.132 + 2 1127845.62 0.095454545 0.098 + 2 52030. 0.03 0.001 + 3 1580000. 1.15 0.001 + 3 1580000. 1.15 0.098 + 3 1580000. 1.15 0.132 + 3 1127845.62 0.095454545 0.098 + 3 52030. 0.03 0.001 + +END TEB ROOF PARAMETERS + +# Road +############################################################################## + +# +# TEB_ALB_ROAD: Surface albedo of ground (road) [ fraction ] +# (sf_urban_physics=4) +# + +TEB_ALB_ROAD: 0.08, 0.08, 0.08 + +# +# TEB_EMIS_ROAD: Surface emissivity of ground (road) [ - ] +# (sf_urban_physics=4) +# + +TEB_EMIS_ROAD: 0.95, 0.95, 0.95 + +TEB ROAD PARAMETERS: +# (sf_urban_physics=4) + +# urban heat thermal thickness +# category capacity conductivity +# [index] [J m-3 K-1] [W/m K] [m] + + 1 1740000. 0.82 0.001 + 1 1740000. 0.82 0.045296296 + 1 1989600. 1.976584 0.092592593 + 1 1640000. 0.5915493 0.27777778 + 1 1400000. 0.4000 0.83333333 + 2 1740000. 0.82 0.001 + 2 1740000. 0.82 0.045296296 + 2 1989600. 1.976584 0.092592593 + 2 1640000. 0.5915493 0.27777778 + 2 1400000. 0.4000 0.83333333 + 3 1740000. 0.82 0.001 + 3 1740000. 0.82 0.045296296 + 3 1989600. 1.976584 0.092592593 + 3 1640000. 0.5915493 0.27777778 + 3 1400000. 0.4000 0.83333333 + +END TEB ROAD PARAMETERS + +# Wall +############################################################################## + +# +# TEB_ALB_WALL: Surface albedo of building wall [ fraction ] +# (sf_urban_physics=1,2,3,4) +# + +TEB_ALB_WALL: 0.25, 0.25, 0.25 + +# +# TEB_EMIS_WALL: Surface emissivity of building wall [-] +# (sf_urban_physics=4) +# + +TEB_EMIS_WALL: 0.92, 0.92, 0.92 + +TEB WALL PARAMETERS: +# (sf_urban_physics=4) + +# urban heat thermal thickness +# category capacity conductivity +# [index] [J m-3 K-1] [W/m K] [m] + + 1 1580000. 1.15 0.001 + 1 1580000. 1.15 0.098 + 1 1580000. 1.15 0.132 + 1 1127845.62 0.095454545 0.098 + 1 52030. 0.03 0.001 + 2 1580000. 1.15 0.001 + 2 1580000. 1.15 0.098 + 2 1580000. 1.15 0.132 + 2 1127845.62 0.095454545 0.098 + 2 52030. 0.03 0.001 + 3 1580000. 1.15 0.001 + 3 1580000. 1.15 0.098 + 3 1580000. 1.15 0.132 + 3 1127845.62 0.095454545 0.098 + 3 52030. 0.03 0.001 + +END TEB WALL PARAMETERS + + +# Floor and internal mass +############################################################################## + +TEB FLOOR PARAMETERS: +# (sf_urban_physics=4) + +# urban heat thermal thickness +# category capacity conductivity +# [index] [J m-3 K-1] [W/m K] [m] + + 1 2016000. 1.95 0.001 + 1 2016000. 1.95 0.0064074074 + 1 2016000. 1.95 0.014814815 + 1 2016000. 1.95 0.044444444 + 1 2016000. 1.95 0.13333333 + 2 2016000. 1.95 0.001 + 2 2016000. 1.95 0.0064074074 + 2 2016000. 1.95 0.014814815 + 2 2016000. 1.95 0.044444444 + 2 2016000. 1.95 0.13333333 + 3 2016000. 1.95 0.001 + 3 2016000. 1.95 0.0064074074 + 3 2016000. 1.95 0.014814815 + 3 2016000. 1.95 0.044444444 + 3 2016000. 1.95 0.13333333 + +END TEB FLOOR PARAMETERS + + +# Anthropogenic heat fluxes +############################################################################## + +# +# TEB_H_TRAFFIC : anthropogenic sensible heat flux dues to traffic (w/m2) +# (sf_urban_physics=4) +# + + TEB_H_TRAFFIC: 8.00, 8.00, 8.00 + +# +# TEB_LE_TRAFFIC : anthropogenic latent heat flux dues to traffic (w/m2) +# (sf_urban_physics=4) +# + + TEB_LE_TRAFFIC: 2.00, 2.00, 2.00 + +# +# TEB_H_INDUSTRY : anthropogenic sensible heat flux dues to industry (w/m2) +# (sf_urban_physics=4) +# + + TEB_H_INDUSTRY: 0.00, 0.00, 0.00 + +# +# TEB_LE_INDUSTRY : latent sensible heat flux dues to industry (w/m2) +# (sf_urban_physics=4) +# + + TEB_LE_INDUSTRY: 0.00, 0.00, 0.00 + + +# Road watering +############################################################################## + +# +# TEB_RD_START_MONTH: start month for watering of roads(included) +# (sf_urban_physics=4) +# + +TEB_RD_START_MONTH: 0. 0. 0. + +# +# TEB_RD_END_MONTH: end month for watering of roads(included) +# (sf_urban_physics=4) +# + +TEB_RD_END_MONTH: 0. 0. 0. + +# +# TEB_RD_START_HOUR: start hour for watering of roads(included) +# (sf_urban_physics=4) +# + +TEB_RD_START_HOUR: 0. 0. 0. + +# +# TEB_RD_END_HOUR: end hour for watering of roads(excluded) +# (sf_urban_physics=4) +# + +TEB_RD_END_HOUR: 0. 0. 0. + +# +# TEB_RD_24H_IRRIG: 24h quantity of water used +# (sf_urban_physics=4) +# + +TEB_RD_24H_IRRIG: 0. 0. 0. + + +# Solar panels +############################################################################## + +# +# TEB_EMIS_PANEL: Emissivity of solar panel +# (sf_urban_physics=4) +# + +TEB_EMIS_PANEL: 0.9 0.9 0.9 + +# +# TEB_ALB_PANEL: albedo of solar panel +# (sf_urban_physics=4) +# + +TEB_ALB_PANEL: 0.1 0.1 0.1 + +# +# TEB_EFF_PANEL: Efficiency of solar panel +# (sf_urban_physics=4) +# + +TEB_EFF_PANEL: 0.14 0.14 0.14 + + +# Buildings' use information +############################################################################## + +# +# TEB_RESIDENTIAL: Fraction of residential use in buildings +# (sf_urban_physics=4) +# + +TEB_RESIDENTIAL: 1.0 1.0 1.0 + +# +# TEB_DT_RES: target temperature change when unoccupied (K) (residential buildings) +# (sf_urban_physics=4) +# + +TEB_DT_RES: 3. 3. 3. + +# +# TEB_DT_OFF: target temperature change when unoccupied (K) (office buildings) +# (sf_urban_physics=4) +# + +TEB_DT_OFF: 3. 3. 3. + + +# Parameters for Building Energy Module (BEM) +############################################################################## + +# Building configuration +############################################################################## + +# +# TEB_FLOOR_HEIGHT: Floor height (m) +# (sf_urban_physics=4) +# + +TEB_FLOOR_HEIGHT: 3.0 3.0 3.0 + +# +# TEB_INF: Infiltration flow rate [AC/H] +# (sf_urban_physics=4) +# + +TEB_INF: 0.5 0.5 0.5 + +# +# TEB_QIN: Internal heat gains [W m-2(floor)] +# (sf_urban_physics=4) +# + +TEB_QIN: 5.8 5.8 5.8 + +# +# TEB_QIN_FRAD: Radiant fraction of internal heat gains +# (sf_urban_physics=4) +# + +TEB_QIN_FRAD: 0.2 0.2 0.2 + +# +# TEB_QIN_FLAT: Latent franction of internal heat gains +# (sf_urban_physics=4) +# + +TEB_QIN_FLAT: 0.2 0.2 0.2 + + +# Windows +############################################################################## + +# +# TEB_GR: Glazing ratio +# (sf_urban_physics=4) +# + +TEB_GR: 0.1 0.1 0.1 + +# +# TEB_SHGC: window solar transmittance +# (sf_urban_physics=4) +# + +TEB_SHGC: 0.763 0.763 0.763 + +# +# TEB_U_WIN: window U-factor [K m W-2] +# (sf_urban_physics=4) +# + +TEB_U_WIN: 2.716 2.716 2.716 + + +# Shading devices +############################################################################## + +# +# TEB_LSHADE: Flag to use shading devices +# (sf_urban_physics=4) +# + +TEB_LSHADE: .FALSE. .FALSE. .FALSE. + +# +# TEB_ZSHADE: flag to activate shading devices (-> REAL for i/o 0. or 1) +# (sf_urban_physics=4) +# + +TEB_ZSHADE: 0. 0. 0. + +# +# TEB_SHGC_SH: window + shading solar heat gain coef. +# (sf_urban_physics=4) +# + +TEB_SHGC_SH: 0.025 0.025 0.025 + + +# Natural ventilation +############################################################################## + +# +# TEB_ZNATVENT: flag to describe surventilation system for i/o (0 for NONE, 1 for MANU and 2 for AUTO) +# (sf_urban_physics=4) +# + +TEB_ZNATVENT: 0. 0. 0. + + +# HVAC system +############################################################################## + +# +# TEB_V_VENT: Ventilation flow rate [AC/H] +# (sf_urban_physics=4) +# + +TEB_V_VENT: 0.0 0.0 0.0 + +# +# TEB_F_WATER_COND: fraction of evaporation for the condensers +# (sf_urban_physics=4) +# + +TEB_F_WATER_COND: 0. 0. 0. + +# +# TEB_F_WASTE_CAN: fraction of waste heat released into the canyon +# (sf_urban_physics=4) +# + +TEB_F_WASTE_CAN: 1.0 1.0 1.0 + + +# Internal target temperatures +############################################################################## + +# +# TEB_TCOOL_TARGET: Cooling setpoint of HVAC system [K] +# (sf_urban_physics=4) +# + +TEB_TCOOL_TARGET: 400.16 400.16 400.16 + +# +# TEB_THEAT_TARGET: Heating setpoint of HVAC system [K] +# (sf_urban_physics=4) +# + +TEB_THEAT_TARGET: 200.16 200.16 200.16 + +# +# TEB_HR_TARGET: Relative humidity setpoint +# (sf_urban_physics=4) +# + +TEB_HR_TARGET: 0.5 0.5 0.5 + + +# Heating system +############################################################################## + +# +# TEB_EFF_HEAT: Efficiency of the heating system +# (sf_urban_physics=4) +# + +TEB_EFF_HEAT: 0.9 0.9 0.9 + +# +# TEB_CAP_SYS_HEAT: Capacity of the heating system [W m-2(bld)] +# (sf_urban_physics=4) +# + +TEB_CAP_SYS_HEAT: 90.0 90.0 90.0 + + +# Cooling system +############################################################################## + +# +# TEB_CAP_SYS_RAT: Rated capacity of the cooling system [W m-2(bld)] +# (sf_urban_physics=4) +# + +TEB_CAP_SYS_RAT: 100.0 100.0 100.0 + +# +# TEB_T_ADP: Apparatus dewpoint temperature of the cooling coil [K] +# (sf_urban_physics=4) +# + +TEB_T_ADP: 285.66 285.66 285.66 + +# +# TEB_M_SYS_RAT: Rated HVAC mass flow rate [kg s-1 m-2(bld)] +# (sf_urban_physics=4) +# + +TEB_M_SYS_RAT: 0.0067 0.0067 0.0067 + +# +# TEB_COP_RAT: Rated COP of the cooling system +# (sf_urban_physics=4) +# + +TEB_COP_RAT: 2.5 2.5 2.5 + + +# Convection coefficients option +############################################################################## + +# +# TEB_CH_BEM: TEB option for building outside conv. coef +# (sf_urban_physics=4) +# + +TEB_CH_BEM: NO-DOE-2 NO-DOE-2 NO-DOE-2 + +# +# TEB_ROUGH_ROOF: roof roughness coef. +# (sf_urban_physics=4) +# + +TEB_ROUGH_ROOF: 1.52 1.52 1.52 + +# +# TEB_ROUGH_WALL: wall roughness coef. +# (sf_urban_physics=4) +# + +TEB_ROUGH_WALL: 1.52 1.52 1.52 + + +# Initialization for first time-step +############################################################################## +# These values correspond to the "capitoul" test case +# These values are NOT a function of urban type, but a function +# of the number of layers. + +# +# TEB_T_ROAD: Initial internal mass layers temperatures [K] +# This is currently NOT a function of urban type, but a function +# of the number of layers. +# (sf_urban_physics=4) +# + +TEB_T_ROAD: 274.07256 274.16793 274.45198 275.21494 277.50383 + +# +# TEB_T_ROOF: Initial internal mass layers temperatures [K] +# This is currently NOT a function of urban type, but a function +# of the number of layers. +# (sf_urban_physics=4) +# + +TEB_T_ROOF: 274.09930 276.95065 283.57500 289.55152 292.12402 + +# +# TEB_T_WALL_A: Initial internal mass layers temperatures [K] +## This is currently NOT a function of urban type, but a function +# of the number of layers. +# (sf_urban_physics=4) +# + +TEB_T_WALL_A: 274.09930 276.95065 283.57500 289.55152 292.12402 + +# +# TEB_T_WALL_B: Initial internal mass layers temperatures [K] +## This is currently NOT a function of urban type, but a function +# of the number of layers. +# (sf_urban_physics=4) +# + +TEB_T_WALL_B: 274.09930 276.95065 283.57500 289.55152 292.12402 + +# +# TEB_T_FLOOR: Initial floor layers temperatures [K] +# This is currently NOT a function of urban type, but a function +# of the number of layers. +# (sf_urban_physics=4) +# + +TEB_T_FLOOR: 292.12 291.89778 291.26111 289.48333 283.84017 + +# +# TEB_T_MASS: Initial internal mass layers temperatures [K] +# This is currently NOT a function of urban type, but a function +# of the number of layers. +# (sf_urban_physics=4) +# + +TEB_T_MASS: 292.15 292.15 292.15 292.15 291.84017 + +# +# TEB_TI_BLD: Initial indoor air temperature +# (sf_urban_physics=4) +# + +TEB_TI_BLD: 292.15 292.15 292.15 + +# +# TEB_T_CANYON: Initial indoor air temperature +# (sf_urban_physics=4) +# + +TEB_T_CANYON: 274.07050 274.07050 274.07050 + + +# +# TEB_T_WIN1: Initial external window temperature +# (sf_urban_physics=4) +# + +TEB_T_WIN1: 275. 275. 275. + +# +# TEB_T_WIN2: Initial internal window temperature +# (sf_urban_physics=4) +# + +TEB_T_WIN2: 292.15 292.15 292.15 + +# +# TEB_Q_CANYON: Initial indoor air temperature +# (sf_urban_physics=4) +# + +TEB_Q_CANYON: 0.0020888011 0.0020888011 0.0020888011 + +# +# TEB_QI_BLD: Initial indoor air specific humidity [kg kg-1] +# (sf_urban_physics=4) +# + +TEB_QI_BLD: 0.0068794074 0.0068794074 0.0068794074 + + +# +# TEB_SNOW_SCH_ROOF: Option for roof snow [ 1: none, 2: D95, 3: 1-L] +# (sf_urban_physics=4) +# + + TEB_SNOW_SCH_ROOF: 3, 3, 3 + +# +# TEB_SNOW_SCH_ROAD: Option for road snow [ 1: none, 2: D95, 3: 1-L] +# (sf_urban_physics=4) +# + + TEB_SNOW_SCH_ROAD: 3, 3, 3 + +# +# TEB_WSNOW_ROOF: Initial snow layers reservoir +# (sf_urban_physics=4) +# + +TEB_WSNOW_ROOF: 0. 0. 0. + +# +# TEB_TSNOW_ROOF: Initial snow layers temperature +# (sf_urban_physics=4) +# + +TEB_TSNOW_ROOF: 0. 0. 0. + +# +# TEB_RSNOW_ROOF: Initial snow layers density +# (sf_urban_physics=4) +# + +TEB_RSNOW_ROOF: 0. 0. 0. + +# +# TEB_ASNOW_ROOF: Initial snow albedo +# (sf_urban_physics=4) +# + +TEB_ASNOW_ROOF: 0.9 0.9 0.9 + +# +# TEB_ESNOW_ROOF: Initial snow emissivity +# (sf_urban_physics=4) +# + +TEB_ESNOW_ROOF: 0.91 0.91 0.91 + +# +# TEB_TSSNOW_ROOF: Initial snow surface temperature +# (sf_urban_physics=4) +# + +TEB_TSSNOW_ROOF: 0. 0. 0. + +# +# TEB_WSNOW_ROAD: Initial snow layers reservoir +# (sf_urban_physics=4) +# + +TEB_WSNOW_ROAD: 0. 0. 0. + +# +# TEB_TSNOW_ROAD: Initial snow layers temperature +# (sf_urban_physics=4) +# + +TEB_TSNOW_ROAD: 0. 0. 0. + +# +# TEB_RSNOW_ROAD: Initial snow layers density +# (sf_urban_physics=4) +# + +TEB_RSNOW_ROAD: 0. 0. 0. + +# +# TEB_ASNOW_ROAD: Initial snow albedo +# (sf_urban_physics=4) +# + +TEB_ASNOW_ROAD: 0.9 0.9 0.9 + +# +# TEB_ESNOW_ROAD: Initial snow emissivity +# (sf_urban_physics=4) +# + +TEB_ESNOW_ROAD: 0.91 0.91 0.91 + +# +# TEB_TSSNOW_ROAD: Initial snow surface temperature +# (sf_urban_physics=4) +# + +TEB_TSSNOW_ROAD: 0. 0. 0. + + +# +# TEB_WS_ROOF: Initial roof water reservoir(kg/m2) +# (sf_urban_physics=4) +# + +TEB_WS_ROOF: 0. 0. 0. + +# +# TEB_WS_ROAD: Initial road water reservoir(kg/m2) +# (sf_urban_physics=4) +# + +TEB_WS_ROAD: 0. 0. 0.