Changeset 4713 for palm/trunk/SOURCE/urban_surface_mod.f90
- Timestamp:
- Sep 29, 2020 12:02:05 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/urban_surface_mod.f90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/urban_surface_mod.f90
r4701 r4713 27 27 ! ----------------- 28 28 ! $Id$ 29 ! - Do not change original fractions in USM energy balance 30 ! - Correct OpenMP parallelization 31 ! Author: J. Resler 32 ! 33 ! 4701 2020-09-27 11:02:15Z maronga 29 34 ! Corrected parameter 33 for building_type 2 (ground floor window emissivity 30 35 ! … … 6763 6768 ! 6764 6769 !-- First, treat horizontal surface elements 6765 !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_surface, lambda_surface_window,&6766 !$OMP& lambda_surface_ green, qv1, rho_cp, rho_lv, drho_l_lv, f_shf,&6767 !$OMP& f_shf_window, f_shf_green, m_total, f1, f2, e_s, e, f3, f_qsws_veg,&6768 !$OMP& q_s, f_qsws_liq, f_qsws, e_s_dt, dq_s_dt, coef_1, coef_window_1,&6769 !$OMP& coef_ green_1, coef_2, coef_window_2, coef_green_2, stend_wall,&6770 !$OMP& stend_w indow, stend_green, tend, m_liq_max)6770 !$OMP PARALLEL PRIVATE (m, i, j, k, frac_win, frac_wall, frac_green, lambda_surface, ueff, & 6771 !$OMP& lambda_surface_window, lambda_surface_green, qv1, rho_cp, rho_lv, & 6772 !$OMP& drho_l_lv, f_shf, f_shf_window, f_shf_green, m_total, f1, f2, e_s, e, & 6773 !$OMP& f3, f_qsws_veg, q_s, f_qsws_liq, f_qsws, e_s_dt, dq_s_dt, coef_1, & 6774 !$OMP& coef_window_1, coef_green_1, coef_2, coef_window_2, coef_green_2, & 6775 !$OMP& stend_wall, stend_window, stend_green, tend, m_liq_max) 6771 6776 !$OMP DO SCHEDULE (STATIC) 6772 6777 DO m = 1, surf%ns … … 6775 6780 !-- Note, this is a temporary fix and needs to be removed later. 6776 6781 IF ( during_spinup ) THEN 6777 frac_win = surf%frac(m,ind_wat_win) 6778 frac_wall = surf%frac(m,ind_veg_wall) 6779 frac_green = surf%frac(m,ind_pav_green) 6780 surf%frac(m,ind_wat_win) = 0.0_wp 6781 surf%frac(m,ind_veg_wall) = 1.0_wp 6782 surf%frac(m,ind_pav_green) = 0.0_wp 6782 frac_win = 0.0_wp 6783 frac_wall = 1.0_wp 6784 frac_green = 0.0_wp 6785 ELSE 6786 frac_win = surf%frac(m,ind_wat_win) 6787 frac_wall = surf%frac(m,ind_veg_wall) 6788 frac_green = surf%frac(m,ind_pav_green) 6783 6789 ENDIF 6784 6790 ! … … 6803 6809 rho_cp = c_p * hyp(k) / ( r_d * surf%pt1(m) * exner(k) ) 6804 6810 6805 IF ( surf%frac(m,ind_pav_green)> 0.0_wp ) THEN6811 IF ( frac_green > 0.0_wp ) THEN 6806 6812 ! 6807 6813 !-- Calculate frequently used parameters … … 6843 6849 !-- to calculation of roughness relative to concrete. Note, wind velocity is limited 6844 6850 !-- to avoid division by zero. The nominator can become <= 0.0 for values z0 < 3*10E-4. 6845 ueff = MAX ( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + 6846 ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + 6847 ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ), 6848 ( ( 4.0_wp + 0.1_wp ) 6849 / ( surf%z0(m) * d_roughness_concrete ) &6850 - 11.8_wp ) / 4.2_wp 6851 ueff = MAX ( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + & 6852 ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + & 6853 ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ), & 6854 ( ( 4.0_wp + 0.1_wp ) & 6855 / ( surf%z0(m) * d_roughness_concrete ) & 6856 - 11.8_wp ) / 4.2_wp & 6851 6857 ) 6852 6858 6853 surf%r_a(m) = rho_cp / ( surf%z0(m) * d_roughness_concrete &6859 surf%r_a(m) = rho_cp / ( surf%z0(m) * d_roughness_concrete & 6854 6860 * ( 11.8_wp + 4.2_wp * ueff ) - 4.0_wp ) 6855 6861 ENDIF … … 6869 6875 f_shf_green = rho_cp / surf%r_a_green(m) 6870 6876 6871 IF ( surf%frac(m,ind_pav_green)> 0.0_wp ) THEN6877 IF ( frac_green > 0.0_wp ) THEN 6872 6878 ! 6873 6879 !-- Adapted from LSM: … … 6876 6882 6877 6883 !-- f1: Correction for incoming shortwave radiation (stomata close at night) 6878 f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) / &6884 f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) / & 6879 6885 (0.81_wp * ( 0.004_wp * surf%rad_sw_in(m) + 1.0_wp ) ) ) 6880 6886 ! … … 6914 6920 !-- obsolete, as r_canopy is not used below. 6915 6921 !-- To do: check for very dry soil -> r_canopy goes to infinity 6916 surf%r_canopy(m) = surf%r_canopy_min(m) / &6922 surf%r_canopy(m) = surf%r_canopy_min(m) / & 6917 6923 ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp ) 6918 6924 ! … … 6968 6974 + f_shf * surf%pt1(m) + lambda_surface * t_wall%val(nzb_wall,m) 6969 6975 6970 IF ( ( .NOT. during_spinup ) .AND. ( surf%frac(m,ind_wat_win)> 0.0_wp ) ) THEN6976 IF ( ( .NOT. during_spinup ) .AND. (frac_win > 0.0_wp ) ) THEN 6971 6977 coef_window_1 = surf%rad_net_l(m) + ( 3.0_wp + 1.0_wp ) & 6972 6978 * surf%emissivity(m,ind_wat_win) * sigma_sb & … … 6974 6980 + lambda_surface_window * t_window%val(nzb_wall,m) 6975 6981 ENDIF 6976 IF ( ( humidity ) .AND. ( surf%frac(m,ind_pav_green)> 0.0_wp ) ) THEN6982 IF ( ( humidity ) .AND. ( frac_green > 0.0_wp ) ) THEN 6977 6983 coef_green_1 = surf%rad_net_l(m) + ( 3.0_wp + 1.0_wp ) & 6978 6984 * surf%emissivity(m,ind_pav_green) * sigma_sb & … … 6990 6996 coef_2 = 4.0_wp * surf%emissivity(m,ind_veg_wall) * sigma_sb * t_surf_wall%val(m)**3 & 6991 6997 + lambda_surface + f_shf / exner(k) 6992 IF ( ( .NOT. during_spinup ) .AND. ( surf%frac(m,ind_wat_win)> 0.0_wp ) ) THEN6998 IF ( ( .NOT. during_spinup ) .AND. ( frac_win > 0.0_wp ) ) THEN 6993 6999 coef_window_2 = 4.0_wp * surf%emissivity(m,ind_wat_win) * sigma_sb * & 6994 7000 t_surf_window%val(m)**3 + lambda_surface_window + f_shf_window / exner(k) 6995 7001 ENDIF 6996 IF ( ( humidity ) .AND. ( surf%frac(m,ind_pav_green)> 0.0_wp ) ) THEN7002 IF ( ( humidity ) .AND. ( frac_green > 0.0_wp ) ) THEN 6997 7003 coef_green_2 = 4.0_wp * surf%emissivity(m,ind_pav_green) * sigma_sb * & 6998 7004 t_surf_green%val(m)**3 + f_qsws * dq_s_dt + lambda_surface_green & … … 7007 7013 * t_surf_wall%val(m) ) & 7008 7014 / ( surf%c_surface(m) + coef_2 * dt_3d * tsc(2) ) 7009 IF ( ( .NOT. during_spinup ) .AND. ( surf%frac(m,ind_wat_win)> 0.0_wp) ) THEN7015 IF ( ( .NOT. during_spinup ) .AND. (frac_win > 0.0_wp) ) THEN 7010 7016 t_surf_window_p%val(m) = ( coef_window_1 * dt_3d * tsc(2) + & 7011 7017 surf%c_surface_window(m) * t_surf_window%val(m) ) / & … … 7027 7033 !-- vpt_surface, which is, due to the lack of moisture on roofs, simply assumed to be the surface 7028 7034 !-- temperature. 7029 surf%pt_surface(m) = ( surf%frac(m,ind_veg_wall) * t_surf_wall_p%val(m)&7030 + surf%frac(m,ind_wat_win) * t_surf_window_p%val(m)&7031 + surf%frac(m,ind_pav_green) * t_surf_green_p%val(m)&7035 surf%pt_surface(m) = ( frac_wall * t_surf_wall_p%val(m) & 7036 + frac_win * t_surf_window_p%val(m) & 7037 + frac_green * t_surf_green_p%val(m) & 7032 7038 ) / exner(k) 7033 7039 … … 7069 7075 !-- Calculate new fluxes 7070 7076 !-- Rad_net_l is never used! 7071 surf%rad_net_l(m) = surf%rad_net_l(m) + surf%frac(m,ind_veg_wall)&7077 surf%rad_net_l(m) = surf%rad_net_l(m) + frac_wall & 7072 7078 * sigma_sb * surf%emissivity(m,ind_veg_wall) & 7073 7079 * ( t_surf_wall_p%val(m)**4 - t_surf_wall%val(m)**4 ) & 7074 + surf%frac(m,ind_wat_win) * sigma_sb&7080 + frac_win * sigma_sb & 7075 7081 * surf%emissivity(m,ind_wat_win) & 7076 7082 * ( t_surf_window_p%val(m)**4 - t_surf_window%val(m)**4 ) & 7077 + surf%frac(m,ind_pav_green) * sigma_sb&7083 + frac_green * sigma_sb & 7078 7084 * surf%emissivity(m,ind_pav_green) & 7079 7085 * ( t_surf_green_p%val(m)**4 - t_surf_green%val(m)**4 ) … … 7088 7094 !-- Ground/wall/roof surface heat flux 7089 7095 surf%wshf_eb(m) = - f_shf * ( surf%pt1(m) - t_surf_wall_p%val(m) / exner(k) ) & 7090 * surf%frac(m,ind_veg_wall) - f_shf_window&7096 * frac_wall - f_shf_window & 7091 7097 * ( surf%pt1(m) - t_surf_window_p%val(m) / exner(k) ) & 7092 * surf%frac(m,ind_wat_win) - f_shf_green&7098 * frac_win - f_shf_green & 7093 7099 * ( surf%pt1(m) - t_surf_green_p%val(m) / exner(k) ) & 7094 * surf%frac(m,ind_pav_green)7100 * frac_green 7095 7101 ! 7096 7102 !-- Store kinematic surface heat fluxes for utilization in other processes diffusion_s, … … 7103 7109 ENDIF 7104 7110 7105 IF ( humidity .AND. surf%frac(m,ind_pav_green)> 0.0_wp ) THEN!7111 IF ( humidity .AND. frac_green > 0.0_wp ) THEN! 7106 7112 !-- Calculate true surface resistance 7107 7113 IF ( upward ) THEN … … 7125 7131 IF ( m_liq_usm_h(l)%val(m) /= m_liq_max ) THEN 7126 7132 surf%qsws_liq(m) = surf%qsws_liq(m) & 7127 + surf%frac(m,ind_pav_green)&7133 + frac_green & 7128 7134 * prr(k+k_off,j+j_off,i+i_off) * hyrho(k+k_off) & 7129 7135 * 0.001_wp * rho_l * l_v … … 7191 7197 ELSE 7192 7198 surf%r_s(m) = 1.0E10_wp 7193 ENDIF7194 !7195 !-- During spinup green and window fraction are set to zero. Here, the original values are7196 !-- restored.7197 IF ( during_spinup ) THEN7198 surf%frac(m,ind_wat_win) = frac_win7199 surf%frac(m,ind_veg_wall) = frac_wall7200 surf%frac(m,ind_pav_green) = frac_green7201 7199 ENDIF 7202 7200
Note: See TracChangeset
for help on using the changeset viewer.
