Changeset 4713 for palm/trunk
- Timestamp:
- Sep 29, 2020 12:02:05 PM (4 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
land_surface_model_mod.f90 (modified) (8 diffs)
-
radiation_model_mod.f90 (modified) (30 diffs)
-
urban_surface_mod.f90 (modified) (18 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/land_surface_model_mod.f90
r4694 r4713 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Optimize calculation in case of calc_soil_moisture=false 28 ! - Improve OMP parallelization 29 ! - Do not calculate surface prognostic temperature for water surfaces (avoid numerical 30 ! instabilities) 31 ! Author: J. Resler 32 ! 33 ! 4694 2020-09-23 15:09:19Z pavelkrc 27 34 ! Fix reading of surface data from MPI restart file 28 35 ! … … 1734 1741 i_off = surf%ioff 1735 1742 1736 !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface, 1743 !$OMP PARALLEL PRIVATE (m, i, j, k, lambda_h_sat, ke, lambda_soil, lambda_surface, ueff, & 1737 1744 !$OMP& c_surface_tmp, f1,m_total, f2, e_s, e, f3, m_min, m_liq_max, q_s, & 1738 1745 !$OMP& f_qsws_veg, f_qsws_soil, f_qsws_liq, f_shf, f_qsws, e_s_dt, dq_s_dt, & … … 2038 2045 !-- Implicit solution when the surface layer has no heat capacity, 2039 2046 !-- otherwise use RK3 scheme. 2040 surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *& 2041 surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2& 2047 IF ( surf%water_surface(m) ) THEN 2048 ! 2049 !-- water temperature is fixed (avoid numerical inaccuracies) 2050 surf_t_surface_p%var_1d(m) = surf_t_surface%var_1d(m) 2051 ELSE 2052 surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & 2053 surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2 & 2042 2054 * dt_3d * tsc(2) ) 2043 2055 ENDIF 2044 2056 ! 2045 2057 !-- Add RK3 term … … 5545 5557 ENDIF 5546 5558 5547 !$OMP PARALLEL PRIVATE (m, k, lambda_temp, lambda_h_sat, ke, tend, gamma_temp, h_vg, m_total) 5559 !$OMP PARALLEL PRIVATE (m, k, lambda_temp, lambda_h_sat, ke, tend, gamma_temp, h_vg, & 5560 !$OMP& m_total, root_extr) 5548 5561 !$OMP DO SCHEDULE (STATIC) 5549 5562 DO m = 1, surf%ns … … 5636 5649 5637 5650 5638 DO k = nzb_soil, nzt_soil5639 5640 !5641 !-- In order to prevent water tranport through paved surfaces,5642 !-- conductivity and diffusivity are set to zero5643 IF ( surf%pavement_surface(m) .AND. &5644 k <= surf%nzt_pavement(m) ) THEN5645 lambda_temp(k) = 0.0_wp5646 gamma_temp(k) = 0.0_wp5647 5648 ELSE5649 5650 !5651 !-- Calculate soil diffusivity at the center of the soil layers5652 lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat &5653 / surf%m_sat(k,m) ) * ( &5654 MAX( surf_m_soil%var_2d(k,m), &5655 surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**( &5656 b_ch + 2.0_wp )5657 5658 !5659 !-- Parametrization of Van Genuchten5660 !-- Calculate the hydraulic conductivity after Van Genuchten (1980)5661 h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) / &5662 ( surf%m_res(k,m) - &5663 MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )&5664 ) &5665 )**( &5666 surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) &5667 ) - 1.0_wp &5668 )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m)5669 5670 gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp + &5671 ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m) &5672 )**( &5673 1.0_wp - 1.0_wp / surf%n_vg(k,m)) - ( &5674 surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m) &5675 - 1.0_wp) )**2 ) &5676 / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg &5677 )**surf%n_vg(k,m) )**( ( 1.0_wp - 1.0_wp &5678 / surf%n_vg(k,m) ) * &5679 ( surf%l_vg(k,m) + 2.0_wp) ) )5680 5681 ENDIF5682 5683 ENDDO5684 5685 5686 5651 IF ( calc_soil_moisture ) THEN 5687 5652 5653 DO k = nzb_soil, nzt_soil 5654 ! 5655 !-- In order to prevent water tranport through paved surfaces, 5656 !-- conductivity and diffusivity are set to zero 5657 IF ( surf%pavement_surface(m) .AND. & 5658 k <= surf%nzt_pavement(m) ) THEN 5659 lambda_temp(k) = 0.0_wp 5660 gamma_temp(k) = 0.0_wp 5661 5662 ELSE 5663 5664 ! 5665 !-- Calculate soil diffusivity at the center of the soil layers 5666 lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat & 5667 / surf%m_sat(k,m) ) * ( & 5668 MAX( surf_m_soil%var_2d(k,m), & 5669 surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**( & 5670 b_ch + 2.0_wp ) 5671 5672 ! 5673 !-- Parametrization of Van Genuchten 5674 !-- Calculate the hydraulic conductivity after Van Genuchten (1980) 5675 h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) / & 5676 ( surf%m_res(k,m) - & 5677 MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )& 5678 ) & 5679 )**( & 5680 surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) & 5681 ) - 1.0_wp & 5682 )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m) 5683 5684 gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp + & 5685 ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m) & 5686 )**( & 5687 1.0_wp - 1.0_wp / surf%n_vg(k,m)) - ( & 5688 surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m) & 5689 - 1.0_wp) )**2 ) & 5690 / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg & 5691 )**surf%n_vg(k,m) )**( ( 1.0_wp - 1.0_wp & 5692 / surf%n_vg(k,m) ) * & 5693 ( surf%l_vg(k,m) + 2.0_wp) ) ) 5694 5695 ENDIF 5696 5697 ENDDO 5688 5698 ! 5689 5699 !-- Prognostic equation for soil moisture content. Only performed, … … 7157 7167 ALLOCATE( t_soil_h(1)%var_2d(nzb_soil:nzt_soil+1, & 7158 7168 1:surf_lsm_h(1)%ns) ) 7159 write(9,*) surf_lsm_h(1)%ns7160 7169 READ ( 13 ) tmp_walltype_h_2d(1)%var_2d 7161 7170 ENDIF … … 7753 7762 ! REAL(wp) :: re_0 !< near-surface roughness Reynolds number 7754 7763 7764 !$OMP DO PRIVATE (m, i, j) SCHEDULE (STATIC) REDUCTION(.OR.:flag_exceed_z0, flag_exceed_z0h) 7755 7765 DO m = 1, surf_lsm_h(0)%ns 7756 7766 !-- only upward facin horizontal surfaces are considered for water surface processing 7757 i = surf_lsm_h(0)%i(m)7758 j = surf_lsm_h(0)%j(m)7759 7767 7760 7768 IF ( surf_lsm_h(0)%water_surface(m) ) THEN 7761 7769 i = surf_lsm_h(0)%i(m) 7770 j = surf_lsm_h(0)%j(m) 7762 7771 ! 7763 7772 !-- Disabled: FLake parameterization. Ideally, the Charnock … … 7813 7822 ENDIF 7814 7823 ENDDO 7824 !$OMP END DO 7815 7825 #if defined( __parallel ) 7816 7826 CALL MPI_ALLREDUCE( MPI_IN_PLACE, flag_exceed_z0, 1, MPI_LOGICAL, & -
palm/trunk/SOURCE/radiation_model_mod.f90
r4708 r4713 28 28 ! ----------------- 29 29 ! $Id$ 30 ! Correct OpenMP parallelization including cycles with cumulative variables (J. Resler) 31 ! 32 ! 4708 2020-09-28 17:42:58Z suehring 30 33 ! - Bugfix, correct mapping of RRTMG heating rates, as well as incoming/outgoing 31 34 ! radiation onto the topography-following grid … … 5967 5970 REAL(wp) :: pabs_pc_lwdifl !< total absorbed LW radiation in plant canopy from sky in local processor (W) 5968 5971 REAL(wp) :: pabs_pc_lwdif !< total absorbed LW radiation in plant canopy from sky in all processors (W) 5969 5972 !- rotation related variables 5970 5973 REAL(wp) :: sun_direct_factor !< factor for direct normal radiation from direct horizontal 5971 5974 REAL(wp) :: sin_rot !< sine of rotation_angle … … 6072 6075 !-- Set up thermal radiation from surfaces 6073 6076 mm = 1 6077 !-- following code depends on the order of the execution 6078 !-- do not parallelize by OpenMP 6074 6079 DO i = nxl, nxr 6075 6080 DO j = nys, nyn … … 6159 6164 6160 6165 IF ( surface_reflections) THEN 6166 !$OMP DO PRIVATE (i, j, k, isvf, isurf, isurfsrc) SCHEDULE (STATIC) 6161 6167 DO isvf = 1, nsvfl 6162 6168 isurf = svfsurf(1, isvf) … … 6173 6179 ENDIF 6174 6180 ENDDO 6181 !$OMP END DO 6175 6182 ENDIF 6176 6183 ! 6177 6184 !-- diffuse radiation using sky view factor 6185 !$OMP DO PRIVATE (i, j, d, isurf) REDUCTION(+:pinswl, pinlwl) SCHEDULE (STATIC) 6178 6186 DO isurf = 1, nsurfl 6179 6187 j = surfl(iy, isurf) … … 6193 6201 ENDIF 6194 6202 ENDDO 6203 !$OMP END DO 6195 6204 ! 6196 6205 !-- MRT diffuse irradiance 6206 !$OMP DO PRIVATE (i, j, imrt) SCHEDULE (STATIC) 6197 6207 DO imrt = 1, nmrtbl 6198 6208 j = mrtbl(iy, imrt) … … 6201 6211 mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i) 6202 6212 ENDDO 6213 !$OMP END DO 6203 6214 ! 6204 6215 !-- Direct radiation … … 6218 6229 isd = dsidir_rev(j, i) 6219 6230 !-- TODO: check if isd = -1 to report that this solar position is not precalculated 6231 !$OMP DO PRIVATE (i, j, d, isurf) REDUCTION(+:pinswl) SCHEDULE (STATIC) 6220 6232 DO isurf = 1, nsurfl 6221 6233 j = surfl(iy, isurf) … … 6227 6239 pinswl = pinswl + surfinswdir(isurf) * facearea(d) 6228 6240 ENDDO 6241 !$OMP END DO 6229 6242 ! 6230 6243 !-- MRT direct irradiance 6244 !$OMP DO PRIVATE (i, j, imrt) SCHEDULE (STATIC) 6231 6245 DO imrt = 1, nmrtbl 6232 6246 j = mrtbl(iy, imrt) … … 6235 6249 * sun_direct_factor / 4.0_wp ! normal to sphere 6236 6250 ENDDO 6251 !$OMP END DO 6237 6252 ENDIF 6238 6253 ! 6239 6254 !-- MRT first pass thermal 6255 !$OMP DO PRIVATE (imrtf, imrt, isurfsrc) SCHEDULE (STATIC) 6240 6256 DO imrtf = 1, nmrtf 6241 6257 imrt = mrtfsurf(1, imrtf) … … 6243 6259 mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc) 6244 6260 ENDDO 6261 !$OMP END DO 6245 6262 ! 6246 6263 !-- Absorption in each local plant canopy grid box from the first atmospheric … … 6252 6269 pcbinlw(:) = 0.0_wp 6253 6270 6271 !$OMP DO PRIVATE (icsf, ipcgb, i, j, k, kk, isurfsrc, pc_abs_frac, pcrad, asrc) & 6272 !$OMP& REDUCTION(+:pinswl, pinlwl, pabslwl, pemitlwl, pabs_pc_lwdifl) SCHEDULE (STATIC) 6254 6273 DO icsf = 1, ncsfl 6255 6274 ipcgb = csfsurf(1, icsf) … … 6309 6328 ENDIF 6310 6329 ENDDO 6330 !$OMP END DO 6311 6331 6312 6332 pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:) … … 6362 6382 ENDIF 6363 6383 6364 !-- Next passes of radiation interactions:6384 !-- Next passes of radiation interactions: 6365 6385 !-- radiation reflections 6366 6386 … … 6404 6424 ! 6405 6425 !-- Reflected radiation 6426 !$OMP DO PRIVATE (isvf, isurf, isurfsrc) SCHEDULE (STATIC) 6406 6427 DO isvf = 1, nsvfl 6407 6428 isurf = svfsurf(1, isvf) … … 6414 6435 ENDIF 6415 6436 ENDDO 6437 !$OMP END DO 6416 6438 ! 6417 6439 !-- NOTE: PC absorbtion and MRT from reflected can both be done at once … … 6421 6443 ! 6422 6444 !-- Radiation absorbed by plant canopy 6445 !$OMP DO PRIVATE (icsf, ipcgb, isurfsrc, asrc) SCHEDULE (STATIC) 6423 6446 DO icsf = 1, ncsfl 6424 6447 ipcgb = csfsurf(1, icsf) … … 6435 6458 ENDIF 6436 6459 ENDDO 6460 !$OMP END DO 6437 6461 ! 6438 6462 !-- MRT reflected 6463 !$OMP DO PRIVATE (imrtf, imrt, isurfsrc) SCHEDULE (STATIC) 6439 6464 DO imrtf = 1, nmrtf 6440 6465 imrt = mrtfsurf(1, imrtf) … … 6443 6468 mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc) 6444 6469 ENDDO 6470 !$OMP END DO 6445 6471 6446 6472 IF ( trace_fluxes_above >= 0.0_wp ) THEN … … 6464 6490 IF ( npcbl > 0 ) THEN 6465 6491 pcm_heating_rate(:,:,:) = 0.0_wp 6492 !$OMP DO PRIVATE (ipcgb, i, j, k, kk) REDUCTION(+:pabsswl) SCHEDULE (STATIC) 6466 6493 DO ipcgb = 1, npcbl 6467 6494 j = pcbl(iy, ipcgb) … … 6476 6503 pabsswl = pabsswl + pcbinsw(ipcgb) 6477 6504 ENDDO 6505 !$OMP END DO 6478 6506 6479 6507 IF ( humidity .AND. plant_canopy_transpiration ) THEN … … 6481 6509 pcm_transpiration_rate(:,:,:) = 0.0_wp 6482 6510 pcm_latent_rate(:,:,:) = 0.0_wp 6511 !$OMP DO PRIVATE (ipcgb, i, j, k, kk) SCHEDULE (STATIC) 6483 6512 DO ipcgb = 1, npcbl 6484 6513 i = pcbl(ix, ipcgb) … … 6490 6519 pcm_transpiration_rate(kk,j,i), & 6491 6520 pcm_latent_rate(kk,j,i) ) 6492 ENDDO 6521 ENDDO 6522 !$OMP END DO 6493 6523 ENDIF 6494 6524 ENDIF … … 6506 6536 ! and claculate relevant radiation model-RTM coupling terms 6507 6537 mm = 1 6538 !-- following code depends on the order of the execution 6539 !-- do not parallelize by OpenMP 6508 6540 DO i = nxl, nxr 6509 6541 DO j = nys, nyn … … 6598 6630 ENDDO 6599 6631 6632 !$OMP DO PRIVATE (i, d) REDUCTION(+:pabsswl, pabslwl, pemitlwl, pabs_surf_lwdifl) & 6633 !$OMP& SCHEDULE (STATIC) 6600 6634 DO i = 1, nsurfl 6601 6635 d = surfl(id, i) … … 6612 6646 emiss_surf(i) * facearea(d) * surfinlwdif(i) 6613 6647 ENDDO 6648 !$OMP END DO 6614 6649 6615 6650 DO l = 0, 1 6651 !$OMP DO PRIVATE (m) SCHEDULE (STATIC) 6616 6652 DO m = 1, surf_usm_h(l)%ns 6617 6653 surf_usm_h(l)%surfhf(m) = surf_usm_h(l)%rad_sw_in(m) + & … … 6620 6656 surf_usm_h(l)%rad_lw_out(m) 6621 6657 ENDDO 6658 !$OMP END DO 6659 !$OMP DO PRIVATE (m) SCHEDULE (STATIC) 6622 6660 DO m = 1, surf_lsm_h(l)%ns 6623 6661 surf_lsm_h(l)%surfhf(m) = surf_lsm_h(l)%rad_sw_in(m) + & … … 6626 6664 surf_lsm_h(l)%rad_lw_out(m) 6627 6665 ENDDO 6666 !$OMP END DO 6628 6667 ENDDO 6668 6629 6669 DO l = 0, 3 6670 !$OMP DO PRIVATE (m) SCHEDULE (STATIC) 6630 6671 !-- urban 6631 6672 DO m = 1, surf_usm_v(l)%ns … … 6635 6676 surf_usm_v(l)%rad_lw_out(m) 6636 6677 ENDDO 6678 !$OMP END DO 6637 6679 !-- land 6680 !$OMP DO PRIVATE (m) SCHEDULE (STATIC) 6638 6681 DO m = 1, surf_lsm_v(l)%ns 6639 6682 surf_lsm_v(l)%surfhf(m) = surf_lsm_v(l)%rad_sw_in(m) + & … … 6643 6686 6644 6687 ENDDO 6688 !$OMP END DO 6645 6689 ENDDO 6646 6690 ! -
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.
