Changeset 1788 for palm/trunk/SOURCE/land_surface_model.f90
- Timestamp:
- Mar 10, 2016 11:01:04 AM (8 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/land_surface_model.f90 (modified) (46 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/land_surface_model.f90
r1784 r1788 19 19 ! Current revisions: 20 20 ! ----------------- 21 ! Bugfix: calculate lambda_surface based on temperature gradient between skin 22 ! layer and soil layer instead of Obukhov length 23 ! Changed: moved calculation of surface specific humidity to energy balance solver 24 ! New: water surfaces are available by using a fixed sea surface temperature. 25 ! The roughness lengths are calculated dynamically using the Charnock 26 ! parameterization. This involves the new roughness length for moisture z0q. 27 ! New: modified solution of the energy balance solver and soil model for 28 ! paved surfaces (i.e. asphalt concrete). 29 ! Syntax layout improved. 30 ! Changed: parameter dewfall removed. 21 31 ! 22 !23 32 ! Former revisions: 24 33 ! ----------------- … … 105 114 !> @todo Consider partial absorption of the net shortwave radiation by the 106 115 !> skin layer. 107 !> @todo Allow for water surfaces116 !> @todo Improve surface water parameterization 108 117 !> @todo Invert indices (running from -3 to 0. Currently: nzb_soil=0, 109 118 !> nzt_soil=3)). … … 118 127 119 128 USE arrays_3d, & 120 ONLY: hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h 129 ONLY: hyp, ol, pt, pt_p, q, q_p, ql, qsws, shf, ts, us, vpt, z0, z0h, & 130 z0q 121 131 122 132 USE cloud_parameters, & … … 171 181 ! 172 182 !-- LSM variables 173 INTEGER(iwp) :: veg_type = 2, & !< vegetation type, 0: user-defined, 1-19: generic (see list) 174 soil_type = 3 !< soil type, 0: user-defined, 1-6: generic (see list) 183 INTEGER(iwp) :: veg_type = 2, & !< NAMELIST veg_type_2d 184 soil_type = 3 !< NAMELIST soil_type_2d 185 186 INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE :: soil_type_2d, & !< soil type, 0: user-defined, 1-7: generic (see list) 187 veg_type_2d !< vegetation type, 0: user-defined, 1-19: generic (see list) 188 189 LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface, & !< flag parameter for water surfaces (classes 14+15) 190 pave_surface, & !< flag parameter for pavements (asphalt etc.) (class 20) 191 building_surface !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet) 175 192 176 193 LOGICAL :: conserve_water_content = .TRUE., & !< open or closed bottom surface for the soil model 177 dewfall = .TRUE., & !< allow/inhibit dewfall178 194 force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls 179 195 land_surface = .FALSE. !< flag parameter indicating wheather the lsm is used … … 200 216 m_total = 0.0_wp, & !< weighted total water content of the soil (m3/m3) 201 217 n_vangenuchten = 9999999.9_wp, & !< NAMELIST n_vg 218 pave_depth = 9999999.9_wp, & !< depth of the pavement 219 pave_heat_capacity = 1.94E6_wp, & !< volumetric heat capacity of pavement (e.g. roads) 220 pave_heat_conductivity = 1.00_wp, & !< heat conductivity for pavements (e.g. roads) 202 221 q_s = 0.0_wp, & !< saturation specific humidity 203 222 residual_moisture = 9999999.9_wp, & !< NAMELIST m_res … … 210 229 wilting_point = 9999999.9_wp, & !< NAMELIST m_wilt 211 230 z0_eb = 9999999.9_wp, & !< NAMELIST z0 (lsm_par) 212 z0h_eb = 9999999.9_wp !< NAMELIST z0h (lsm_par) 231 z0h_eb = 9999999.9_wp, & !< NAMELIST z0h (lsm_par) 232 z0q_eb = 9999999.9_wp !< NAMELIST z0q (lsm_par) 213 233 214 234 REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: & … … 295 315 296 316 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: & 297 lambda_h, & !< heat conductivity of soil ( ?)317 lambda_h, & !< heat conductivity of soil (W/m/K) 298 318 lambda_w, & !< hydraulic diffusivity of soil (?) 299 gamma_w, & !< hydraulic conductivity of soil ( ?)319 gamma_w, & !< hydraulic conductivity of soil (W/m/K) 300 320 rho_c_total !< volumetric heat capacity of the actual soil matrix (?) 301 321 … … 327 347 ! 328 348 !-- Predefined Land surface classes (veg_type) 329 CHARACTER(26), DIMENSION(0: 19), PARAMETER :: veg_type_name = (/ &349 CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ & 330 350 'user defined ', & ! 0 331 351 'crops, mixed farming ', & ! 1 … … 347 367 'deciduous shrubs ', & ! 17 348 368 'mixed forest/woodland ', & ! 18 349 'interrupted forest ' & ! 19 369 'interrupted forest ', & ! 19 370 'pavements/roads ' & ! 20 350 371 /) 351 372 … … 368 389 !-- Land surface parameters I 369 390 !-- r_canopy_min, lai, c_veg, g_d 370 REAL(wp), DIMENSION(0:3,1: 19), PARAMETER :: veg_pars = RESHAPE( (/ &391 REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ & 371 392 180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp, & ! 1 372 393 110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp, & ! 2 … … 387 408 225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp, & ! 17 388 409 250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp, & ! 18 389 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp & ! 19 390 /), (/ 4, 19 /) ) 391 392 ! 393 !-- Land surface parameters II z0, z0h 394 REAL(wp), DIMENSION(0:1,1:19), PARAMETER :: roughness_par = RESHAPE( (/ & 395 0.25_wp, 0.25E-2_wp, & ! 1 396 0.20_wp, 0.20E-2_wp, & ! 2 397 2.00_wp, 2.00_wp, & ! 3 398 2.00_wp, 2.00_wp, & ! 4 399 2.00_wp, 2.00_wp, & ! 5 400 2.00_wp, 2.00_wp, & ! 6 401 0.47_wp, 0.47E-2_wp, & ! 7 402 0.013_wp, 0.013E-2_wp, & ! 8 403 0.034_wp, 0.034E-2_wp, & ! 9 404 0.5_wp, 0.50E-2_wp, & ! 10 405 0.17_wp, 0.17E-2_wp, & ! 11 406 1.3E-3_wp, 1.3E-4_wp, & ! 12 407 0.83_wp, 0.83E-2_wp, & ! 13 408 0.00_wp, 0.00E-2_wp, & ! 14 409 0.00_wp, 0.00E-2_wp, & ! 15 410 0.10_wp, 0.10E-2_wp, & ! 16 411 0.25_wp, 0.25E-2_wp, & ! 17 412 2.00_wp, 2.00E-2_wp, & ! 18 413 1.10_wp, 1.10E-2_wp & ! 19 414 /), (/ 2, 19 /) ) 410 175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp, & ! 19 411 0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp & ! 20 412 /), (/ 4, 20 /) ) 413 414 ! 415 !-- Land surface parameters II z0, z0h, z0q 416 REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ & 417 0.25_wp, 0.25E-2_wp, 0.25E-2_wp, & ! 1 418 0.20_wp, 0.20E-2_wp, 0.20E-2_wp, & ! 2 419 2.00_wp, 2.00_wp, 2.00_wp, & ! 3 420 2.00_wp, 2.00_wp, 2.00_wp, & ! 4 421 2.00_wp, 2.00_wp, 2.00_wp, & ! 5 422 2.00_wp, 2.00_wp, 2.00_wp, & ! 6 423 0.47_wp, 0.47E-2_wp, 0.47E-2_wp, & ! 7 424 0.013_wp, 0.013E-2_wp, 0.013E-2_wp, & ! 8 425 0.034_wp, 0.034E-2_wp, 0.034E-2_wp, & ! 9 426 0.5_wp, 0.50E-2_wp, 0.50E-2_wp, & ! 10 427 0.17_wp, 0.17E-2_wp, 0.17E-2_wp, & ! 11 428 1.3E-3_wp, 1.3E-4_wp, 1.3E-4_wp, & ! 12 429 0.83_wp, 0.83E-2_wp, 0.83E-2_wp, & ! 13 430 0.00_wp, 0.00_wp, 0.00_wp, & ! 14 431 0.00_wp, 0.00_wp, 0.00_wp, & ! 15 432 0.10_wp, 0.10E-2_wp, 0.10E-2_wp, & ! 16 433 0.25_wp, 0.25E-2_wp, 0.25E-2_wp, & ! 17 434 2.00_wp, 2.00E-2_wp, 2.00E-2_wp, & ! 18 435 1.10_wp, 1.10E-2_wp, 1.10E-2_wp, & ! 19 436 1.0E-4_wp, 1.0E-5_wp, 1.0E-5_wp & ! 20 437 /), (/ 3, 20 /) ) 415 438 416 439 ! 417 440 !-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in 418 REAL(wp), DIMENSION(0:2,1: 19), PARAMETER :: surface_pars = RESHAPE( (/ &441 REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ & 419 442 10.0_wp, 10.0_wp, 0.05_wp, & ! 1 420 443 10.0_wp, 10.0_wp, 0.05_wp, & ! 2 … … 430 453 58.0_wp, 58.0_wp, 0.00_wp, & ! 12 431 454 10.0_wp, 10.0_wp, 0.05_wp, & ! 13 432 1.0E 20_wp, 1.0E20_wp, 0.00_wp, & ! 14433 1.0E 20_wp, 1.0E20_wp, 0.00_wp, & ! 15455 1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 14 456 1.0E10_wp, 1.0E10_wp, 0.00_wp, & ! 15 434 457 10.0_wp, 10.0_wp, 0.05_wp, & ! 16 435 458 10.0_wp, 10.0_wp, 0.05_wp, & ! 17 436 459 20.0_wp, 15.0_wp, 0.03_wp, & ! 18 437 20.0_wp, 15.0_wp, 0.03_wp & ! 19 438 /), (/ 3, 19 /) ) 460 20.0_wp, 15.0_wp, 0.03_wp, & ! 19 461 0.0_wp, 0.0_wp, 0.00_wp & ! 20 462 /), (/ 3, 20 /) ) 439 463 440 464 ! 441 465 !-- Root distribution (sum = 1) level 1, level 2, level 3, level 4, 442 REAL(wp), DIMENSION(0:3,1: 19), PARAMETER :: root_distribution = RESHAPE( (/ &466 REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: root_distribution = RESHAPE( (/ & 443 467 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 1 444 468 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & ! 2 … … 459 483 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 17 460 484 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 18 461 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp & ! 19 462 /), (/ 4, 19 /) ) 485 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 19 486 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp & ! 20 487 /), (/ 4, 20 /) ) 463 488 464 489 ! … … 499 524 !-- Public parameters, constants and initial values 500 525 PUBLIC alpha_vangenuchten, c_surface, canopy_resistance_coefficient, & 501 conserve_water_content, dewfall, field_capacity,&526 conserve_water_content, field_capacity, & 502 527 f_shortwave_incoming, hydraulic_conductivity, init_lsm, & 503 528 init_lsm_arrays, lambda_surface_stable, lambda_surface_unstable, & 504 529 land_surface, leaf_area_index, lsm_energy_balance, lsm_soil_model, & 505 530 lsm_swap_timelevel, l_vangenuchten, min_canopy_resistance, & 506 min_soil_resistance, n_vangenuchten, residual_moisture, rho_cp, & 531 min_soil_resistance, n_vangenuchten, pave_heat_capacity, & 532 pave_depth, pave_heat_conductivity, residual_moisture, rho_cp, & 507 533 rho_lv, root_fraction, saturation_moisture, skip_time_do_lsm, & 508 534 soil_moisture, soil_temperature, soil_type, soil_type_name, & 509 535 vegetation_coverage, veg_type, veg_type_name, wilting_point, z0_eb, & 510 z0h_eb 536 z0h_eb, z0q_eb 511 537 512 538 ! … … 586 612 !-- Allocate 2D vegetation model arrays 587 613 ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) ) 614 ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) ) 588 615 ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) ) 589 616 ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) ) … … 601 628 ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) ) 602 629 ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) ) 630 ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) ) 603 631 ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) ) 604 632 ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) ) … … 613 641 ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) ) 614 642 ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) ) 643 ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) ) 644 ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) ) 645 ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) ) 615 646 616 647 #if ! defined( __nopointer ) … … 650 681 ! 651 682 !-- If no cloud physics is used, rho_surface has not been calculated before 652 IF ( .NOT.cloud_physics ) THEN683 IF ( .NOT. cloud_physics ) THEN 653 684 rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn ) 654 685 ENDIF … … 762 793 ENDIF 763 794 795 ! 796 !-- Map values to the respective 2D arrays 764 797 alpha_vg = alpha_vangenuchten 765 798 l_vg = l_vangenuchten … … 796 829 797 830 ! 798 !-- Set artifical values for ts and us so that r_a has its initial value for799 !-- the first time step831 !-- Set artifical values for ts and us so that r_a has its initial value 832 !-- for the first time step 800 833 DO i = nxlg, nxrg 801 834 DO j = nysg, nyng … … 864 897 z0h_eb = roughness_par(1,veg_type) 865 898 ENDIF 866 z0h_factor = z0h_eb / z0_eb 899 IF ( z0q_eb == 9999999.9_wp ) THEN 900 z0q_eb = roughness_par(2,veg_type) 901 ENDIF 902 z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp ) 867 903 868 904 IF ( ANY( root_fraction == 9999999.9_wp ) ) THEN … … 881 917 z0h_eb = z0_eb * z0h_factor 882 918 ENDIF 919 IF ( z0q_eb == 9999999.9_wp ) THEN 920 z0q_eb = z0_eb * z0h_factor 921 ENDIF 883 922 884 923 ENDIF 885 924 886 925 ! 887 !-- Initialize vegetation 926 !-- For surfaces covered with pavement, set depth of the pavement (with dry 927 !-- soil below). The depth must be greater than the first soil layer depth 928 IF ( veg_type == 20 ) THEN 929 IF ( pave_depth == 9999999.9_wp ) THEN 930 pave_depth = zs(nzb_soil) 931 ELSE 932 pave_depth = MAX( zs(nzb_soil), pave_depth ) 933 ENDIF 934 ENDIF 935 936 ! 937 !-- Map vegetation and soil types to 2D array to allow for heterogeneous 938 !-- surfaces via user interface see below 939 veg_type_2d = veg_type 940 soil_type_2d = soil_type 941 942 ! 943 !-- Map vegetation parameters to the respective 2D arrays 888 944 r_canopy_min = min_canopy_resistance 889 945 lai = leaf_area_index … … 895 951 z0 = z0_eb 896 952 z0h = z0h_eb 953 z0q = z0q_eb 897 954 898 955 ! … … 900 957 CALL user_init_land_surface 901 958 959 ! 960 !-- Set flag parameter if vegetation type was set to a water surface. Also 961 !-- set temperature to a constant value in all "soil" layers. 962 DO i = nxlg, nxrg 963 DO j = nysg, nyng 964 IF ( veg_type_2d(j,i) == 14 .OR. veg_type_2d(j,i) == 15 ) THEN 965 water_surface(j,i) = .TRUE. 966 ELSEIF ( veg_type_2d(j,i) == 20 ) THEN 967 pave_surface(j,i) = .TRUE. 968 m_soil(:,j,i) = 0.0_wp 969 ENDIF 970 971 ENDDO 972 ENDDO 973 974 ! 975 !-- Calculate new roughness lengths (for water surfaces only) 976 CALL calc_z0_water_surface 902 977 903 978 t_soil_p = t_soil … … 935 1010 INTEGER(iwp) :: k, ks !< running index 936 1011 937 REAL(wp) :: f1, & !< resistance correction term 1 1012 REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface 1013 f1, & !< resistance correction term 1 938 1014 f2, & !< resistance correction term 2 939 1015 f3, & !< resistance correction term 3 … … 965 1041 966 1042 ! 967 !-- Set lambda_surface according to stratification 968 IF ( ol(j,i) >= 0.0_wp ) THEN 969 lambda_surface = lambda_surface_s(j,i) 1043 !-- Set lambda_surface according to stratification between skin layer and soil 1044 IF ( .NOT. pave_surface(j,i) ) THEN 1045 1046 c_surface_tmp = c_surface 1047 1048 IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i)) THEN 1049 lambda_surface = lambda_surface_s(j,i) 1050 ELSE 1051 lambda_surface = lambda_surface_u(j,i) 1052 ENDIF 970 1053 ELSE 971 lambda_surface = lambda_surface_u(j,i) 1054 1055 c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp 1056 lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil) 1057 972 1058 ENDIF 973 1059 … … 1016 1102 ENDDO 1017 1103 1018 IF ( m_total > m_wilt(j,i) .AND.m_total < m_fc(j,i) ) THEN1104 IF ( m_total > m_wilt(j,i) .AND. m_total < m_fc(j,i) ) THEN 1019 1105 f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) ) 1020 1106 ELSEIF ( m_total >= m_fc(j,i) ) THEN … … 1050 1136 !-- Third step: calculate bare soil resistance r_soil. The Clapp & 1051 1137 !-- Hornberger parametrization does not consider c_veg. 1052 IF ( soil_type /= 7 ) THEN1138 IF ( soil_type_2d(j,i) /= 7 ) THEN 1053 1139 m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) * & 1054 1140 m_res(j,i) … … 1064 1150 1065 1151 ! 1066 !-- Calculate fraction of liquid water reservoir 1067 m_liq_eb_max = m_max_depth * lai(j,i) 1068 c_liq(j,i) = MIN(1.0_wp, m_liq_eb(j,i) / (m_liq_eb_max+1.0E-20_wp)) 1069 1152 !-- Calculate the maximum possible liquid water amount on plants and 1153 !-- bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is 1154 !-- assumed, while paved surfaces might hold up 1 mm of water. The 1155 !-- liquid water fraction for paved surfaces is calculated after 1156 !-- Noilhan & Planton (1989), while the ECMWF formulation is used for 1157 !-- vegetated surfaces and bare soils. 1158 IF ( pave_surface(j,i) ) THEN 1159 m_liq_eb_max = m_max_depth * 5.0_wp 1160 c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 ) 1161 ELSE 1162 m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i) & 1163 + (1.0_wp - c_veg(j,i)) ) 1164 c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max ) 1165 ENDIF 1070 1166 1071 1167 ! … … 1076 1172 !-- In case of dewfall, set evapotranspiration to zero 1077 1173 !-- All super-saturated water is then removed from the air 1078 IF ( humidity .AND.q_s <= qv1 ) THEN1174 IF ( humidity .AND. q_s <= qv1 ) THEN 1079 1175 r_canopy(j,i) = 0.0_wp 1080 1176 r_soil(j,i) = 0.0_wp … … 1083 1179 ! 1084 1180 !-- Calculate coefficients for the total evapotranspiration 1085 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1086 (r_a(j,i) + r_canopy(j,i)) 1087 1088 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1181 !-- In case of water surface, set vegetation and soil fluxes to zero. 1182 !-- For pavements, only evaporation of liquid water is possible. 1183 IF ( water_surface(j,i) ) THEN 1184 f_qsws_veg = 0.0_wp 1185 f_qsws_soil = 0.0_wp 1186 f_qsws_liq = rho_lv / r_a(j,i) 1187 ELSEIF ( pave_surface (j,i) ) THEN 1188 f_qsws_veg = 0.0_wp 1189 f_qsws_soil = 0.0_wp 1190 f_qsws_liq = rho_lv * c_liq(j,i) / r_a(j,i) 1191 ELSE 1192 f_qsws_veg = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/ & 1193 (r_a(j,i) + r_canopy(j,i)) 1194 f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) + & 1089 1195 r_soil(j,i)) 1090 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1091 1092 1196 f_qsws_liq = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i) 1197 ENDIF 1093 1198 ! 1094 1199 !-- If soil moisture is below wilting point, plants do no longer … … 1098 1203 ! ENDIF 1099 1204 1100 !1101 !-- If dewfall is deactivated, vegetation, soil and liquid water1102 !-- reservoir are not allowed to take up water from the super-saturated1103 !-- air.1104 IF ( humidity .AND. q_s <= qv1 .AND. .NOT. dewfall ) THEN1105 f_qsws_veg = 0.0_wp1106 f_qsws_soil = 0.0_wp1107 f_qsws_liq = 0.0_wp1108 ENDIF1109 1110 1205 f_shf = rho_cp / r_a(j,i) 1111 1206 f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq … … 1123 1218 rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i) 1124 1219 1220 ! 1221 !-- Calculate new skin temperature 1125 1222 IF ( humidity ) THEN 1126 1223 #if defined ( __rrtmg ) … … 1186 1283 !-- Implicit solution when the surface layer has no heat capacity, 1187 1284 !-- otherwise use RK3 scheme. 1188 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface * & 1189 t_surface(j,i) ) / ( c_surface + coef_2 * dt_3d& 1190 * tsc(2) ) 1285 t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp * & 1286 t_surface(j,i) ) / ( c_surface_tmp + coef_2 & 1287 * dt_3d * tsc(2) ) 1288 1191 1289 ! 1192 1290 !-- Add RK3 term 1193 IF ( c_surface /= 0.0_wp ) THEN1291 IF ( c_surface_tmp /= 0.0_wp ) THEN 1194 1292 1195 1293 t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3) & … … 1211 1309 ENDIF 1212 1310 ENDIF 1213 1214 1311 ENDIF 1215 1312 … … 1223 1320 !-- often reached, when no oscillations would occur (causes immense 1224 1321 !-- computing time for the radiation code). 1225 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND.&1322 IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp .AND. & 1226 1323 unscheduled_radiation_calls ) THEN 1227 1324 force_radiation_call_l = .TRUE. … … 1249 1346 #endif 1250 1347 1251 1252 1348 ghf_eb(j,i) = lambda_surface * (t_surface_p(j,i) & 1253 1349 - t_soil(nzb_soil,j,i)) … … 1275 1371 * t_surface_p(j,i) ) 1276 1372 ENDIF 1373 1277 1374 ! 1278 1375 !-- Calculate the true surface resistance … … 1280 1377 r_s(j,i) = 1.0E10_wp 1281 1378 ELSE 1282 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt &1379 r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt & 1283 1380 * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) ) & 1284 1381 / qsws_eb(j,i) - r_a(j,i) … … 1297 1394 ! 1298 1395 !-- Add precipitation to liquid water reservoir, if possible. 1299 !-- Otherwise, add the water to bare soil fraction. 1396 !-- Otherwise, add the water to soil. In case of 1397 !-- pavements, the exceeding water amount is implicitely removed 1398 !-- as runoff as qsws_soil_eb is then not used in the soil model 1300 1399 IF ( m_liq_eb(j,i) /= m_liq_eb_max ) THEN 1301 1400 qsws_liq_eb(j,i) = qsws_liq_eb(j,i) & … … 1376 1475 !-- Make a logical OR for all processes. Force radiation call if at 1377 1476 !-- least one processor reached the threshold change in skin temperature 1378 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count&1477 IF ( unscheduled_radiation_calls .AND. intermediate_timestep_count & 1379 1478 == intermediate_timestep_count_max-1 ) THEN 1380 1479 #if defined( __parallel ) … … 1388 1487 ENDIF 1389 1488 1489 ! 1490 !-- Calculate surface specific humidity 1491 IF ( humidity ) THEN 1492 CALL calc_q_surface 1493 ENDIF 1494 1495 ! 1496 !-- Calculate new roughness lengths (for water surfaces only) 1497 CALL calc_z0_water_surface 1498 1390 1499 1391 1500 END SUBROUTINE lsm_energy_balance … … 1415 1524 DO i = nxlg, nxrg 1416 1525 DO j = nysg, nyng 1417 DO k = nzb_soil, nzt_soil 1418 ! 1419 !-- Calculate volumetric heat capacity of the soil, taking into 1420 !-- account water content 1421 rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) & 1422 + rho_c_water * m_soil(k,j,i)) 1423 1424 ! 1425 !-- Calculate soil heat conductivity at the center of the soil 1426 !-- layers 1427 lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) * & 1428 lambda_h_water ** m_soil(k,j,i) 1429 1430 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) 1431 1432 lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + & 1433 lambda_h_dry 1434 1435 ENDDO 1436 1437 ! 1438 !-- Calculate soil heat conductivity (lambda_h) at the _stag level 1439 !-- using linear interpolation 1440 DO k = nzb_soil, nzt_soil-1 1441 lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp 1442 ENDDO 1443 lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil) 1444 1445 ! 1446 !-- Prognostic equation for soil temperature t_soil 1447 tend(:) = 0.0_wp 1448 tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * & 1449 ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) & 1450 - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1) & 1451 + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil) 1452 1453 1454 1455 1456 DO k = nzb_soil+1, nzt_soil 1457 tend(k) = (1.0_wp/rho_c_total(k,j,i)) & 1458 * ( lambda_h(k,j,i) & 1459 * ( t_soil(k+1,j,i) - t_soil(k,j,i) ) & 1460 * ddz_soil(k+1) & 1461 - lambda_h(k-1,j,i) & 1462 * ( t_soil(k,j,i) - t_soil(k-1,j,i) ) & 1463 * ddz_soil(k) & 1464 ) * ddz_soil_stag(k) 1465 ENDDO 1466 1467 t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i) & 1468 + dt_3d * ( tsc(2) & 1469 * tend(:) + tsc(3) & 1470 * tt_soil_m(:,j,i) ) 1471 1472 ! 1473 !-- Calculate t_soil tendencies for the next Runge-Kutta step 1474 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1475 IF ( intermediate_timestep_count == 1 ) THEN 1476 DO k = nzb_soil, nzt_soil 1477 tt_soil_m(k,j,i) = tend(k) 1478 ENDDO 1479 ELSEIF ( intermediate_timestep_count < & 1480 intermediate_timestep_count_max ) THEN 1481 DO k = nzb_soil, nzt_soil 1482 tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & 1483 * tt_soil_m(k,j,i) 1484 ENDDO 1485 ENDIF 1486 ENDIF 1487 1488 1489 DO k = nzb_soil, nzt_soil 1490 1491 ! 1492 !-- Calculate soil diffusivity at the center of the soil layers 1493 lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat & 1494 / m_sat(j,i) ) * ( MAX(m_soil(k,j,i), & 1495 m_wilt(j,i)) / m_sat(j,i) )**(b_ch + 2.0_wp) 1496 1497 ! 1498 !-- Parametrization of Van Genuchten 1499 IF ( soil_type /= 7 ) THEN 1500 ! 1501 !-- Calculate the hydraulic conductivity after Van Genuchten 1502 !-- (1980) 1503 h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & 1504 MAX(m_soil(k,j,i),m_wilt(j,i)) ) )**(n_vg(j,i) & 1505 / (n_vg(j,i)-1.0_wp)) - 1.0_wp & 1506 )**(1.0_wp/n_vg(j,i)) / alpha_vg(j,i) 1507 1508 1509 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & 1510 (alpha_vg(j,i)*h_vg)**n_vg(j,i))**(1.0_wp & 1511 -1.0_wp/n_vg(j,i)) - (alpha_vg(j,i)*h_vg & 1512 )**(n_vg(j,i)-1.0_wp))**2 ) & 1513 / ( (1.0_wp + (alpha_vg(j,i)*h_vg & 1514 )**n_vg(j,i))**((1.0_wp - 1.0_wp/n_vg(j,i)) & 1515 *(l_vg(j,i) + 2.0_wp)) ) 1516 1517 ! 1518 !-- Parametrization of Clapp & Hornberger 1519 ELSE 1520 gamma_temp(k) = gamma_w_sat(j,i) * (m_soil(k,j,i) & 1521 / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp) 1522 ENDIF 1523 1524 ENDDO 1525 1526 1527 IF ( humidity ) THEN 1528 ! 1529 !-- Calculate soil diffusivity (lambda_w) at the _stag level 1530 !-- using linear interpolation. To do: replace this with 1531 !-- ECMWF-IFS Eq. 8.81 1526 1527 IF ( pave_surface(j,i) ) THEN 1528 rho_c_total(nzb_soil,j,i) = pave_heat_capacity 1529 lambda_temp(nzb_soil) = pave_heat_conductivity 1530 ENDIF 1531 1532 IF ( .NOT. water_surface(j,i) ) THEN 1533 DO k = nzb_soil, nzt_soil 1534 1535 1536 IF ( pave_surface(j,i) .AND. zs(k) <= pave_depth ) THEN 1537 1538 rho_c_total(k,j,i) = pave_heat_capacity 1539 lambda_temp(k) = pave_heat_conductivity 1540 1541 ELSE 1542 ! 1543 !-- Calculate volumetric heat capacity of the soil, taking 1544 !-- into account water content 1545 rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) & 1546 + rho_c_water * m_soil(k,j,i)) 1547 1548 ! 1549 !-- Calculate soil heat conductivity at the center of the soil 1550 !-- layers 1551 lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) * & 1552 lambda_h_water ** m_soil(k,j,i) 1553 1554 ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i) / m_sat(j,i))) 1555 1556 lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) + & 1557 lambda_h_dry 1558 ENDIF 1559 1560 ENDDO 1561 1562 ! 1563 !-- Calculate soil heat conductivity (lambda_h) at the _stag level 1564 !-- using linear interpolation. For pavement surface, the 1565 !-- true pavement depth is considered 1532 1566 DO k = nzb_soil, nzt_soil-1 1533 1534 lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & 1535 * 0.5_wp 1536 gamma_w(k,j,i) = ( gamma_temp(k+1) + gamma_temp(k) ) & 1537 * 0.5_wp 1538 1567 IF ( pave_surface(j,i) .AND. zs(k) < pave_depth & 1568 .AND. zs(k+1) > pave_depth ) THEN 1569 lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1) & 1570 * lambda_temp(k) & 1571 + ( 1.0_wp - ( pave_depth - zs(k) ) & 1572 / dz_soil(k+1) ) * lambda_temp(k+1) 1573 ELSE 1574 lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & 1575 * 0.5_wp 1576 ENDIF 1539 1577 ENDDO 1540 1541 ! 1542 ! 1543 !-- In case of a closed bottom (= water content is conserved), set 1544 !-- hydraulic conductivity to zero to that no water will be lost 1545 !-- in the bottom layer. 1546 IF ( conserve_water_content ) THEN 1547 gamma_w(nzt_soil,j,i) = 0.0_wp 1548 ELSE 1549 gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil) 1550 ENDIF 1551 1552 !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v)) 1553 !-- ensures the mass conservation for water. The transpiration of 1554 !-- plants equals the cumulative withdrawals by the roots in the 1555 !-- soil. The scheme takes into account the availability of water 1556 !-- in the soil layers as well as the root fraction in the 1557 !-- respective layer. Layer with moisture below wilting point will 1558 !-- not contribute, which reflects the preference of plants to 1559 !-- take water from moister layers. 1560 1561 ! 1562 !-- Calculate the root extraction (ECMWF 7.69, the sum of root_extr 1563 !-- = 1). The energy balance solver guarantees a positive 1564 !-- transpiration, so that there is no need for an additional check. 1565 DO k = nzb_soil, nzt_soil 1566 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1567 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) 1568 ENDIF 1569 ENDDO 1570 1571 IF ( m_total > 0.0_wp ) THEN 1572 DO k = nzb_soil, nzt_soil 1573 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1574 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 1575 ELSE 1576 root_extr(k) = 0.0_wp 1577 ENDIF 1578 ENDDO 1579 ENDIF 1580 1581 ! 1582 !-- Prognostic equation for soil water content m_soil 1578 lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil) 1579 1580 1581 1582 1583 ! 1584 !-- Prognostic equation for soil temperature t_soil 1583 1585 tend(:) = 0.0_wp 1584 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( & 1585 m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )&1586 * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - (&1587 root_extr(nzb_soil) * qsws_veg_eb(j,i)&1588 + qsws_soil_eb(j,i) ) * drho_l_lv ) &1589 * ddz_soil_stag(nzb_soil) 1590 1591 DO k = nzb_soil+1, nzt_soil-11592 tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)&1593 - m_soil(k,j,i) ) * ddz_soil(k+1) - gamma_w(k,j,i)&1594 - lambda_w(k-1,j,i) * (m_soil(k,j,i) -&1595 m_soil(k-1,j,i)) * ddz_soil(k)&1596 + gamma_w(k-1,j,i) - (root_extr(k)&1597 * qsws_veg_eb(j,i) * drho_l_lv)&1598 ) * ddz_soil_stag(k)1586 1587 tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) * & 1588 ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i) & 1589 - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1) & 1590 + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil) 1591 1592 DO k = nzb_soil+1, nzt_soil 1593 tend(k) = (1.0_wp/rho_c_total(k,j,i)) & 1594 * ( lambda_h(k,j,i) & 1595 * ( t_soil(k+1,j,i) - t_soil(k,j,i) ) & 1596 * ddz_soil(k+1) & 1597 - lambda_h(k-1,j,i) & 1598 * ( t_soil(k,j,i) - t_soil(k-1,j,i) ) & 1599 * ddz_soil(k) & 1600 ) * ddz_soil_stag(k) 1599 1601 1600 1602 ENDDO 1601 tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i) & 1602 - lambda_w(nzt_soil-1,j,i) & 1603 * (m_soil(nzt_soil,j,i) & 1604 - m_soil(nzt_soil-1,j,i)) & 1605 * ddz_soil(nzt_soil) & 1606 + gamma_w(nzt_soil-1,j,i) - ( & 1607 root_extr(nzt_soil) & 1608 * qsws_veg_eb(j,i) * drho_l_lv ) & 1609 ) * ddz_soil_stag(nzt_soil) 1610 1611 m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)& 1612 + dt_3d * ( tsc(2) * tend(:) & 1613 + tsc(3) * tm_soil_m(:,j,i) ) 1614 1615 ! 1616 !-- Account for dry soils (find a better solution here!) 1617 DO k = nzb_soil, nzt_soil 1618 IF ( m_soil_p(k,j,i) < 0.0_wp ) m_soil_p(k,j,i) = 0.0_wp 1619 ENDDO 1620 1621 1622 1623 1624 ! 1625 !-- Calculate m_soil tendencies for the next Runge-Kutta step 1603 1604 t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)& 1605 + dt_3d * ( tsc(2) & 1606 * tend(nzb_soil:nzt_soil) & 1607 + tsc(3) & 1608 * tt_soil_m(:,j,i) ) 1609 1610 ! 1611 !-- Calculate t_soil tendencies for the next Runge-Kutta step 1626 1612 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1627 1613 IF ( intermediate_timestep_count == 1 ) THEN 1628 1614 DO k = nzb_soil, nzt_soil 1629 t m_soil_m(k,j,i) = tend(k)1615 tt_soil_m(k,j,i) = tend(k) 1630 1616 ENDDO 1631 1617 ELSEIF ( intermediate_timestep_count < & 1632 1618 intermediate_timestep_count_max ) THEN 1633 1619 DO k = nzb_soil, nzt_soil 1634 t m_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp &1635 * tm_soil_m(k,j,i)1620 tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp & 1621 * tt_soil_m(k,j,i) 1636 1622 ENDDO 1637 1623 ENDIF 1638 1624 ENDIF 1639 1625 1626 1627 DO k = nzb_soil, nzt_soil 1628 1629 ! 1630 !-- Calculate soil diffusivity at the center of the soil layers 1631 lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat & 1632 / m_sat(j,i) ) * ( MAX( m_soil(k,j,i), & 1633 m_wilt(j,i) ) / m_sat(j,i) )**( & 1634 b_ch + 2.0_wp ) 1635 1636 ! 1637 !-- Parametrization of Van Genuchten 1638 IF ( soil_type /= 7 ) THEN 1639 ! 1640 !-- Calculate the hydraulic conductivity after Van Genuchten 1641 !-- (1980) 1642 h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) - & 1643 MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**( & 1644 n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp & 1645 )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i) 1646 1647 1648 gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp + & 1649 ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**( & 1650 1.0_wp - 1.0_wp / n_vg(j,i) ) - ( & 1651 alpha_vg(j,i) * h_vg )**( n_vg(j,i) & 1652 - 1.0_wp) )**2 ) & 1653 / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg & 1654 )**n_vg(j,i) )**( ( 1.0_wp - 1.0_wp & 1655 / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) ) 1656 1657 ! 1658 !-- Parametrization of Clapp & Hornberger 1659 ELSE 1660 gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i) & 1661 / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp) 1662 ENDIF 1663 1664 ENDDO 1665 1666 ! 1667 !-- Prognostic equation for soil moisture content. Only performed, 1668 !-- when humidity is enabled in the atmosphere and the surface type 1669 !-- is not pavement (implies dry soil below). 1670 IF ( humidity .AND. .NOT. pave_surface(j,i) ) THEN 1671 ! 1672 !-- Calculate soil diffusivity (lambda_w) at the _stag level 1673 !-- using linear interpolation. To do: replace this with 1674 !-- ECMWF-IFS Eq. 8.81 1675 DO k = nzb_soil, nzt_soil-1 1676 1677 lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) ) & 1678 * 0.5_wp 1679 gamma_w(k,j,i) = ( gamma_temp(k+1) + gamma_temp(k) ) & 1680 * 0.5_wp 1681 1682 ENDDO 1683 1684 ! 1685 ! 1686 !-- In case of a closed bottom (= water content is conserved), set 1687 !-- hydraulic conductivity to zero to that no water will be lost 1688 !-- in the bottom layer. 1689 IF ( conserve_water_content ) THEN 1690 gamma_w(nzt_soil,j,i) = 0.0_wp 1691 ELSE 1692 gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil) 1693 ENDIF 1694 1695 !-- The root extraction (= root_extr * qsws_veg_eb / (rho_l * l_v)) 1696 !-- ensures the mass conservation for water. The transpiration of 1697 !-- plants equals the cumulative withdrawals by the roots in the 1698 !-- soil. The scheme takes into account the availability of water 1699 !-- in the soil layers as well as the root fraction in the 1700 !-- respective layer. Layer with moisture below wilting point will 1701 !-- not contribute, which reflects the preference of plants to 1702 !-- take water from moister layers. 1703 1704 ! 1705 !-- Calculate the root extraction (ECMWF 7.69, the sum of root_extr 1706 !-- = 1). The energy balance solver guarantees a positive 1707 !-- transpiration, so that there is no need for an additional check. 1708 DO k = nzb_soil, nzt_soil 1709 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1710 m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i) 1711 ENDIF 1712 ENDDO 1713 1714 IF ( m_total > 0.0_wp ) THEN 1715 DO k = nzb_soil, nzt_soil 1716 IF ( m_soil(k,j,i) > m_wilt(j,i) ) THEN 1717 root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i) / m_total 1718 ELSE 1719 root_extr(k) = 0.0_wp 1720 ENDIF 1721 ENDDO 1722 ENDIF 1723 1724 ! 1725 !-- Prognostic equation for soil water content m_soil. 1726 tend(:) = 0.0_wp 1727 1728 tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * ( & 1729 m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) ) & 1730 * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( & 1731 root_extr(nzb_soil) * qsws_veg_eb(j,i) & 1732 + qsws_soil_eb(j,i) ) * drho_l_lv ) & 1733 * ddz_soil_stag(nzb_soil) 1734 1735 DO k = nzb_soil+1, nzt_soil-1 1736 tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i) & 1737 - m_soil(k,j,i) ) * ddz_soil(k+1) & 1738 - gamma_w(k,j,i) & 1739 - lambda_w(k-1,j,i) * (m_soil(k,j,i) - & 1740 m_soil(k-1,j,i)) * ddz_soil(k) & 1741 + gamma_w(k-1,j,i) - (root_extr(k) & 1742 * qsws_veg_eb(j,i) * drho_l_lv) & 1743 ) * ddz_soil_stag(k) 1744 1745 ENDDO 1746 tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i) & 1747 - lambda_w(nzt_soil-1,j,i) & 1748 * (m_soil(nzt_soil,j,i) & 1749 - m_soil(nzt_soil-1,j,i)) & 1750 * ddz_soil(nzt_soil) & 1751 + gamma_w(nzt_soil-1,j,i) - ( & 1752 root_extr(nzt_soil) & 1753 * qsws_veg_eb(j,i) * drho_l_lv ) & 1754 ) * ddz_soil_stag(nzt_soil) 1755 1756 m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)& 1757 + dt_3d * ( tsc(2) * tend(:) & 1758 + tsc(3) * tm_soil_m(:,j,i) ) 1759 1760 ! 1761 !-- Account for dry soils (find a better solution here!) 1762 DO k = nzb_soil, nzt_soil 1763 IF ( m_soil_p(k,j,i) < 0.0_wp ) m_soil_p(k,j,i) = 0.0_wp 1764 ENDDO 1765 1766 ! 1767 !-- Calculate m_soil tendencies for the next Runge-Kutta step 1768 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1769 IF ( intermediate_timestep_count == 1 ) THEN 1770 DO k = nzb_soil, nzt_soil 1771 tm_soil_m(k,j,i) = tend(k) 1772 ENDDO 1773 ELSEIF ( intermediate_timestep_count < & 1774 intermediate_timestep_count_max ) THEN 1775 DO k = nzb_soil, nzt_soil 1776 tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp& 1777 * tm_soil_m(k,j,i) 1778 ENDDO 1779 ENDIF 1780 ENDIF 1781 1782 ENDIF 1783 1640 1784 ENDIF 1641 1785 1642 1786 ENDDO 1643 1787 ENDDO 1644 1645 !1646 !-- Calculate surface specific humidity1647 IF ( humidity ) THEN1648 CALL calc_q_surface1649 ENDIF1650 1651 1788 1652 1789 END SUBROUTINE lsm_soil_model … … 1656 1793 ! Description: 1657 1794 ! ------------ 1658 !> Calculation of specific humidity of the surface layer (surface) 1795 !> Calculation of roughness length for open water (lakes, ocean). The 1796 !> parameterization follows Charnock (1955). Two different implementations 1797 !> are available: as in ECMWF-IFS (Beljaars 1994) or as in FLake (Subin et al. 1798 !> 2012) 1799 !------------------------------------------------------------------------------! 1800 SUBROUTINE calc_z0_water_surface 1801 1802 USE control_parameters, & 1803 ONLY: g, kappa, molecular_viscosity 1804 1805 IMPLICIT NONE 1806 1807 INTEGER :: i !< running index 1808 INTEGER :: j !< running index 1809 1810 REAL(wp), PARAMETER :: alpha_ch = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF 1811 ! REAL(wp), PARAMETER :: pr_number = 0.71_wp !< molecular Prandtl number in the Charnock parameterization (differs from prandtl_number) 1812 ! REAL(wp), PARAMETER :: sc_number = 0.66_wp !< molecular Schmidt number in the Charnock parameterization 1813 ! REAL(wp) :: re_0 !< near-surface roughness Reynolds number 1814 1815 1816 DO i = nxlg, nxrg 1817 DO j = nysg, nyng 1818 IF ( water_surface(j,i) ) THEN 1819 1820 ! 1821 !-- Disabled: FLake parameterization. Ideally, the Charnock 1822 !-- coefficient should depend on the water depth and the fetch 1823 !-- length 1824 ! re_0 = z0(j,i) * us(j,i) / molecular_viscosity 1825 ! 1826 ! z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i), & 1827 ! alpha_ch * us(j,i) / g ) 1828 ! 1829 ! z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) ) 1830 ! z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) ) 1831 1832 ! 1833 !-- Set minimum roughness length for u* > 0.2 1834 ! IF ( us(j,i) > 0.2_wp ) THEN 1835 ! z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) ) 1836 ! z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) ) 1837 ! ENDIF 1838 1839 ! 1840 !-- ECMWF IFS model parameterization after Beljaars (1994). At low 1841 !-- wind speed, the sea surface becomes aerodynamically smooth and 1842 !-- the roughness scales with the viscosity. At high wind speed, the 1843 !-- Charnock relation is used. 1844 z0(j,i) = ( 0.11_wp * molecular_viscosity / us(j,i) ) & 1845 + ( alpha_ch * us(j,i)**2 / g ) 1846 1847 z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i) 1848 z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i) 1849 1850 ENDIF 1851 ENDDO 1852 ENDDO 1853 1854 END SUBROUTINE calc_z0_water_surface 1855 1856 1857 !------------------------------------------------------------------------------! 1858 ! Description: 1859 ! ------------ 1860 !> Calculation of specific humidity of the skin layer (surface). It is assumend 1861 !> that the skin is always saturated. 1659 1862 !------------------------------------------------------------------------------! 1660 1863 SUBROUTINE calc_q_surface … … 1665 1868 INTEGER :: j !< running index 1666 1869 INTEGER :: k !< running index 1870 1667 1871 REAL(wp) :: resistance !< aerodynamic and soil resistance term 1668 1872 … … 1701 1905 END SUBROUTINE calc_q_surface 1702 1906 1907 1703 1908 !------------------------------------------------------------------------------! 1704 1909 ! Description:
Note: See TracChangeset
for help on using the changeset viewer.
