diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 50af9c465c..03f234edec 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -754,6 +754,7 @@ state real slope ij misc 1 - rdu "SLOP state real slp_azi ij misc 1 - rdu "SLP_AZI" "ELEVATION SLOPE AZIMUTH" "rad" state real shdmax ij misc 1 - i012rhd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "SHDMAX" "ANNUAL MAX VEG FRACTION" "" state real shdmin ij misc 1 - i012rhd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "SHDMIN" "ANNUAL MIN VEG FRACTION" "" +state real shdavg ij misc 1 - i012rhd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "SHDAVG" "ANNUAL AVG VEG FRACTION" "" state real snoalb ij misc 1 - i012rhd "SNOALB" "ANNUAL MAX SNOW ALBEDO IN FRACTION" "" state real toposoil ij misc 1 - i12 "SOILHGT" "ELEVATION OF LSM DATA" "m" state real landusef iuj misc 1 Z i012rdu "LANDUSEF" "LANDUSE FRACTION BY CATEGORY" "" @@ -807,14 +808,20 @@ state real DZR l em - Z r "DZR" state real DZB l em - Z r "DZB" "THICKNESSES OF WALL LAYERS" "m" state real DZG l em - Z r "DZG" "THICKNESSES OF ROAD LAYERS" "m" state real URB_PARAM i{urb}j misc 1 - i1 "URB_PARAM" "NUDAPT_NBSD Urban Parameters" "parameter" -state real LP_URB2D ij misc 1 - ir "BUILD_AREA_FRACTION" "BUILDING PLAN AREA DENSITY" "dimensionless" +state real LP_URB2D ij misc 1 - i01r "BUILD_AREA_FRACTION" "BUILDING PLAN AREA DENSITY" "dimensionless" state real HI_URB2D i{uhi}j misc 1 Z ir "HEIGHT_HISTOGRAMS" "DISTRIBUTION OF BUILDING HEIGHTS" "dimensionless" state real LB_URB2D ij misc 1 - ir "BUILD_SURF_RATIO" "BUILDING SURFACE AREA TO PLAN AREA RATIO" "dimensionless" state real HGT_URB2D ij misc 1 - ir "BUILD_HEIGHT" "AVERAGE BUILDING HEIGHT WEIGHTED BY BUILDING PLAN AREA" "m" -state real MH_URB2D ij misc 1 - ir "MH_URB2D" "Mean Building Height" "m" +state real MH_URB2D ij misc 1 - i01r "MH_URB2D" "Mean Building Height" "m" state real STDH_URB2D ij misc 1 - ir "STDH_URB2D" "Standard Deviation of Building Height" "m2" state real LF_URB2D i{udr}j misc 1 Z ir "LF_URB2D" "Frontal Area Index" "dimensionless" +state real ZD_URB2D ij misc 1 - i1 "ZD_URB2D" "Zero-plane Displacement" "m" +state real Z0_URB2D ij misc 1 - i01r "Z0_URB2D" "Roughness length for momentum" "m" +state real LF_URB2D_S ij misc 1 - i01r "LF_URB2D_S" "Frontal area index (no wind directional dependency)" "" +# AHE with month and hour dimension flattened to one dimension, Jan = (0:23), Feb = (24:47) +state real AHE i{m_hr}j misc 1 - i01r "AHE" "Anthropogenic heat emission" "W m-2" + # lsm State Variables state real SMOIS ilj - 1 Z i02rhd=(interp_mask_field:lu_index,iswater)u=(copy_fcnm) "SMOIS" "SOIL MOISTURE" "m3 m-3" @@ -2515,6 +2522,8 @@ rconfig integer ishallow namelist,physics 1 0 rconfig real convtrans_avglen_m namelist,physics 1 30 rh "convtrans_avglen_m" "averaging time for convective transport output variables (minutes)" "" rconfig integer num_land_cat namelist,physics 1 21 - "num_land_cat" "" "" rconfig integer use_wudapt_lcz namelist,physics 1 0 - "use_wudapt_lcz" "" "" +rconfig logical use_distributed_aerodynamics namelist,physics 1 .false. rh "use_distributed_aerodynamics" "" "" +rconfig integer distributed_ahe_opt namelist,physics 1 1 rh "distributed_ahe_opt" "AHE handling: 1=add to first level temperature tendency, 2=add to surface sensible heat flux" "" rconfig integer num_soil_cat namelist,physics 1 16 - "num_soil_cat" "" "" rconfig integer mp_zero_out namelist,physics 1 0 - "mp_zero_out" "microphysics fields set to zero 0=no action taken, 1=all fields but Qv, 2=all fields including Qv" "flag" rconfig real mp_zero_out_thresh namelist,physics 1 1.e-8 - "mp_zero_out_thresh" "minimum threshold for non-Qv moist fields, below are set to zero" "kg/kg" diff --git a/Registry/registry.dimspec b/Registry/registry.dimspec index 041bb2fefa..6761de8f7b 100644 --- a/Registry/registry.dimspec +++ b/Registry/registry.dimspec @@ -140,3 +140,5 @@ endif # Dimensions for PSU-DENG SCP dimspec nsh 2 constant=100 z nsh +# Dimensions for AHE +dimspec m_hr 2 constant=(0:287) z month_hour diff --git a/dyn_em/module_first_rk_step_part1.F b/dyn_em/module_first_rk_step_part1.F index f5eb26734d..7472ab395a 100644 --- a/dyn_em/module_first_rk_step_part1.F +++ b/dyn_em/module_first_rk_step_part1.F @@ -644,6 +644,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,BL_PBL_PHYSICS=config_flags%bl_pbl_physics & & ,SF_SURFACE_PHYSICS=config_flags%sf_surface_physics ,SH2O=grid%sh2o & & ,SHDMAX=grid%shdmax ,SHDMIN=grid%shdmin ,SMOIS=grid%smois & + & ,SHDAVG=grid%shdavg & & ,SMSTAV=grid%smstav ,SMSTOT=grid%smstot ,SNOALB=grid%snoalb & & ,SNOW=grid%snow ,SNOWC=grid%snowc ,SNOWH=grid%snowh & & ,SMCREL=grid%smcrel & @@ -810,6 +811,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,urban_map_zgrd = config_flags%urban_map_zgrd & !multi-layer urban & ,NUM_URBAN_HI=config_flags%num_urban_hi & !multi-layer urban & ,use_wudapt_lcz=config_flags%use_wudapt_lcz & !wudapt + & ,use_distributed_aerodynamics=config_flags%use_distributed_aerodynamics & !SLUCM & ,TSK_RURAL=grid%tsk_rural & !multi-layer urban & ,TRB_URB4D=grid%trb_urb4d,TW1_URB4D=grid%tw1_urb4d & !multi-layer urban & ,TW2_URB4D=grid%tw2_urb4d,TGB_URB4D=grid%tgb_urb4d & !multi-layer urban @@ -846,6 +848,7 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & & ,LB_URB2D=grid%lb_urb2d,HGT_URB2D=grid%hgt_urb2d & !multi-layer urban & ,MH_URB2D=grid%mh_urb2d,STDH_URB2D=grid%stdh_urb2d & !SLUCM & ,LF_URB2D=grid%lf_urb2d & + & ,lf_urb2d_s=grid%lf_urb2d_s, z0_urb2d=grid%z0_urb2d & & ,GMT=grid%gmt,XLAT=grid%xlat,XLONG=grid%xlong,JULDAY=grid%julday & & ,A_U_BEP=grid%a_u_bep,A_V_BEP=grid%a_v_bep,A_T_BEP=grid%a_t_bep & & ,A_Q_BEP=grid%a_q_bep & @@ -1155,6 +1158,8 @@ SUBROUTINE first_rk_step_part1 ( grid , config_flags & ! Bep changes end ! add tke_pbl, and turbulent fluxes & ,TKE_PBL=grid%tke_pbl,EL_PBL=grid%el_pbl,WU_TUR=grid%wu_tur & + & , gmt=grid%gmt, xtime=grid%xtime, julday=grid%julday, julyr=grid%julyr + & , ahe=grid%ahe, distributed_ahe_opt=grid%distributed_ahe_opt & ,WV_tur=grid%wv_tur,WT_tur=grid%wt_tur,WQ_tur=grid%wq_tur & & ,DISS_PBL=grid%diss_pbl,TPE_PBL=grid%tpe_pbl & & ,TKE_ADV=scalar(ims,kms,jms,P_tke_adv) & diff --git a/dyn_em/module_initialize_real.F b/dyn_em/module_initialize_real.F index 96629232bf..0604afb142 100644 --- a/dyn_em/module_initialize_real.F +++ b/dyn_em/module_initialize_real.F @@ -617,6 +617,15 @@ SUBROUTINE init_domain_rk ( grid & END IF END IF + IF (config_flags%use_distributed_aerodynamics) THEN + CALL wrf_message('Adding zero-plane displacement height to topography') + DO j = jts, MIN(jde - 1, jte) + DO i = its, MIN(ide - 1, ite) + IF (grid%zd_urb2d(i, j) > 0) grid%ht_gc(i, j) = grid%ht_gc(i, j) + grid%zd_urb2d(i, j) + END DO + END DO + END IF + ! Is there any vertical interpolation to do? The "old" data comes in on the correct ! vertical locations already. @@ -1330,6 +1339,11 @@ SUBROUTINE init_domain_rk ( grid & ims , ime , jms , jme , kms , kme , & its , ite , jts , jte , kts , kte ) + CALL monthly_avg ( grid%greenfrac , grid%shdavg , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte ) + ! The model expects the green-ness and vegetation fraction values to be in percent, not fraction. DO j = jts, MIN(jte,jde-1) @@ -1338,6 +1352,7 @@ SUBROUTINE init_domain_rk ( grid & grid%vegfra(i,j) = grid%vegfra(i,j) * 100. grid%shdmax(i,j) = grid%shdmax(i,j) * 100. grid%shdmin(i,j) = grid%shdmin(i,j) * 100. + grid%shdavg(i,j) = grid%shdavg(i,j) * 100. END DO END DO @@ -3089,6 +3104,16 @@ SUBROUTINE init_domain_rk ( grid & ! Split NUDAPT Urban Parameters + distributed_aerodynamics_if: IF (config_flags%sf_urban_physics == 1 .AND. config_flags%use_distributed_aerodynamics) THEN + DO j = jts , MIN(jde-1,jte) + DO i = its , MIN(ide-1,ite) + IF (grid%ivgtyp(i, j) == model_config_rec%isurban(grid%id)) THEN + grid%frc_urb2d(i, j) = MAX(0.1, MIN(0.9, 1 - grid%shdavg(i, j) / 100.)) + END IF + END DO + END DO + ELSE + IF ( ( config_flags%sf_urban_physics == 1 ) .OR. ( config_flags%sf_urban_physics == 2 ) .OR. ( config_flags%sf_urban_physics == 3 ) ) THEN DO j = jts , MIN(jde-1,jte) DO i = its , MIN(ide-1,ite) @@ -3145,6 +3170,8 @@ SUBROUTINE init_domain_rk ( grid & END DO END DO + END IF distributed_aerodynamics_if + END IF ! Adjustments for the seaice field PRIOR to the grid%tslb computations. This is @@ -7922,6 +7949,34 @@ SUBROUTINE monthly_min_max ( field_in , field_min , field_max , & END SUBROUTINE monthly_min_max +!--------------------------------------------------------------------- + + SUBROUTINE monthly_avg ( field_in , field_avg , & + ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte ) + + IMPLICIT NONE + + INTEGER , INTENT(IN) :: ids , ide , jds , jde , kds , kde , & + ims , ime , jms , jme , kms , kme , & + its , ite , jts , jte , kts , kte + + REAL , DIMENSION(ims:ime,12,jms:jme) , INTENT(IN) :: field_in + REAL , DIMENSION(ims:ime, jms:jme) , INTENT(OUT) :: field_avg + + ! Local vars + + INTEGER :: i , j + + DO j = jts , MIN(jde-1,jte) + DO i = its , MIN(ide-1,ite) + field_avg(i, j) = SUM(field_in(i, :, j)) / 12 + END DO + END DO + + END SUBROUTINE monthly_avg + !--------------------------------------------------------------------- SUBROUTINE monthly_interp_to_date ( field_in , date_str , field_out , & diff --git a/phys/module_pbl_driver.F b/phys/module_pbl_driver.F index fd2075f45b..24ab5bcf42 100644 --- a/phys/module_pbl_driver.F +++ b/phys/module_pbl_driver.F @@ -103,6 +103,9 @@ SUBROUTINE pbl_driver( & ,tke_adv,diss_adv,tpe_adv & ,pr_pbl,el_pbl & ,wu_tur,wv_tur,wt_tur,wq_tur & +! variables added for AHE + , gmt, xtime, julday, julyr, ahe & + , distributed_ahe_opt & ! variables for GBM PBL ,exch_tke, rthraten & ,a_e_bep,b_e_bep,dlg_bep,dl_u_bep & @@ -199,6 +202,7 @@ SUBROUTINE pbl_driver( & USE module_bl_fogdes USE module_wind_fitch #endif + use module_ra_gfdleta, only: cal_mon_day ! This driver calls subroutines for the PBL parameterizations. ! @@ -610,6 +614,11 @@ SUBROUTINE pbl_driver( & REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(OUT) :: EL_PBL + REAL, INTENT(IN) :: gmt, xtime + INTEGER, INTENT(IN) :: julday, julyr + REAL, OPTIONAL, DIMENSION( ims:ime, 0:287, jms:jme ), INTENT(IN) :: ahe + INTEGER, INTENT(IN) :: distributed_ahe_opt + REAL , INTENT(IN ) :: u_frame, & v_frame ! @@ -820,6 +829,7 @@ SUBROUTINE pbl_driver( & integer iu_bep,iurb,idiff real seamask,thsk,zzz,unew,vnew,tnew,qnew,umom,vmom REAL :: z0,z1,z2,w1,w2 + INTEGER :: ihour, jmonth, jday ! ! FASDAS ! @@ -2203,6 +2213,26 @@ SUBROUTINE pbl_driver( & ,ITS=its,ITE=ite,JTS=jts,JTE=jte,KTS=kts,KTE=kte) ENDIF + + IF (PRESENT(ahe)) THEN + call cal_mon_day(julday, julyr, jmonth, jday) + ihour = (jmonth - 1) * 24 + MOD(INT(gmt + xtime / 60.0), 24) + IF (distributed_ahe_opt == 1) THEN + DO j = jts, jte + DO i = its, ite + ! Volumetric heat capacity of air = 1200 J/(K m3) + RTHBLTEN(i, 1, j) = RTHBLTEN(i, 1, j) + ahe(i, ihour, j) / 1200 / DZ8W(i, 1, j) + END DO + END DO + ELSE IF (distributed_ahe_opt == 2) THEN + DO j = jts, jte + DO i = its, ite + HFX(i, j) = HFX(i, j) + ahe(i, ihour, j) + END DO + END DO + END IF + END IF + ENDDO !$OMP END PARALLEL DO diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 7b1487ed56..ac59b17ac9 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -3302,7 +3302,8 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & IF ( PRESENT( FRC_URB2D ) .AND. PRESENT( UTYPE_URB2D )) THEN CALL urban_param_init(DZR,DZB,DZG,num_soil_layers, & !urban - sf_urban_physics,config_flags%use_wudapt_lcz) !urban + sf_urban_physics,config_flags%use_wudapt_lcz, & + config_flags%use_distributed_aerodynamics) CALL urban_var_init(ISURBAN,TSK,TSLB,TMN,IVGTYP, & !urban @@ -3351,7 +3352,6 @@ SUBROUTINE bl_init(STEPBL,BLDT,DT,RUBLTEN,RVBLTEN,RTHBLTEN, & DL_U_BEP,SF_BEP,VL_BEP, & !multi-layer urban FRC_URB2D, UTYPE_URB2D,config_flags%use_wudapt_lcz) !urban - max_utype_urb2d = maxval(UTYPE_URB2D)*1.0 #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) max_utype_urb2d = wrf_dm_max_real(max_utype_urb2d) @@ -3445,7 +3445,8 @@ 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,config_flags%use_wudapt_lcz) + sf_urban_physics,config_flags%use_wudapt_lcz, & + config_flags%use_distributed_aerodynamics) CALL urban_var_init(ISURBAN,TSK,TSLB,TMN,IVGTYP, & !urban ims,ime,jms,jme,kms,kme,num_soil_layers, & !urban LCZ_1_TABLE,LCZ_2_TABLE,LCZ_3_TABLE,LCZ_4_TABLE, & !urban diff --git a/phys/module_sf_clm.F b/phys/module_sf_clm.F index 6d11ac7857..366cd333ad 100644 --- a/phys/module_sf_clm.F +++ b/phys/module_sf_clm.F @@ -59346,6 +59346,11 @@ subroutine clmdrv(zgcmxy ,forc_qxy ,ps ,forc_txy ,tsxy & real :: mh_urb,stdh_urb,lp_urb,hgt_urb,frc_urb,lb_urb,check real, dimension(4) :: lf_urb +! Distributed aerodynamics parameters + real :: lf_urb_s + real :: z0_urb + real :: vegfrac + logical, external :: wrf_dm_on_monitor ! ---------------------------------------------------------------------- @@ -60318,6 +60323,10 @@ subroutine clmdrv(zgcmxy ,forc_qxy ,ps ,forc_txy ,tsxy & enddo frc_urb = FRC_URB2D(I,J) check = 0. +! Distributed aerodynamics + lf_urb_s = 0 + z0_urb = 0 + vegfrac = 0 ! ! Call urban @@ -60346,7 +60355,8 @@ subroutine clmdrv(zgcmxy ,forc_qxy ,ps ,forc_txy ,tsxy & hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H TGRL_URB,SMR_URB,CMGR_URB, CHGR_URB, jmonth, & ! H DRELR_URB,DRELB_URB, & ! H - DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB) + DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, & + lf_urb_s, z0_urb, vegfrac) !sw-- TS_URB2D(I,J) = TS_URB diff --git a/phys/module_sf_noahdrv.F b/phys/module_sf_noahdrv.F index 5c7df673a2..e52d740854 100644 --- a/phys/module_sf_noahdrv.F +++ b/phys/module_sf_noahdrv.F @@ -110,6 +110,7 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d,& !RMS lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, & !H multi-layer urban mh_urb2d,stdh_urb2d,lf_urb2d, & !SLUCM + lf_urb2d_s, z0_urb2d, & !SLUCM th_phy,rho,p_phy,ust, & !I multi-layer urban gmt,julday,xlong,xlat, & !I multi-layer urban a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban @@ -595,6 +596,10 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & REAL :: lp_urb REAL :: hgt_urb REAL, DIMENSION(4) :: lf_urb +! Distributed aerodynamics parameters + REAL :: lf_urb_s + REAL :: z0_urb + REAL :: vegfrac ! Variables for multi-layer UCM (Martilli et al. 2002) REAL, OPTIONAL, INTENT(IN ) :: GMT INTEGER, OPTIONAL, INTENT(IN ) :: JULDAY @@ -655,6 +660,8 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: mh_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: stdh_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN) :: lf_urb2d + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lf_urb2d_s + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: z0_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature @@ -1416,6 +1423,10 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & if (I.eq.73.and.J.eq.125)THEN check = 1 end if +! Distributed aerodynamics + lf_urb_s = lf_urb2d_s(I, J) + z0_urb = z0_urb2d(I, J) + vegfrac = vegfra(I, J) / 100 ! ! Call urban CALL cal_mon_day(julian,julyr,jmonth,jday) @@ -1439,7 +1450,8 @@ SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, & hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H DRELR_URB,DRELB_URB, & ! H - DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB) + DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, & + lf_urb_s, z0_urb, vegfrac) #if 0 IF(IPRINT) THEN @@ -2345,7 +2357,7 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & myj,frpcpn, & SH2O,SNOWH, & !H U_PHY,V_PHY, & !I - SNOALB,SHDMIN,SHDMAX, & !I + SNOALB,SHDMIN,SHDMAX,SHDAVG, & !I SNOTIME, & !? ACSNOM,ACSNOW, & !O SNOPCX, & !O @@ -2415,6 +2427,7 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & urban_map_zgrd, & !I multi-layer urban num_urban_hi, & !I multi-layer urban use_wudapt_lcz, & !I wudapt + use_distributed_aerodynamics, & !I slucm tsk_rural_bep, & !H multi-layer urban trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, & !H multi-layer urban tlev_urb3d,qlev_urb3d, & !H multi-layer urban @@ -2430,6 +2443,7 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d,& !RMS lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, & !H multi-layer urban mh_urb2d,stdh_urb2d,lf_urb2d, & !SLUCM + lf_urb2d_s, z0_urb2d, & !SLUCM th_phy,rho,p_phy,ust, & !I multi-layer urban gmt,julday,xlong,xlat, & !I multi-layer urban a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban @@ -2607,6 +2621,7 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & VEGFRA, & SHDMIN, & SHDMAX, & + SHDAVG, & SNOALB, & GSW, & SWDOWN, & !added 10 jan 2007 @@ -2911,6 +2926,10 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & REAL :: lp_urb REAL :: hgt_urb REAL, DIMENSION(4) :: lf_urb +! Distributed aerodynamics parameters + REAL :: lf_urb_s + REAL :: z0_urb + REAL :: vegfrac ! Variables for multi-layer UCM (Martilli et al. 2002) REAL, OPTIONAL, INTENT(IN ) :: GMT INTEGER, OPTIONAL, INTENT(IN ) :: JULDAY @@ -2928,6 +2947,7 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & INTEGER, INTENT(IN ) :: urban_map_zgrd INTEGER, INTENT(IN ) :: NUM_URBAN_HI INTEGER, INTENT(IN ) :: use_wudapt_lcz + LOGICAL, INTENT(IN ) :: use_distributed_aerodynamics REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: tsk_rural_bep REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d @@ -2971,6 +2991,8 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: mh_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: stdh_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN) :: lf_urb2d + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lf_urb2d_s + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: z0_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature @@ -3721,7 +3743,12 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & ! UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial) ! this need to be changed in the mosaic danli - IF (use_wudapt_lcz == 1) THEN + IF (use_distributed_aerodynamics) THEN + IF (IVGTYP(I, J) == ISURBAN) THEN + UTYPE_URB = 2 + FRC_URB2D(I, J) = MAX(0.1, MIN(0.9, 1 - shdavg(I, J) / 100.)) + END IF + ELSE IF (use_wudapt_lcz == 1) THEN IF(IVGTYP(I,J)==ISURBAN) UTYPE_URB=5 IF(IVGTYP(I,J)==LCZ_1) UTYPE_URB=1 IF(IVGTYP(I,J)==LCZ_2) UTYPE_URB=2 @@ -3871,6 +3898,10 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & if (I.eq.73.and.J.eq.125)THEN check = 1 end if + ! Distributed aerodynamics + lf_urb_s = lf_urb2d_s(I, J) + z0_urb = z0_urb2d(I, J) + vegfrac = vegfra(I, J) / 100 ! ! Call urban CALL cal_mon_day(julian,julyr,jmonth,jday) @@ -3894,7 +3925,8 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H DRELR_URB,DRELB_URB, & ! H - DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB) + DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, & + lf_urb_s, z0_urb, vegfrac) #if 0 IF(IPRINT) THEN @@ -3936,7 +3968,9 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & ! Convert QSFC back to mixing ratio QSFC(I,J)= Q1/(1.0-Q1) UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s] + IF (.not. use_distributed_aerodynamics) THEN ZNT(I,J)= EXP(FRC_URB2D(I,J)*ALOG(ZNT_URB)+(1-FRC_URB2D(I,J))* ALOG(ZNT(I,J))) ! ADD BY DAN + END IF #if 0 IF(IPRINT)THEN @@ -4743,6 +4777,10 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & if (I.eq.73.and.J.eq.125)THEN check = 1 end if + ! Distributed aerodynamics + lf_urb_s = lf_urb2d_s(I, J) + z0_urb = z0_urb2d(I, J) + vegfrac = vegfra(I, J) / 100 ! ! Call urban CALL cal_mon_day(julian,julyr,jmonth,jday) @@ -4766,7 +4804,8 @@ SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, & hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H DRELR_URB,DRELB_URB, & ! H - DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB) + DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, & + lf_urb_s, z0_urb, vegfrac) #if 0 IF(IPRINT) THEN diff --git a/phys/module_sf_urban.F b/phys/module_sf_urban.F index 6a7eec58cc..58019acbe4 100644 --- a/phys/module_sf_urban.F +++ b/phys/module_sf_urban.F @@ -9,6 +9,8 @@ MODULE module_sf_urban #define WRITE_MESSAGE(M) call wrf_message( M ) #endif +USE module_model_constants, ONLY : piconst + !=============================================================================== ! Single-Layer Urban Canopy Model for WRF Noah-LSM ! Original Version: 2002/11/06 by Hiroyuki Kusaka @@ -88,6 +90,7 @@ MODULE module_sf_urban INTEGER :: IMP_SCHEME, IRI_SCHEME INTEGER :: alhoption ! anthropogenic latent heat option INTEGER :: groption ! anthropogenic latent heat option + LOGICAL :: distributed_aerodynamics_option REAL :: fgr ! green roof fraction REAL :: oasis ! urban oasis parameter REAL, DIMENSION(1:4) :: DZGR ! Layer depth of green roof @@ -316,7 +319,8 @@ SUBROUTINE urban(LSOLAR, & ! L U10,V10,TH2,Q2,UST,mh_urb,stdh_urb,lf_urb, & ! O lp_urb,hgt_urb,frc_urb,lb_urb,zo_check, & ! O CMCR,TGR,TGRL,SMR,CMGR_URB,CHGR_URB,jmonth, & ! H - DRELR,DRELB,DRELG,FLXHUMR,FLXHUMB,FLXHUMG) + DRELR,DRELB,DRELG,FLXHUMR,FLXHUMB,FLXHUMG, & + lf_urb_s, z0_urb, vegfrac_in) IMPLICIT NONE @@ -397,6 +401,13 @@ SUBROUTINE urban(LSOLAR, & ! L REAL, INTENT(INOUT), DIMENSION(4) :: lf_urb ! frontal area index [-] REAL, INTENT(INOUT) :: zo_check ! check for printing ZOC +!------------------------------------------------------------------------------- +! I: Distributed aerodynamics parameters +!------------------------------------------------------------------------------- + REAL, INTENT(IN) :: lf_urb_s ! frontal area index [-] + REAL, INTENT(IN) :: z0_urb ! roughness length [m] + REAL, INTENT(IN) :: vegfrac_in ! vegetation fraction (0 to 1) [-] + !------------------------------------------------------------------------------- ! O: output variables from Urban to LSM !------------------------------------------------------------------------------- @@ -544,10 +555,12 @@ SUBROUTINE urban(LSOLAR, & ! L REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10 REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST + REAL :: TSP, CHS_LOCAL, CHS2_LOCAL REAL :: WDR,HGT2,BW,DHGT REAL, parameter :: VonK = 0.4 REAL :: lambda_f,alpha_macd,beta_macd,lambda_fr + REAL :: lambda_p, vegfrac INTEGER :: iteration, K, NUDAPT INTEGER :: tloc, tloc2, Kalh @@ -595,6 +608,14 @@ SUBROUTINE urban(LSOLAR, & ! L integer,parameter :: IMPB = 2 integer,parameter :: IMPG = 3 + SHADOW = .false. +! SHADOW = .true. + + IF (distributed_aerodynamics_option .and. groption == 1) THEN + FATAL_ERROR("use_distributed_aerodynamics is not compatible with groption") + END IF + + !------------------------------------------------------------------------------- ! Set parameters !------------------------------------------------------------------------------- @@ -626,7 +647,7 @@ SUBROUTINE urban(LSOLAR, & ! L ! Glotfelty, 2012/07/05, NUDAPT Modification - if(mh_urb.gt.0.0)THEN + if (mh_urb.gt.0.0 .and. .not. distributed_aerodynamics_option) THEN !write(mesg,*) 'Mean Height NUDAPT',mh_urb !WRITE_MESSAGE(mesg) !write(mesg,*) 'Mean Height Table',ZR @@ -783,6 +804,25 @@ SUBROUTINE urban(LSOLAR, & ! L endif if(alhoption==1) ALH = ALH*alhdiuprf(tloc2)*alhseason(Kalh) + IF (distributed_aerodynamics_option) THEN + ZDC = 0. + IF (Z0_URB > MH_URB) THEN + FATAL_ERROR("Z0_URB is larger than MH_URB") + END IF + ZR = MAX(MH_URB, 3.5) + Z0C = MAX(Z0_URB, 0.1) + lambda_p = MAX(0.05, MIN(1.0, LP_URB)) + lambda_f = MAX(0.05, MIN(1.0, LF_URB_S)) + + R = lambda_p + RW = 1 - R + SVF = kanda_kawai_svf(lambda_p, lambda_f) + + vegfrac = MIN(0.9, MAX(0.1, vegfrac_in)) + + HGT = lambda_f + END IF + IF( ZDC+Z0C+2. >= ZA) THEN FATAL_ERROR("ZDC + Z0C + 2m is larger than the 1st WRF level - Stop in subroutine urban - change ZDC and Z0C" ) END IF @@ -818,6 +858,8 @@ SUBROUTINE urban(LSOLAR, & ! L TCP=TC QCP=QC + TSP = (TR * R + TB * W + TG * RW) / (R + RW + W) + !===Yang,2014/10/08, urban hydrological variables for single layer UCM=== FLXHUMRP = FLXHUMR FLXHUMBP = FLXHUMB @@ -865,9 +907,6 @@ SUBROUTINE urban(LSOLAR, & ! L ! Net Short Wave Radiation at roof, wall, and road !------------------------------------------------------------------------------- - SHADOW = .false. -! SHADOW = .true. - IF (SSG > 0.0) THEN IF(.NOT.SHADOW) THEN ! no shadow effects model @@ -968,9 +1007,22 @@ SUBROUTINE urban(LSOLAR, & ! L ! note that CHR_URB contains UA (=CHR_MOS*UA) RLMO_URB=0.0 + IF (distributed_aerodynamics_option) THEN + T1VC = TSP* (1.0+ 0.61 * QA) + CALL SFCDIF_URB (ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC,Z0HC,vegfrac) + CHC = CHC_URB / UA ! canopy bulk transfer coef. + ALPHAC = RHO * CP * CHC_URB + CHR = CHC * R / (R + W + RW) ! local bulk transfer coef for roof + CHB = CHC * W / (R + W + RW) ! local bulk transfer coef for building wall + CHG = CHC * RW / (R + W + RW) ! local bulk transfer coef for floor + ALPHAR = RHO * CP * CHR * UA + ALPHAB = RHO * CP * CHB * UA + ALPHAG = RHO * CP * CHG * UA + ELSE CALL SFCDIF_URB (ZA,Z0R,T1VR,TH2V,UA,AKANDA_URBAN,CMR_URB,CHR_URB,RLMO_URB,CDR) ALPHAR = RHO*CP*CHR_URB CHR=ALPHAR/RHO/CP/UA + END IF ! Yang, 03/12/2014 -- LH for impervious roof surface RAIN1 = RAIN * 0.001 /3600 ! CONVERT FROM mm/hr to m/s @@ -1190,6 +1242,8 @@ SUBROUTINE urban(LSOLAR, & ! L ! CALL mos(XXXC,ALPHAC,CDC,BHC,RIBC,Z,Z0C,UA,TA,TCP,RHO) ! Virtual temperatures needed by SFCDIF routine from Noah + IF (.not. distributed_aerodynamics_option) THEN + T1VC = TCP* (1.0+ 0.61 * QA) RLMO_URB=0.0 CALL SFCDIF_URB(ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC) @@ -1219,6 +1273,8 @@ SUBROUTINE urban(LSOLAR, & ! L CHB=ALPHAB/RHO/CP/UC CHG=ALPHAG/RHO/CP/UC + END IF + !Yang 10/10/2013 -- LH from impervious wall and ground IF (IMP_SCHEME==1) then BETB=0.0 @@ -1319,6 +1375,23 @@ SUBROUTINE urban(LSOLAR, & ! L DRGDTB=DRGDTB1+DRGDTB2 DRGDTG=DRGDTG1+DRGDTG2 + IF (distributed_aerodynamics_option) THEN + HB=RHO*CP*CHB*UA*(TBP-TA)*100. + HG=RHO*CP*CHG*UA*(TGP-TA)*100. + + DHBDTB=RHO*CP*CHB*UA*100. + DHBDTG=0. + DHGDTG=RHO*CP*CHG*UA*100. + DHGDTB=0. + + ELEB=RHO*EL*CHB*UA*BETB*(QS0B-QA)*100. + ELEG=RHO*EL*CHG*UA*BETG*(QS0G-QA)*100. + + DELEBDTB=RHO*EL*CHB*UA*BETB*DQS0BDTB*100. + DELEBDTG=0. + DELEGDTG=RHO*EL*CHG*UA*BETG*DQS0GDTG*100. + DELEGDTB=0. + ELSE HB=RHO*CP*CHB*UC*(TBP-TCP)*100. HG=RHO*CP*CHG*UC*(TGP-TCP)*100. @@ -1340,6 +1413,7 @@ SUBROUTINE urban(LSOLAR, & ! L DELEBDTG=RHO*EL*CHB*UC*BETB*(0.-DQCDTG)*100. DELEGDTG=RHO*EL*CHG*UC*BETG*(DQS0GDTG-DQCDTG)*100. DELEGDTB=RHO*EL*CHG*UC*BETG*(0.-DQCDTB)*100. + ENDIF G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.) G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.) @@ -1366,6 +1440,9 @@ SUBROUTINE urban(LSOLAR, & ! L TBP = TB TGP = TG + IF (distributed_aerodynamics_option) THEN + DTC = 0.0 + ELSE TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP TC=TC2/TC1 @@ -1377,6 +1454,7 @@ SUBROUTINE urban(LSOLAR, & ! L DTC=TCP - TC TCP=TC QCP=QC + END IF IF( ABS(F) < 0.000001 .AND. ABS(DTB) < 0.000001 & .AND. ABS(GF) < 0.000001 .AND. ABS(DTG) < 0.000001 & @@ -1484,7 +1562,11 @@ SUBROUTINE urban(LSOLAR, & ! L else FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG ) endif + IF (distributed_aerodynamics_option) THEN + FLXUV = CDC * UA * UA + ELSE FLXUV = ( R*CDR + RW*CDC )*UA*UA + END IF FLXG = ( R*G0R + W*G0B + RW*G0G ) LNET = R*RR + W*RB + RW*RG endif @@ -1532,14 +1614,20 @@ SUBROUTINE urban(LSOLAR, & ! L GZ1OZ0 = ALOG(Z/Z0) CD = 0.4**2./(ALOG(Z/Z0)-PSIM)**2. + CHS_LOCAL = 0.4 * UST / (ALOG(Z / Z0H) - PSIH) ! !m CH = 0.4**2./(ALOG(Z/Z0)-PSIM)/(ALOG(Z/Z0H)-PSIH) !m CHS = 0.4*UST/(ALOG(Z/Z0H)-PSIH) !m TS = TA + FLXTH/CH/UA ! surface potential temp (flux temp) !m QS = QA + FLXHUM/CH/UA ! surface humidity ! + IF (distributed_aerodynamics_option) THEN + TS = TA + FLXTH / (ALPHAC / (RHO * CP)) ! surface potential temp (flux temp) + QS = QA + FLXHUM / (ALPHAC / (RHO * CP)) ! surface humidity + ELSE TS = TA + FLXTH/CHS ! surface potential temp (flux temp) QS = QA + FLXHUM/CHS ! surface humidity + END IF !------------------------------------------------------- ! diagnostic GRID AVERAGED U10 V10 TH2 Q2 --> WRF @@ -1589,9 +1677,14 @@ SUBROUTINE urban(LSOLAR, & ! L ! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! potential temp at 2 m [K] ! TH2 = TS + (TA-TS)*(PSIT2/PSIT) ! Fei: this seems to be temp (not potential) at 2 m [K] !Fei: consistant with M-O theory + IF (distributed_aerodynamics_option) THEN + CHS2_LOCAL = 0.4 * UST / (ALOG(2. / Z0H) - PSIH2) + TH2 = TS + (TA - TS) * (CHS_LOCAL / CHS2_LOCAL) + Q2 = QS + (QA - QS) * (CHS_LOCAL / CHS2_LOCAL) + ELSE TH2 = TS + (TA-TS) *(CHS/CHS2) - Q2 = QS + (QA-QS)*(PSIT2/PSIT) ! humidity at 2 m [-] + END IF ! TS = (LW/SIG_SI/0.88)**0.25 ! Radiative temperature [K] @@ -1947,7 +2040,7 @@ END SUBROUTINE read_param ! !=============================================================================== SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & - sf_urban_physics,use_wudapt_lcz) + sf_urban_physics,use_wudapt_lcz, use_distributed_aerodynamics) ! num_roof_layers,num_wall_layers,num_road_layers) IMPLICIT NONE @@ -1962,6 +2055,7 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZG INTEGER, INTENT(IN) :: SF_URBAN_PHYSICS INTEGER, INTENT(IN) :: USE_WUDAPT_LCZ !AndreaLCZ + LOGICAL, INTENT(IN) :: use_distributed_aerodynamics INTEGER :: LC, K INTEGER :: IOSTATUS, ALLOCATE_STATUS @@ -1999,6 +2093,8 @@ SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, & ICATE=0 + distributed_aerodynamics_option = use_distributed_aerodynamics + if(USE_WUDAPT_LCZ.eq.0)then !AndreaLCZ OPEN (UNIT=11, & @@ -2647,6 +2743,14 @@ SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, !m !FS FRC_URB2D(I,J)=0. UTYPE_URB2D(I,J)=0 + + distributed_aerodynamics_check: IF (distributed_aerodynamics_option) THEN + IF (IVGTYP(I, J) == ISURBAN) THEN + UTYPE_URB2D(I, J) = 2 + ELSE + UTYPE_URB2D(I, J) = 0 + END IF + ELSE SWITCH_URB=1 IF( IVGTYP(I,J) == ISURBAN) THEN @@ -2729,6 +2833,7 @@ SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP, ENDDO ENDIF ENDIF + END IF distributed_aerodynamics_check QC_URB2D(I,J)=0.01 @@ -3007,7 +3112,7 @@ SUBROUTINE bisection(TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ,TS END SUBROUTINE bisection !=========================================================================== -SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD) +SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD,ZT_OUT,VEGFRAC) ! ---------------------------------------------------------------------- ! SUBROUTINE SFCDIF_URB (Urban version of SFCDIF_off) @@ -3026,6 +3131,8 @@ SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD) REAL SFCSPD, AKANDA, AKMS, AKHS, ZU, ZT, RDZ, CXCH REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4 + REAL, INTENT(OUT), OPTIONAL :: ZT_OUT + REAL, INTENT(IN), OPTIONAL :: VEGFRAC !CC ......REAL ZTFC REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, & @@ -3107,7 +3214,12 @@ SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD) ! ---------------------------------------------------------------------- ! KCL/TL Try Kanda approach instead (Kanda et al. 2007, JAMC) ! ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0 + IF (PRESENT(VEGFRAC)) THEN + ! Kawai et al. (2009) JAMC + ZT = EXP (2.0-(AKANDA-0.9*VEGFRAC**0.29)*(SQVISC**2 * USTAR * Z0)**0.25)* Z0 + ELSE ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0 + END IF ZSLU = ZLM + ZU @@ -3176,7 +3288,12 @@ SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD) USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST) !KCL/TL !ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0 + IF (PRESENT(VEGFRAC)) THEN + ! Kawai et al. (2009) JAMC + ZT = EXP (2.0-(AKANDA-0.9*VEGFRAC**0.29)*(SQVISC**2 * USTAR * Z0)**0.25)* Z0 + ELSE ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0 + END IF ZSLT = ZLM + ZT RLOGT = log (ZSLT / ZT) USTARK = USTAR * VKRM @@ -3200,6 +3317,8 @@ SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD) END DO CD = USTAR*USTAR/SFCSPD**2 + + IF (PRESENT(ZT_OUT)) ZT_OUT = ZT ! ---------------------------------------------------------------------- END SUBROUTINE SFCDIF_URB ! ---------------------------------------------------------------------- @@ -4055,5 +4174,17 @@ SUBROUTINE TDFCND (DF, SMC, QZ, SMCMAX) ! ---------------------------------------------------------------------- END SUBROUTINE TDFCND ! ---------------------------------------------------------------------- + + FUNCTION kanda_kawai_svf(lp, lf) RESULT (svf) + IMPLICIT NONE + real, intent(in) :: lp, lf + real :: hovl, vloc, vmod, svf + + hovl = lf * lp ** (-0.5) / (1. - lp ** 0.5) + vloc = cos(atan(2. * hovl)) * (2. - 4. / piconst * atan(cos(atan(2. * hovl)))) + vmod = 0.1120 * lp * vloc - 0.4817 * lp + 0.0246 * vloc + 0.9570 + svf = vloc * vmod + END FUNCTION kanda_kawai_svf + !=========================================================================== END MODULE module_sf_urban diff --git a/phys/module_surface_driver.F b/phys/module_surface_driver.F index f1592a1f00..557f2f606b 100644 --- a/phys/module_surface_driver.F +++ b/phys/module_surface_driver.F @@ -32,7 +32,7 @@ SUBROUTINE surface_driver( & & ,zs & & ,albsi, icedepth,snowsi & & ,xicem,isice,iswater,ct,tke_pbl & - & ,albbck,embck,lh,sh2o,shdmax,shdmin,z0 & + & ,albbck,embck,lh,sh2o,shdmax,shdmin,shdavg,z0 & & ,flqc,flhc,psfc,sst,sst_input,sstsk,dtw,sst_update,sst_skin & & ,scm_force_skintemp,scm_force_flux,t2,emiss & & ,sf_sfclay_physics,sf_surface_physics,ra_lw_physics & @@ -270,6 +270,7 @@ SUBROUTINE surface_driver( & & ,urban_map_zgrd & !multi-layer urban & ,num_urban_hi & !multi-layer urban & ,use_wudapt_lcz & !wudapt + & ,use_distributed_aerodynamics & !SLUCM & ,tsk_rural & !multi-layer urban & ,trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d & !multi-layer urban & ,tlev_urb3d,qlev_urb3d & !multi-layer urban @@ -288,6 +289,7 @@ SUBROUTINE surface_driver( & & ,swddir,swddif & !gl & ,lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d & !multi-layer urban & ,mh_urb2d,stdh_urb2d,lf_urb2d & + & ,lf_urb2d_s, z0_urb2d & & ,a_u_bep,a_v_bep,a_t_bep,a_q_bep & & ,b_u_bep,b_v_bep,b_t_bep,b_q_bep & & ,sf_bep,vl_bep & @@ -870,6 +872,7 @@ SUBROUTINE surface_driver( & REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), INTENT(INOUT):: SH2O REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDMAX REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDMIN + REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDAVG REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT ):: Z0 ! NoahMP specific fields @@ -953,6 +956,7 @@ SUBROUTINE surface_driver( & INTEGER , INTENT(IN) :: urban_map_zgrd INTEGER, INTENT(IN ):: NUM_URBAN_HI INTEGER, INTENT(IN ):: use_wudapt_lcz + LOGICAL, INTENT(IN ):: use_distributed_aerodynamics REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: tsk_rural REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d @@ -997,6 +1001,8 @@ SUBROUTINE surface_driver( & REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: mh_urb2d !urban REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: stdh_urb2d!urban REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN) :: lf_urb2d !urban + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lf_urb2d_s !urban + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: z0_urb2d !urban REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature @@ -1883,7 +1889,7 @@ SUBROUTINE surface_driver( & decided = .TRUE. END IF - IF ( run_param ) then + run_param_if: IF ( run_param ) then radiation = .false. frpcpn = .false. @@ -2664,7 +2670,7 @@ SUBROUTINE surface_driver( & myj,frpcpn, & sh2o,snowh, & !h u_phy,v_phy, & !I - snoalb,shdmin,shdmax, & !i + snoalb,shdmin,shdmax,shdavg, & !i snotime, & !o acsnom,acsnow, & !o snopcx, & !o @@ -2733,6 +2739,7 @@ SUBROUTINE surface_driver( & urban_map_zgrd, & !I multi-layer urban num_urban_hi, & !I multi-layer urban use_wudapt_lcz, & !I wudapt + use_distributed_aerodynamics, & !I SLUCM tsk_rural, & !H multi-layer urban trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, & !H multi-layer urban tlev_urb3d,qlev_urb3d, & !H multi-layer urban @@ -2748,6 +2755,7 @@ SUBROUTINE surface_driver( & lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d,& !GRZ lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, & !H multi-layer urban mh_urb2d,stdh_urb2D,lf_urb2d, & !SLUCM + lf_urb2d_s, z0_urb2d, & !SLUCM th_phy,rho,p_phy,ust, & !I multi-layer urban gmt,julday,xlong,xlat, & !I multi-layer urban a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban @@ -2880,6 +2888,7 @@ SUBROUTINE surface_driver( & lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d,& !GRZ lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, & !H multi-layer urban mh_urb2d,stdh_urb2D,lf_urb2d, & !SLUCM + lf_urb2d_s, z0_urb2d, & !SLUCM th_phy,rho,p_phy,ust, & !I multi-layer urban gmt,julday,xlong,xlat, & !I multi-layer urban a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban @@ -3217,6 +3226,7 @@ SUBROUTINE surface_driver( & dgr_urb3d, dg_urb3d, lfr_urb3d, lfg_urb3d, & !GRZ lp_urb2d, hi_urb2d, lb_urb2d, hgt_urb2d, & !H multi-layer urban mh_urb2d, stdh_urb2d, lf_urb2d, & !SLUCM + lf_urb2d_s, z0_urb2d, vegfra, & !SLUCM th_phy, rho, p_phy, ust, & !I multi-layer urban gmt, julday, xlong, xlat, & !I multi-layer urban a_u_bep, a_v_bep, a_t_bep, a_q_bep, & !O multi-layer urban @@ -4468,8 +4478,7 @@ SUBROUTINE surface_driver( & ENDIF ENDIF - - ENDIF + ENDIF run_param_if END SUBROUTINE surface_driver diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index 1acb3bda82..a3c2cb1ba8 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -504,6 +504,30 @@ END FUNCTION bep_bem_ngr_u END IF ENDDO +!----------------------------------------------------------------------- +! Check that only compatible options are set when use_distributed_aerodynamics is set +!----------------------------------------------------------------------- + IF (model_config_rec % use_distributed_aerodynamics) THEN + + IF (model_config_rec % use_wudapt_lcz) THEN + wrf_err_message = '--- ERROR: use_distributed_aerodynamics cannot work with use_wudapt_lcz' + CALL wrf_message ( wrf_err_message ) + count_fatal_error = count_fatal_error + 1 + END IF + + DO i = 1, model_config_rec % max_dom + IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE + IF ( (model_config_rec % sf_urban_physics(i) /= 0 .AND. & + model_config_rec % sf_urban_physics(i) /= 1 ) ) THEN + wrf_err_message = '--- ERROR: use_distributed_aerodynamics only works with urban options 1' + CALL wrf_message ( wrf_err_message ) + count_fatal_error = count_fatal_error + 1 + END IF + END DO + + END IF + + !----------------------------------------------------------------------- ! Check that channel irrigation is run with Noah !----------------------------------------------------------------------- diff --git a/share/output_wrf.F b/share/output_wrf.F index 1d07dcf97a..1e8aa825dd 100644 --- a/share/output_wrf.F +++ b/share/output_wrf.F @@ -668,6 +668,8 @@ SUBROUTINE output_wrf ( fid , grid , config_flags, switch , ierr ) CALL wrf_put_dom_ti_integer ( fid, 'SF_URBAN_PHYSICS', sf_urban_physics , 1 , ierr ) CALL wrf_put_dom_ti_integer ( fid, 'SF_SURFACE_MOSAIC', config_flags%sf_surface_mosaic , 1 , ierr ) CALL wrf_put_dom_ti_integer ( fid, 'SF_OCEAN_PHYSICS', config_flags%sf_ocean_physics , 1 , ierr ) + CALL wrf_put_dom_ti_logical ( fid, 'USE_DISTRIBUTED_AERODYNAMICS', config_flags%use_distributed_aerodynamics, 1, ierr) + CALL wrf_put_dom_ti_integer ( fid, 'DISTRIBUTED_AHE_OPT', config_flags%distributed_ahe_opt, 1, ierr) #endif IF ( switch .EQ. history_only ) THEN