Changeset 1691 for palm/trunk/SOURCE/land_surface_model.f90
- Timestamp:
- Oct 26, 2015 4:17:44 PM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/land_surface_model.f90 (modified) (54 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/land_surface_model.f90
r1683 r1691 14 14 ! PALM. If not, see <http://www.gnu.org/licenses/>. 15 15 ! 16 ! Copyright 1997-201 4Leibniz Universitaet Hannover16 ! Copyright 1997-2015 Leibniz Universitaet Hannover 17 17 !--------------------------------------------------------------------------------! 18 18 ! 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! 21 ! Added skip_time_do_lsm to allow for spin-ups without LSM. Various bugfixes: 22 ! Soil temperatures are now defined at the edges of the layers, calculation of 23 ! shb_eb corrected, prognostic equation for skin temperature corrected. Surface 24 ! fluxes are now directly transfered to atmosphere 22 25 ! 23 26 ! Former revisions: … … 73 76 !> DALES and UCLA-LES models. 74 77 !> 75 !> @todo Check dewfall parametrization for fog simulations. 76 !> @todo Consider partial absorption of the net shortwave radiation by the surface layer. 77 !> @todo Allow for water surfaces, check performance for bare soils. 78 !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, nzt_soil=3)). 79 !> @todo Implement surface runoff model (required when performing long-term LES with 80 !> considerable precipitation. 78 !> @todo Consider partial absorption of the net shortwave radiation by the 79 !> surface layer. 80 !> @todo Allow for water surfaces 81 !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 82 !> nzt_soil=3)). 83 !> @todo Implement surface runoff model (required when performing long-term LES 84 !> with considerable precipitation. 81 85 !> 82 !> @note No time step criterion is required as long as the soil layers do not become83 !> too thin.86 !> @note No time step criterion is required as long as the soil layers do not 87 !> become too thin. 84 88 !------------------------------------------------------------------------------! 85 89 MODULE land_surface_model_mod 86 90 87 USE arrays_3d,&88 ONLY: pt, pt_p, q_p, qsws, rif, shf, ts, us, z0, z0h89 90 USE cloud_parameters,&91 ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v92 93 USE control_parameters,&94 ONLY: cloud_physics, dt_3d, humidity, intermediate_timestep_count,&95 initializing_actions, intermediate_timestep_count_max,&96 max_masks, precipitation, pt_surface, rho_surface,&97 roughness_length, surface_pressure, timestep_scheme, tsc,&98 99 100 USE indices,&101 102 103 91 USE arrays_3d, & 92 ONLY: hyp, pt, pt_p, q, q_p, ql, qsws, rif, shf, ts, us, vpt, z0, z0h 93 94 USE cloud_parameters, & 95 ONLY: cp, hyrho, l_d_cp, l_d_r, l_v, prr, pt_d_t, rho_l, r_d, r_v 96 97 USE control_parameters, & 98 ONLY: cloud_physics, dt_3d, humidity, intermediate_timestep_count, & 99 initializing_actions, intermediate_timestep_count_max, & 100 max_masks, precipitation, pt_surface, rho_surface, & 101 roughness_length, surface_pressure, timestep_scheme, tsc, & 102 z0h_factor, time_since_reference_point 103 104 USE indices, & 105 ONLY: nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner 106 107 USE kinds 104 108 105 109 USE netcdf_control, & 106 110 ONLY: dots_label, dots_num, dots_unit 107 111 108 USE radiation_model_mod, & 109 ONLY: radiation_scheme, rad_net, rad_sw_in, sigma_sb 110 111 USE statistics, & 112 ONLY: hom, statistic_regions 112 USE pegrid 113 114 USE radiation_model_mod, & 115 ONLY: force_radiation_call, radiation_scheme, rad_net, rad_sw_in, & 116 rad_lw_out, sigma_sb 117 118 #if defined ( __rrtmg ) 119 USE radiation_model_mod, & 120 ONLY: rrtm_lwuflx_dt, rrtm_idrv 121 #endif 122 123 USE statistics, & 124 ONLY: hom, statistic_regions 113 125 114 126 IMPLICIT NONE … … 144 156 soil_type = 3 !< soil type, 0: user-defined, 1-6: generic (see list) 145 157 146 LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model 147 dewfall = .TRUE., & !< allow/inhibit dewfall 148 land_surface = .FALSE. !< flag parameter indicating wheather the lsm is used 158 LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model 159 dewfall = .TRUE., & !< allow/inhibit dewfall 160 force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls 161 land_surface = .FALSE. !< flag parameter indicating wheather the lsm is used 149 162 150 163 ! value 9999999.9_wp -> generic available or user-defined value must be set … … 160 173 hydraulic_conductivity = 9999999.9_wp, & !< NAMELIST gamma_w_sat 161 174 ke = 0.0_wp, & !< Kersten number 175 lambda_h_sat = 0.0_wp, & !< heat conductivity for saturated soil 162 176 lambda_surface_stable = 9999999.9_wp, & !< NAMELIST lambda_surface_s 163 177 lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u … … 174 188 rd_d_rv, & !< r_d / r_v 175 189 saturation_moisture = 9999999.9_wp, & !< NAMELIST m_sat 190 skip_time_do_lsm = 0.0_wp, & !< LSM is not called before this time 176 191 vegetation_coverage = 9999999.9_wp, & !< NAMELIST c_veg 177 192 wilting_point = 9999999.9_wp, & !< NAMELIST m_wilt 178 193 z0_eb = 9999999.9_wp, & !< NAMELIST z0 (lsm_par) 179 z0h_eb = 9999999.9_wp !< NAMELIST z0h (lsm_par)194 z0h_eb = 9999999.9_wp !< NAMELIST z0h (lsm_par) 180 195 181 196 REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: & 182 ddz_soil, & !< 1/dz_soil 183 ddz_soil_stag, & !< 1/dz_soil_stag 184 dz_soil, & !< soil grid spacing (center-center) 185 dz_soil_stag, & !< soil grid spacing (edge-edge) 186 root_extr = 0.0_wp, & !< root extraction 197 ddz_soil_stag, & !< 1/dz_soil_stag 198 dz_soil_stag, & !< soil grid spacing (center-center) 199 root_extr = 0.0_wp, & !< root extraction 187 200 root_fraction = (/9999999.9_wp, 9999999.9_wp, & 188 201 9999999.9_wp, 9999999.9_wp /), & !< distribution of root surface area to the individual soil layers … … 191 204 192 205 REAL(wp), DIMENSION(nzb_soil:nzt_soil+1) :: & 193 soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp, 283.0_wp, & 194 283.0_wp /) !< soil temperature (K) 206 soil_temperature = (/290.0_wp, 287.0_wp, 285.0_wp, 283.0_wp, & !< soil temperature (K) 207 283.0_wp /), & 208 ddz_soil, & !< 1/dz_soil 209 dz_soil !< soil grid spacing (edge-edge) 195 210 196 211 #if defined( __nopointer ) … … 218 233 ! 219 234 !-- Energy balance variables 220 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: 235 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: & 221 236 alpha_vg, & !< coef. of Van Genuchten 222 237 c_liq, & !< liquid water coverage (of vegetated area) … … 232 247 lai, & !< leaf area index 233 248 lai_av, & !< average of lai 234 lambda_h_sat, & !< heat conductivity for dry soil 235 lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type) 236 lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type) 249 lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type) 250 lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type) 237 251 l_vg, & !< coef. of Van Genuchten 238 252 m_fc, & !< soil moisture at field capacity (m3/m3) … … 253 267 r_a_av, & !< avergae of r_a 254 268 r_canopy, & !< canopy resistance 255 r_soil, & !< soil resi tance269 r_soil, & !< soil resistance 256 270 r_soil_min, & !< minimum soil resistance 257 271 r_s, & !< total surface resistance (combination of r_soil and r_canopy) 258 r_s_av, & !< aver gae of r_s272 r_s_av, & !< average of r_s 259 273 r_canopy_min, & !< minimum canopy (stomatal) resistance 260 274 shf_eb, & !< surface flux of sensible heat 261 275 shf_eb_av !< average of shf_eb 276 262 277 263 278 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & … … 472 487 lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance, & 473 488 min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp, & 474 rho_lv, root_fraction, saturation_moisture, soil_moisture, & 475 soil_temperature, soil_type, soil_type_name, vegetation_coverage, & 476 veg_type, veg_type_name, wilting_point, z0_eb, z0h_eb 489 rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm, & 490 soil_moisture, soil_temperature, soil_type, soil_type_name, & 491 vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, & 492 z0h_eb 477 493 478 494 ! … … 483 499 zs 484 500 485 486 501 ! 487 502 !-- Public 2D output variables … … 490 505 qsws_soil_eb, qsws_soil_eb_av, qsws_veg_eb, qsws_veg_eb_av, & 491 506 r_a, r_a_av, r_s, r_s_av, shf_eb, shf_eb_av 492 493 507 494 508 ! … … 565 579 ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) ) 566 580 ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) ) 567 ALLOCATE ( lambda_h_sat(nysg:nyng,nxlg:nxrg) )568 581 ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) ) 569 582 ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) ) … … 612 625 INTEGER(iwp) :: k !< running index 613 626 627 REAL(wp) :: pt_1 !< potential temperature at first grid level 628 614 629 615 630 ! 616 631 !-- Calculate Exner function 617 exn 632 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 618 633 619 634 … … 674 689 ! 675 690 !-- Calculate grid spacings. Temperature and moisture are defined at 676 !-- the center of the soil layers, whereas gradients/fluxes are defined677 !-- at the edges (_stag)678 dz_soil _stag(nzb_soil) = zs(nzb_soil)679 680 DO k = nzb_soil+1, nzt_soil681 dz_soil _stag(k) = zs(k) - zs(k-1)691 !-- the edges of the soil layers (_stag), whereas gradients/fluxes are defined 692 !-- at the centers 693 dz_soil(nzb_soil) = zs(nzb_soil) 694 695 DO k = nzb_soil+1, nzt_soil 696 dz_soil(k) = zs(k) - zs(k-1) 682 697 ENDDO 683 684 DO k = nzb_soil, nzt_soil-1 685 dz_soil(k) = 0.5 * (dz_soil_stag(k+1) + dz_soil_stag(k)) 698 dz_soil(nzt_soil+1) = dz_soil(nzt_soil) 699 700 DO k = nzb_soil, nzt_soil-1 701 dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k)) 686 702 ENDDO 687 dz_soil (nzt_soil) = dz_soil_stag(nzt_soil)703 dz_soil_stag(nzt_soil) = dz_soil(nzt_soil) 688 704 689 705 ddz_soil = 1.0_wp / dz_soil … … 752 768 !-- (make sure that the soil moisture does not drop below the permanent 753 769 !-- wilting point) -> problems with devision by zero) 754 DO k = nzb_soil, nzt_soil770 DO k = nzb_soil, nzt_soil 755 771 t_soil(k,:,:) = soil_temperature(k) 756 772 m_soil(k,:,:) = MAX(soil_moisture(k),m_wilt(:,:)) … … 761 777 ! 762 778 !-- Calculate surface temperature 763 t_surface = pt_surface * exn 779 t_surface = pt_surface * exn 780 t_surface_p = t_surface 764 781 765 782 ! … … 770 787 k = nzb_s_inner(j,i) 771 788 789 IF ( cloud_physics ) THEN 790 pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 791 ELSE 792 pt_1 = pt(k+1,j,i) 793 ENDIF 794 772 795 ! 773 796 !-- Assure that r_a cannot be zero at model start 774 IF ( pt(k+1,j,i) == pt(k,j,i) ) pt(k+1,j,i) = pt(k+1,j,i) & 775 + 1.0E-10_wp 776 777 us(j,i) = 0.1_wp 778 ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i) 797 IF ( pt_1 == pt(k,j,i) ) pt_1 = pt_1 + 1.0E-10_wp 798 799 us(j,i) = 0.1_wp 800 ts(j,i) = (pt_1 - pt(k,j,i)) / r_a(j,i) 779 801 shf(j,i) = - us(j,i) * ts(j,i) 780 802 ENDDO … … 794 816 ENDIF 795 817 796 ! 797 !-- Calculate saturation soil heat conductivity 798 lambda_h_sat(:,:) = lambda_h_sm ** (1.0_wp - m_sat(:,:)) * & 799 lambda_h_water ** m_sat(:,:) 800 801 802 803 804 DO k = nzb_soil, nzt_soil 818 DO k = nzb_soil, nzt_soil 805 819 root_fr(k,:,:) = root_fraction(k) 806 820 ENDDO … … 838 852 839 853 IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN 840 DO k = nzb_soil, nzt_soil854 DO k = nzb_soil, nzt_soil 841 855 root_fr(k,:,:) = root_distribution(k,veg_type) 842 856 root_fraction(k) = root_distribution(k,veg_type) … … 886 900 887 901 ! 888 !-- Calculate humidity at the surface889 IF ( humidity ) THEN890 CALL calc_q_surface891 ENDIF892 893 894 895 !896 902 !-- Add timeseries for land surface model 897 898 dots_label(dots_num+1) = "ghf_eb"899 dots_label(dots_num+2) = "shf_eb"900 dots_label(dots_num+3) = "qsws_eb"901 dots_label(dots_num+4) = "qsws_liq_eb"902 dots_label(dots_num+5) = "qsws_soil_eb"903 dots_label(dots_num+6) = "qsws_veg_eb"904 dots_unit(dots_num+1:dots_num+6) = "W/m2"905 906 dots_label(dots_num+7) = "r_a"907 dots_label(dots_num+8) = "r_s"908 dots_unit(dots_num+7:dots_num+8) = "s/m"909 910 903 dots_soil = dots_num + 1 911 904 dots_num = dots_num + 8 905 906 dots_label(dots_soil) = "ghf_eb" 907 dots_label(dots_soil+1) = "shf_eb" 908 dots_label(dots_soil+2) = "qsws_eb" 909 dots_label(dots_soil+3) = "qsws_liq_eb" 910 dots_label(dots_soil+4) = "qsws_soil_eb" 911 dots_label(dots_soil+5) = "qsws_veg_eb" 912 dots_label(dots_soil+6) = "r_a" 913 dots_label(dots_soil+7) = "r_s" 914 915 dots_unit(dots_soil:dots_soil+5) = "W/m2" 916 dots_unit(dots_soil+6:dots_soil+7) = "s/m" 912 917 913 918 RETURN … … 921 926 ! ------------ 922 927 !> Solver for the energy balance at the surface. 923 !> @note Surface fluxes are calculated in the land surface model, but these are924 !> not used in the atmospheric code. The fluxes are calculated afterwards in925 !> prandtl_fluxes using the surface values of temperature and humidity as926 !> provided by the land surface model. In this way, the fluxes in the land927 !> surface model are not equal to the ones calculated in prandtl_fluxes928 928 !------------------------------------------------------------------------------! 929 929 SUBROUTINE lsm_energy_balance … … 940 940 f3, & !< resistance correction term 3 941 941 m_min, & !< minimum soil moisture 942 T_1, & !< actual temperature at first grid point943 942 e, & !< water vapour pressure 944 943 e_s, & !< water vapour saturation pressure … … 954 953 f_shf, & !< factor for shf_eb 955 954 lambda_surface, & !< Current value of lambda_surface 956 m_liq_eb_max !< maxmimum value of the liq. water reservoir 957 955 m_liq_eb_max, & !< maxmimum value of the liq. water reservoir 956 pt_1, & !< potential temperature at first grid level 957 qv_1 !< specific humidity at first grid level 958 958 959 959 ! … … 961 961 exn = ( surface_pressure / 1000.0_wp )**0.286_wp 962 962 963 964 965 DO i = nxlg, nxrg 966 DO j = nysg, nyng 963 DO i = nxlg, nxrg 964 DO j = nysg, nyng 967 965 k = nzb_s_inner(j,i) 968 966 … … 980 978 !-- time step is used here. Note that this formulation is the 981 979 !-- equivalent to the ECMWF formulation using drag coefficients 982 r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + & 983 1.0E-20_wp) 980 IF ( cloud_physics ) THEN 981 pt_1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i) 982 qv_1 = q(k+1,j,i) - ql(k+1,j,i) 983 ELSE 984 pt_1 = pt(k+1,j,i) 985 qv_1 = q(k+1,j,i) 986 ENDIF 987 988 r_a(j,i) = (pt_1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp) 984 989 985 990 ! … … 1002 1007 !-- f2 = 0 for very dry soils 1003 1008 m_total = 0.0_wp 1004 DO ks = nzb_soil, nzt_soil1009 DO ks = nzb_soil, nzt_soil 1005 1010 m_total = m_total + root_fr(ks,j,i) & 1006 1011 * MAX(m_soil(ks,j,i),m_wilt(j,i)) 1007 1012 ENDDO 1008 1013 1009 IF ( m_total .GT. m_wilt(j,i) .AND. m_total .LT.m_fc(j,i) ) THEN1014 IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN 1010 1015 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 1011 ELSEIF ( m_total .GE.m_fc(j,i) ) THEN1016 ELSEIF ( m_total >= m_fc(j,i) ) THEN 1012 1017 f2 = 1.0_wp 1013 1018 ELSE … … 1025 1030 ! 1026 1031 !-- Calculate vapour pressure 1027 e = q _p(k+1,j,i)* surface_pressure / 0.622_wp1032 e = qv_1 * surface_pressure / 0.622_wp 1028 1033 f3 = EXP ( -g_d(j,i) * (e_s - e) ) 1029 1034 ELSE … … 1067 1072 !-- In case of dewfall, set evapotranspiration to zero 1068 1073 !-- All super-saturated water is then removed from the air 1069 IF ( humidity .AND. dewfall .AND. q_s .LE. q_p(k+1,j,i)) THEN1074 IF ( humidity .AND. q_s <= qv_1 ) THEN 1070 1075 r_canopy(j,i) = 0.0_wp 1071 1076 r_soil(j,i) = 0.0_wp 1072 ENDIF 1077 ENDIF 1073 1078 1074 1079 ! … … 1085 1090 !-- If soil moisture is below wilting point, plants do no longer 1086 1091 !-- transpirate. 1087 ! IF ( m_soil(k,j,i) .LT.m_wilt(j,i) ) THEN1092 ! IF ( m_soil(k,j,i) < m_wilt(j,i) ) THEN 1088 1093 ! f_qsws_veg = 0.0_wp 1089 1094 ! ENDIF … … 1094 1099 !-- air. 1095 1100 IF ( humidity ) THEN 1096 IF ( q_s .LE. q_p(k+1,j,i)) THEN1101 IF ( q_s <= qv_1 ) THEN 1097 1102 IF ( .NOT. dewfall ) THEN 1098 1103 f_qsws_veg = 0.0_wp … … 1114 1119 dq_s_dt = 0.622_wp * e_s_dt / surface_pressure 1115 1120 1116 T_1 = pt_p(k+1,j,i) * exn1117 1118 1121 ! 1119 1122 !-- Add LW up so that it can be removed in prognostic equation 1120 rad_net_l(j,i) = rad_net(j,i) + sigma_sb * t_surface(j,i) ** 4 1123 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1124 1121 1125 1122 1126 IF ( humidity ) THEN 1123 1127 #if defined ( __rrtmg ) 1124 1128 ! 1125 1129 !-- Numerator of the prognostic equation 1126 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4& 1127 + f_shf / exn * T_1 + f_qsws * ( q_p(k+1,j,i) - q_s & 1130 coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb) & 1131 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1132 + f_shf * pt_1 + f_qsws * ( qv_1 - q_s & 1128 1133 + dq_s_dt * t_surface(j,i) ) + lambda_surface & 1129 1134 * t_soil(nzb_soil,j,i) … … 1131 1136 ! 1132 1137 !-- Denominator of the prognostic equation 1133 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws & 1134 * dq_s_dt + lambda_surface + f_shf / exn 1135 1138 coef_2 = rrtm_lwuflx_dt(0,nzb) + f_qsws * dq_s_dt & 1139 + lambda_surface + f_shf / exn 1140 #else 1141 1142 ! 1143 !-- Numerator of the prognostic equation 1144 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1145 * t_surface(j,i) ** 4 & 1146 + f_shf * pt_1 + f_qsws * ( qv_1 & 1147 - q_s + dq_s_dt * t_surface(j,i) ) & 1148 + lambda_surface * t_soil(nzb_soil,j,i) 1149 1150 ! 1151 !-- Denominator of the prognostic equation 1152 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 + f_qsws & 1153 * dq_s_dt + lambda_surface + f_shf / exn 1154 1155 #endif 1136 1156 ELSE 1137 1157 1158 #if defined ( __rrtmg ) 1138 1159 ! 1139 1160 !-- Numerator of the prognostic equation 1140 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb * t_surface(j,i) ** 4& 1141 + f_shf / exn * T_1 + lambda_surface & 1161 coef_1 = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb) & 1162 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1163 + f_shf * pt_1 + lambda_surface & 1142 1164 * t_soil(nzb_soil,j,i) 1165 1166 ! 1167 !-- Denominator of the prognostic equation 1168 coef_2 = rrtm_lwuflx_dt(0,nzb) + lambda_surface + f_shf / exn 1169 #else 1170 1171 ! 1172 !-- Numerator of the prognostic equation 1173 coef_1 = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1174 * t_surface(j,i) ** 4 + f_shf * pt_1 & 1175 + lambda_surface * t_soil(nzb_soil,j,i) 1143 1176 1144 1177 ! … … 1146 1179 coef_2 = 4.0_wp * sigma_sb * t_surface(j,i) ** 3 & 1147 1180 + lambda_surface + f_shf / exn 1148 1149 ENDIF 1181 #endif 1182 ENDIF 1183 1150 1184 1151 1185 tend = 0.0_wp … … 1159 1193 ! 1160 1194 !-- Add RK3 term 1161 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & 1162 * tt_surface_m(j,i) 1163 1164 ! 1165 !-- Calculate true tendency 1166 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & 1167 * tt_surface_m(j,i)) / (dt_3d * tsc(2)) 1168 ! 1169 !-- Calculate t_surface tendencies for the next Runge-Kutta step 1170 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1171 IF ( intermediate_timestep_count == 1 ) THEN 1172 tt_surface_m(j,i) = tend 1173 ELSEIF ( intermediate_timestep_count < & 1174 intermediate_timestep_count_max ) THEN 1175 tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1176 * tt_surface_m(j,i) 1195 IF ( c_surface /= 0.0_wp ) THEN 1196 1197 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & 1198 * tt_surface_m(j,i) 1199 1200 ! 1201 !-- Calculate true tendency 1202 tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3) & 1203 * tt_surface_m(j,i)) / (dt_3d * tsc(2)) 1204 ! 1205 !-- Calculate t_surface tendencies for the next Runge-Kutta step 1206 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1207 IF ( intermediate_timestep_count == 1 ) THEN 1208 tt_surface_m(j,i) = tend 1209 ELSEIF ( intermediate_timestep_count < & 1210 intermediate_timestep_count_max ) THEN 1211 tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1212 * tt_surface_m(j,i) 1213 ENDIF 1177 1214 ENDIF 1178 ENDIF 1179 1180 pt_p(k,j,i) = t_surface_p(j,i) / exn 1215 1216 ENDIF 1217 1218 ! 1219 !-- In case of fast changes in the skin temperature, it is required to 1220 !-- update the radiative fluxes in order to keep the solution stable 1221 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp ) THEN 1222 force_radiation_call_l = .TRUE. 1223 ENDIF 1224 1225 pt(k,j,i) = t_surface_p(j,i) / exn 1226 1181 1227 ! 1182 1228 !-- Calculate fluxes 1183 rad_net_l(j,i) = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1184 * t_surface(j,i)**4 - 4.0_wp * sigma_sb & 1185 * t_surface(j,i)**3 * t_surface_p(j,i) 1229 #if defined ( __rrtmg ) 1230 rad_net_l(j,i) = rad_net_l(j,i) + rrtm_lwuflx_dt(0,nzb) & 1231 * t_surface(j,i) - rad_lw_out(nzb,j,i) & 1232 - rrtm_lwuflx_dt(0,nzb) * t_surface_p(j,i) 1233 1234 IF ( rrtm_idrv == 1 ) THEN 1235 rad_net(j,i) = rad_net_l(j,i) 1236 rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) & 1237 + rrtm_lwuflx_dt(0,nzb) & 1238 * ( t_surface_p(j,i) - t_surface(j,i) ) 1239 ENDIF 1240 #else 1241 rad_net_l(j,i) = rad_net_l(j,i) + 3.0_wp * sigma_sb & 1242 * t_surface(j,i)**4 - 4.0_wp * sigma_sb & 1243 * t_surface(j,i)**3 * t_surface_p(j,i) 1244 rad_net(j,i) = rad_net_l(j,i) 1245 #endif 1246 1186 1247 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1187 1248 - t_soil(nzb_soil,j,i)) 1188 shf_eb(j,i) = - f_shf * ( pt_p(k+1,j,i) - pt_p(k,j,i) ) 1249 1250 shf_eb(j,i) = - f_shf * ( pt_1 - pt(k,j,i) ) 1251 1252 shf(j,i) = shf_eb(j,i) / rho_cp 1189 1253 1190 1254 IF ( humidity ) THEN 1191 qsws_eb(j,i) = - f_qsws * ( q _p(k+1,j,i) - q_s + dq_s_dt&1255 qsws_eb(j,i) = - f_qsws * ( qv_1 - q_s + dq_s_dt & 1192 1256 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) 1193 1257 1194 qsws_veg_eb(j,i) = - f_qsws_veg * ( q_p(k+1,j,i) - q_s & 1258 qsws(j,i) = qsws_eb(j,i) / rho_lv 1259 1260 qsws_veg_eb(j,i) = - f_qsws_veg * ( qv_1 - q_s & 1195 1261 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1196 1262 * t_surface_p(j,i) ) 1197 1263 1198 qsws_soil_eb(j,i) = - f_qsws_soil * ( q _p(k+1,j,i) - q_s&1264 qsws_soil_eb(j,i) = - f_qsws_soil * ( qv_1 - q_s & 1199 1265 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1200 1266 * t_surface_p(j,i) ) 1201 1267 1202 qsws_liq_eb(j,i) = - f_qsws_liq * ( q _p(k+1,j,i) - q_s&1268 qsws_liq_eb(j,i) = - f_qsws_liq * ( qv_1 - q_s & 1203 1269 + dq_s_dt * t_surface(j,i) - dq_s_dt & 1204 1270 * t_surface_p(j,i) ) 1205 1271 ENDIF 1206 1207 1272 ! 1208 1273 !-- Calculate the true surface resistance 1209 IF ( qsws_eb(j,i) .EQ.0.0_wp ) THEN1274 IF ( qsws_eb(j,i) == 0.0_wp ) THEN 1210 1275 r_s(j,i) = 1.0E10_wp 1211 1276 ELSE 1212 r_s(j,i) = - rho_lv * ( q _p(k+1,j,i) - q_s + dq_s_dt&1277 r_s(j,i) = - rho_lv * ( qv_1 - q_s + dq_s_dt & 1213 1278 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1214 1279 / qsws_eb(j,i) - r_a(j,i) … … 1228 1293 !-- Add precipitation to liquid water reservoir, if possible. 1229 1294 !-- Otherwise, add the water to bare soil fraction. 1230 IF ( m_liq_eb(j,i) .EQ.m_liq_eb_max ) THEN1295 IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN 1231 1296 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & 1232 + c_veg(j,i) * pr ecipitation_rate(j,i)&1297 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1233 1298 * 0.001_wp * rho_l * l_v 1234 1299 ELSE 1235 1300 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1236 + c_veg(j,i) * pr ecipitation_rate(j,i)&1301 + c_veg(j,i) * prr(k,j,i) * hyrho(k) & 1237 1302 * 0.001_wp * rho_l * l_v 1238 1303 ENDIF … … 1241 1306 !-- coverage. 1242 1307 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) * (1.0_wp & 1243 - c_veg(j,i)) * pr ecipitation_rate(j,i)&1308 - c_veg(j,i)) * prr(k,j,i) * hyrho(k) & 1244 1309 * 0.001_wp * rho_l * l_v 1245 1310 ENDIF 1311 1246 1312 ! 1247 1313 !-- If the air is saturated, check the reservoir water level 1248 IF ( q_s .LE. q_p(k+1,j,i)) THEN 1314 IF ( qsws_eb(j,i) < 0.0_wp ) THEN 1315 1249 1316 ! 1250 1317 !-- Check if reservoir is full (avoid values > m_liq_eb_max) … … 1252 1319 !-- case qsws_veg_eb is zero anyway (because c_liq = 1), 1253 1320 !-- so that tend is zero and no further check is needed 1254 IF ( m_liq_eb(j,i) .EQ.m_liq_eb_max ) THEN1321 IF ( m_liq_eb(j,i) == m_liq_eb_max ) THEN 1255 1322 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) & 1256 1323 + qsws_liq_eb(j,i) … … 1262 1329 !-- let the water enter the liquid water reservoir as dew on the 1263 1330 !-- plant 1264 IF ( qsws_veg_eb(j,i) .LT.0.0_wp ) THEN1331 IF ( qsws_veg_eb(j,i) < 0.0_wp ) THEN 1265 1332 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i) 1266 1333 qsws_veg_eb(j,i) = 0.0_wp … … 1301 1368 ENDDO 1302 1369 1370 #if defined( __parallel ) 1371 ! 1372 !-- Make a logical OR for all processes. Force radiation call if at 1373 !-- least one processor 1374 IF ( intermediate_timestep_count == intermediate_timestep_count_max-1 )& 1375 THEN 1376 1377 IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 1378 CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call, & 1379 1, MPI_LOGICAL, MPI_LOR, comm2d, ierr ) 1380 #else 1381 force_radiation_call = force_radiation_call_l 1382 #endif 1383 force_radiation_call_l = .FALSE. 1384 ENDIF 1385 1386 1303 1387 END SUBROUTINE lsm_energy_balance 1304 1388 … … 1322 1406 1323 1407 REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp, & !< temp. gamma 1324 lambda_temp, & !< temp. lambda1325 tend !< tendency1326 1327 DO i = nxlg, nxrg1328 DO j = nysg, nyng1329 DO k = nzb_soil, nzt_soil1408 lambda_temp, & !< temp. lambda 1409 tend !< tendency 1410 1411 DO i = nxlg, nxrg 1412 DO j = nysg, nyng 1413 DO k = nzb_soil, nzt_soil 1330 1414 ! 1331 1415 !-- Calculate volumetric heat capacity of the soil, taking into … … 1335 1419 1336 1420 ! 1337 !-- Calculate soil heat conductivity at the center of the soil 1421 !-- Calculate soil heat conductivity at the center of the soil 1338 1422 !-- layers 1423 lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) * & 1424 lambda_h_water ** m_soil(k,j,i) 1425 1339 1426 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) 1340 lambda_temp(k) = ke * (lambda_h_sat(j,i) + lambda_h_dry) + & 1427 1428 lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + & 1341 1429 lambda_h_dry 1342 1430 … … 1346 1434 !-- Calculate soil heat conductivity (lambda_h) at the _stag level 1347 1435 !-- using linear interpolation 1348 DO k = nzb_soil, nzt_soil-1 1349 1350 lambda_h(k,j,i) = lambda_temp(k) + & 1351 ( lambda_temp(k+1) - lambda_temp(k) ) & 1352 * 0.5 * dz_soil_stag(k) * ddz_soil(k+1) 1353 1436 DO k = nzb_soil, nzt_soil-1 1437 lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp 1354 1438 ENDDO 1355 1439 lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil) … … 1358 1442 !-- Prognostic equation for soil temperature t_soil 1359 1443 tend(:) = 0.0_wp 1360 tend( 0) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *&1444 tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * & 1361 1445 ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) & 1362 - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil )&1446 - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1) & 1363 1447 + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil) 1364 1448 1365 DO k = 1, nzt_soil 1449 1450 1451 1452 DO k = nzb_soil+1, nzt_soil 1366 1453 tend(k) = (1.0_wp/rho_c_total(k,j,i)) & 1367 1454 * ( lambda_h(k,j,i) & 1368 1455 * ( t_soil(k+1,j,i) - t_soil(k,j,i) ) & 1369 * ddz_soil(k )&1456 * ddz_soil(k+1) & 1370 1457 - lambda_h(k-1,j,i) & 1371 1458 * ( t_soil(k,j,i) - t_soil(k-1,j,i) ) & 1372 * ddz_soil(k -1)&1459 * ddz_soil(k) & 1373 1460 ) * ddz_soil_stag(k) 1374 1461 ENDDO … … 1396 1483 1397 1484 1398 DO k = nzb_soil, nzt_soil1485 DO k = nzb_soil, nzt_soil 1399 1486 1400 1487 ! … … 1414 1501 / (n_vg(j,i)-1.0_wp)) - 1.0_wp & 1415 1502 )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i) 1503 1416 1504 1417 1505 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & … … 1438 1526 !-- using linear interpolation. To do: replace this with 1439 1527 !-- ECMWF-IFS Eq. 8.81 1440 DO k = nzb_soil, nzt_soil-11528 DO k = nzb_soil, nzt_soil-1 1441 1529 1442 lambda_w(k,j,i) = lambda_temp(k) + & 1443 ( lambda_temp(k+1) - lambda_temp(k) ) & 1444 * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1) 1445 gamma_w(k,j,i) = gamma_temp(k) + & 1446 ( gamma_temp(k+1) - gamma_temp(k) ) & 1447 * 0.5_wp * dz_soil_stag(k) * ddz_soil(k+1) 1530 lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & 1531 * 0.5_wp 1532 gamma_w(k,j,i) = ( gamma_temp(k+1) + gamma_temp(k) ) & 1533 * 0.5_wp 1448 1534 1449 1535 ENDDO … … 1457 1543 gamma_w(nzt_soil,j,i) = 0.0_wp 1458 1544 ELSE 1459 gamma_w(nzt_soil,j,i) = lambda_temp(nzt_soil)1545 gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil) 1460 1546 ENDIF 1461 1547 … … 1470 1556 1471 1557 ! 1472 !-- Calculate the root extraction (ECMWF 7.69, modified so that the 1473 !-- sum of root_extr = 1). The energy balance solver guarantees a 1474 !-- positive transpiration, so that there is no need for an 1475 !-- additional check. 1476 1477 m_total = 0.0_wp 1478 DO k = nzb_soil, nzt_soil 1479 IF ( m_soil(k,j,i) .GT. m_wilt(j,i) ) THEN 1558 !-- Calculate the root extraction (ECMWF 7.69, the sum of root_extr 1559 !-- = 1). The energy balance solver guarantees a positive 1560 !-- transpiration, so that there is no need for an additional check. 1561 DO k = nzb_soil, nzt_soil 1562 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1480 1563 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) 1481 1564 ENDIF 1482 1565 ENDDO 1483 1566 1484 DO k = nzb_soil, nzt_soil 1485 IF ( m_soil(k,j,i) .GT. m_wilt(j,i) ) THEN 1486 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 1487 ELSE 1488 root_extr(k) = 0.0_wp 1489 ENDIF 1490 ENDDO 1567 IF ( m_total > 0.0_wp ) THEN 1568 DO k = nzb_soil, nzt_soil 1569 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1570 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 1571 ELSE 1572 root_extr(k) = 0.0_wp 1573 ENDIF 1574 ENDDO 1575 ENDIF 1491 1576 1492 1577 ! … … 1495 1580 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( & 1496 1581 m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) ) & 1497 * ddz_soil(nzb_soil ) - gamma_w(nzb_soil,j,i) - (&1582 * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( & 1498 1583 root_extr(nzb_soil) * qsws_veg_eb(j,i) & 1499 1584 + qsws_soil_eb(j,i) ) * drho_l_lv ) & … … 1502 1587 DO k = nzb_soil+1, nzt_soil-1 1503 1588 tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) & 1504 - m_soil(k,j,i) ) * ddz_soil(k) - gamma_w(k,j,i)&1505 - lambda_w(k-1,j,i) * (m_soil(k,j,i) -&1506 m_soil(k-1,j,i)) * ddz_soil(k-1)&1507 + gamma_w(k-1,j,i) - (root_extr(k)&1508 * qsws_veg_eb(j,i) * drho_l_lv)&1589 - m_soil(k,j,i) ) * ddz_soil(k+1) - gamma_w(k,j,i)& 1590 - lambda_w(k-1,j,i) * (m_soil(k,j,i) - & 1591 m_soil(k-1,j,i)) * ddz_soil(k) & 1592 + gamma_w(k-1,j,i) - (root_extr(k) & 1593 * qsws_veg_eb(j,i) * drho_l_lv) & 1509 1594 ) * ddz_soil_stag(k) 1510 1595 … … 1514 1599 * (m_soil(nzt_soil,j,i) & 1515 1600 - m_soil(nzt_soil-1,j,i)) & 1516 * ddz_soil(nzt_soil -1)&1601 * ddz_soil(nzt_soil) & 1517 1602 + gamma_w(nzt_soil-1,j,i) - ( & 1518 1603 root_extr(nzt_soil) & … … 1573 1658 REAL(wp) :: resistance !< aerodynamic and soil resistance term 1574 1659 1575 DO i = nxlg, nxrg1576 DO j = nysg, nyng1660 DO i = nxlg, nxrg 1661 DO j = nysg, nyng 1577 1662 k = nzb_s_inner(j,i) 1578 1663 1579 1664 ! 1580 1665 !-- Calculate water vapour pressure at saturation 1581 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i) & 1582 - 273.16_wp ) / ( t_surface(j,i) & 1583 - 35.86_wp ) ) 1666 e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i) & 1667 - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) ) 1584 1668 1585 1669 ! … … 1587 1671 q_s = 0.622_wp * e_s / surface_pressure 1588 1672 1589 1590 1673 resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i)) 1591 1674 1592 1675 ! 1593 1676 !-- Calculate specific humidity at surface 1594 q_p(k,j,i) = resistance * q_s + (1.0_wp - resistance) & 1595 * q_p(k+1,j,i) 1677 IF ( cloud_physics ) THEN 1678 q(k,j,i) = resistance * q_s + (1.0_wp - resistance) & 1679 * ( q(k+1,j,i) - ql(k+1,j,i) ) 1680 ELSE 1681 q(k,j,i) = resistance * q_s + (1.0_wp - resistance) & 1682 * q(k+1,j,i) 1683 ENDIF 1684 1685 ! 1686 !-- Update virtual potential temperature 1687 vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) ) 1596 1688 1597 1689 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
