Changeset 1500 for palm/trunk
- Timestamp:
- Dec 3, 2014 5:42:41 PM (9 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
Makefile (modified) (4 diffs)
-
check_parameters.f90 (modified) (1 diff)
-
land_surface_model.f90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/Makefile
r1497 r1500 20 20 # Current revisions: 21 21 # ------------------ 22 # 22 # Bugfix: missing adjustments for land surface model and radiation model 23 23 # 24 24 # Former revisions: … … 298 298 check_open.o: modules.o mod_kinds.o mod_particle_attributes.o 299 299 check_parameters.o: modules.o mod_kinds.o subsidence.o land_surface_model.o\ 300 plant_canopy_model.o 300 plant_canopy_model.o radiation_model.o 301 301 close_file.o: modules.o mod_kinds.o 302 302 compute_vpt.o: modules.o mod_kinds.o … … 398 398 poisfft.o: modules.o cpulog.o fft_xy.o mod_kinds.o tridia_solver.o 399 399 poismg.o: modules.o cpulog.o mod_kinds.o 400 prandtl_fluxes.o: modules.o mod_kinds.o 400 prandtl_fluxes.o: modules.o mod_kinds.o land_surface_model.o 401 401 pres.o: modules.o cpulog.o mod_kinds.o poisfft.o 402 402 print_1d.o: modules.o cpulog.o mod_kinds.o … … 422 422 sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o 423 423 surface_coupler.o: modules.o cpulog.o mod_kinds.o 424 swap_timelevel.o: modules.o cpulog.o mod_kinds.o 424 swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model.o 425 425 temperton_fft.o: modules.o mod_kinds.o 426 426 time_integration.o: modules.o advec_ws.o buoyancy.o calc_mean_profile.o \ -
palm/trunk/SOURCE/check_parameters.f90
r1497 r1500 958 958 ENDIF 959 959 960 ! Dirichlet boundary conditions are allowed at the moment for testing 961 ! purposes 962 ! IF ( bc_pt_b == 'dirichlet' .OR. bc_q_b == 'dirichlet' ) THEN 963 ! message_string = 'lsm requires setting of'// & 964 ! 'bc_pt_b = "neumann" and '// & 965 ! 'bc_q_b = "neumann"' 966 ! CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 ) 967 ! ENDIF 960 ! 961 !-- Dirichlet boundary conditions are required as the surface fluxes are 962 !-- calculated from the temperature/humidity gradients in the land surface 963 !-- model 964 IF ( bc_pt_b == 'neumann' .OR. bc_q_b == 'neumann' ) THEN 965 message_string = 'lsm requires setting of'// & 966 'bc_pt_b = "dirichlet" and '// & 967 'bc_q_b = "dirichlet"' 968 CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 ) 969 ENDIF 968 970 969 971 IF ( .NOT. prandtl_layer ) THEN -
palm/trunk/SOURCE/land_surface_model.f90
r1497 r1500 20 20 ! Current revisions: 21 21 ! ----------------- 22 ! 22 ! Corrected calculation of aerodynamic resistance (r_a). 23 ! Precipitation is now added to liquid water reservoir using LE_liq. 24 ! Added support for dry runs. 23 25 ! 24 26 ! Former revisions: … … 37 39 ! H-TESSEL. The implementation is based on the formulation implemented in the 38 40 ! DALES model. 41 ! 42 ! To do list: 43 ! ----------- 44 ! - Add support for binary I/O support 45 ! - Add support for lsm data output 46 ! - Check for time step criterion 47 ! - Check use with RK-2 and Euler time-stepping 48 ! - Adaption for use with cloud physics (liquid water potential temperature) 49 ! - Check reaction of plants at wilting point and at atmospheric saturation 50 ! - Consider partial absorption of the net shortwave radiation by the skin layer 51 ! - Allow for water surfaces, check performance for bare soils 39 52 !------------------------------------------------------------------------------! 40 53 USE arrays_3d, & … … 42 55 43 56 USE cloud_parameters, & 44 ONLY: cp, l_d_r, l_v, rho_l, r_d, r_v57 ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v 45 58 46 59 USE control_parameters, & 47 60 ONLY: dt_3d, humidity, intermediate_timestep_count, & 48 intermediate_timestep_count_max, p t_surface, rho_surface,&49 surface_pressure, timestep_scheme, tsc61 intermediate_timestep_count_max, precipitation, pt_surface, & 62 rho_surface, surface_pressure, timestep_scheme, tsc 50 63 51 64 USE indices, & … … 119 132 root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers 120 133 soil_level = (/0.07_wp, 0.28_wp, 1.00_wp, 2.89_wp/), & !: soil layer depths (m) 121 soil_moisture = 0.0_wp!: soil moisture content (m3/m3)134 soil_moisture = 0.0_wp !: soil moisture content (m3/m3) 122 135 123 136 REAL(wp), DIMENSION(0:soil_layers) :: & … … 532 545 G = 0.0_wp 533 546 H = rho_cp * shf 534 LE = rho_l * l_v * qsws 547 548 IF ( humidity ) THEN 549 LE = rho_l * l_v * qsws 550 ELSE 551 LE = 0.0_wp 552 ENDIF 553 535 554 LE_veg = 0.0_wp 536 555 LE_soil = LE … … 658 677 DO j = nysg, nyng 659 678 k = nzb_s_inner(j,i) 679 ! 680 !-- Assure that r_a cannot be zero at model start 681 IF ( pt(k+1,j,i) == pt(k,j,i) ) pt(k+1,j,i) = pt(k+1,j,i) + 1.0E-10_wp 682 660 683 us(j,i) = 0.1_wp 661 684 ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i) … … 717 740 DO i = nxlg, nxrg 718 741 DO j = nysg, nyng 719 742 k = nzb_s_inner(j,i) 720 743 721 744 ! … … 726 749 lambda_skin = lambda_skin_u(j,i) 727 750 ENDIF 728 ! 729 !-- First step: calculate aerodyamic resistance. As pt(0), us, ts 730 !-- are not available for the current time step, data from the last 731 !-- time step is used here. 732 k = nzb_s_inner(j,i) 733 734 ! r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20) 735 r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / - (shf(j,i) + 1.0E-20) 751 752 ! 753 !-- First step: calculate aerodyamic resistance. As pt, us, ts 754 !-- are not available for the prognostic time step, data from the last 755 !-- time step is used here. Note that this formulation is the 756 !-- equivalent to the ECMWF formulation using drag coefficients 757 r_a(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20) 736 758 737 759 ! … … 793 815 794 816 q_s = 0.622_wp * e_s / surface_pressure 795 IF ( q_s .LE. q_p(k+1,j,i)) THEN 796 ! PRINT*, "dew fall at (before)", time_since_reference_point 797 r_canopy(j,i) = 0.0_wp 798 r_soil(j,i) = 0.0_wp 817 818 ! 819 !-- In case of dew fall, set resistances to zero. 820 !-- To do: what does that physically reasoning behind this? 821 IF ( humidity ) THEN 822 IF ( q_s .LE. q_p(k+1,j,i) ) THEN 823 r_canopy(j,i) = 0.0_wp 824 r_soil(j,i) = 0.0_wp 825 ENDIF 799 826 ENDIF 800 827 … … 808 835 809 836 810 ! Plant cannot transpirate below wilting point. here, r_canopy 811 ! should go to infinity 812 ! IF ( m_soil(k,j,i) .LT. m_wilt(j,i) ) THEN 813 ! f_LE_veg(j,i) = 0.0 814 ! ENDIF 837 ! 838 !-- If soil moisture is below wilting point, plants do no longer 839 !-- transpirate. 840 IF ( m_soil(k,j,i) .LT. m_wilt(j,i) ) THEN 841 f_LE_veg = 0.0 842 ENDIF 815 843 816 844 f_H = rho_cp / r_a(j,i) … … 831 859 Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4 832 860 833 ! 834 !-- Numerator of the prognostic equation 835 coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H / exn & 836 * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT * T_0(j,i) & 837 ) + lambda_skin * T_soil(0,j,i) 838 839 ! 840 !-- Denominator of the prognostic equation 841 coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT + & 842 lambda_skin + f_H / exn 861 IF ( humidity ) THEN 862 863 ! 864 !-- Numerator of the prognostic equation 865 coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H & 866 / exn * T_1 + f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT & 867 * T_0(j,i) ) + lambda_skin * T_soil(0,j,i) 868 869 ! 870 !-- Denominator of the prognostic equation 871 coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 + f_LE * dq_s_dT & 872 + lambda_skin + f_H / exn 873 874 ELSE 875 876 ! 877 !-- Numerator of the prognostic equation 878 coef_1 = Rn(j,i) + 3.0_wp * sigma_SB * T_0(j,i) ** 4 + f_H / & 879 exn * T_1 + lambda_skin * T_soil(0,j,i) 880 881 ! 882 !-- Denominator of the prognostic equation 883 coef_2 = 4.0_wp * sigma_SB * T_0(j,i) ** 3 & 884 + lambda_skin + f_H / exn 885 886 ENDIF 843 887 844 888 tend = 0.0_wp … … 877 921 G(j,i) = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i)) 878 922 H(j,i) = - f_H * ( pt_p(k+1,j,i) - pt_p(k,j,i) ) 879 LE(j,i) = - f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT * & 880 T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 881 882 LE_veg(j,i) = - f_LE_veg * ( q_p(k+1,j,i) - q_s + dq_s_dT * & 883 T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 884 LE_soil(j,i) = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT * & 885 T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 886 LE_liq(j,i) = - f_LE_liq * ( q_p(k+1,j,i) - q_s + dq_s_dT * & 887 T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 888 923 924 IF ( humidity ) THEN 925 LE(j,i) = - f_LE * ( q_p(k+1,j,i) - q_s + dq_s_dT & 926 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 927 928 LE_veg(j,i) = - f_LE_veg * ( q_p(k+1,j,i) - q_s + dq_s_dT & 929 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 930 LE_soil(j,i) = - f_LE_soil * ( q_p(k+1,j,i) - q_s + dq_s_dT & 931 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 932 LE_liq(j,i) = - f_LE_liq * ( q_p(k+1,j,i) - q_s + dq_s_dT & 933 * T_0(j,i) - dq_s_dT * T_0_p(j,i) ) 934 ENDIF 889 935 890 936 ! IF ( i == 1 .AND. j == 1 ) THEN … … 898 944 ! ENDIF 899 945 900 946 ! 947 !-- Calculate the true surface resistance 901 948 IF ( LE(j,i) .EQ. 0.0 ) THEN 902 ! PRINT*, "+++ Evapotranspiration -> 0"903 949 r_s(j,i) = 1.0E10 904 950 ELSE … … 908 954 909 955 ! 910 !-- Calculate change in liquid water reservoir due to dew fall or911 !-- evaporation of liquid water (to do: add interception from rainfall)912 IF ( q_s .LE. q_p(k+1,j,i)) THEN913 !914 !-- Check if reservoir is full (avoid values > m_liq_max)915 !-- In that case, LE_liq goes to LE_soil. In this case916 !-- LE_veg is zero anyway (because c_liq = 1), so that tend is917 !-- zero and no further check is needed918 IF ( m_liq(j,i) .EQ. m_liq_max ) THEN919 LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)920 LE_liq(j,i) = 0.0_wp921 ENDIF922 923 !924 !-- In case LE_veg becomes negative (unphysical behavior), let925 !-- the water enter the liquid water reservoir as dew on the926 !-- plant927 IF ( LE_veg(j,i) .LT. 0.0_wp ) THEN928 LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)929 LE_veg(j,i) = 0.0_wp930 ENDIF931 ENDIF932 933 tend = - LE_liq(j,i) * drho_l_lv934 935 936 m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend &937 + tsc(3) * tm_liq_m(j,i) )938 939 !940 !-- Check if reservoir is overfull -> reduce to maximum941 !-- (conservation of water is violated here)942 m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max)943 944 !945 !-- Check if reservoir is empty (avoid values < 0.0)946 !-- (conservation of water is violated here)947 m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp)948 949 950 !951 !-- Calculate m_liq tendencies for the next Runge-Kutta step952 IF ( timestep_scheme(1:5) == 'runge' ) THEN953 IF ( intermediate_timestep_count == 1 ) THEN954 tm_liq_m(j,i) = tend955 ELSEIF ( intermediate_timestep_count < &956 intermediate_timestep_count_max ) THEN957 tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tm_liq_m(j,i)958 ENDIF959 ENDIF960 961 !962 956 !-- Calculate fluxes in the atmosphere 963 957 shf(j,i) = H(j,i) / rho_cp 964 qsws(j,i) = LE(j,i) / rho_lv 965 966 ENDDO 958 959 ! 960 !-- Calculate change in liquid water reservoir due to dew fall or 961 !-- evaporation of liquid water 962 IF ( humidity ) THEN 963 ! 964 !-- If precipitation is activated, add rain water to LE_liq. 965 !-- precipitation_rate is given in mm. 966 IF ( precipitation ) THEN 967 LE_liq(j,i) = LE_liq(j,i) + precipitation_rate(j,i) & 968 * 0.001_wp * rho_l * l_v 969 ENDIF 970 ! 971 !-- If the air is saturated, check the reservoir water level 972 IF ( q_s .LE. q_p(k+1,j,i)) THEN 973 ! 974 !-- Check if reservoir is full (avoid values > m_liq_max) 975 !-- In that case, LE_liq goes to LE_soil. In this case 976 !-- LE_veg is zero anyway (because c_liq = 1), so that tend is 977 !-- zero and no further check is needed 978 IF ( m_liq(j,i) .EQ. m_liq_max ) THEN 979 LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i) 980 LE_liq(j,i) = 0.0_wp 981 ENDIF 982 983 ! 984 !-- In case LE_veg becomes negative (unphysical behavior), let 985 !-- the water enter the liquid water reservoir as dew on the 986 !-- plant 987 IF ( LE_veg(j,i) .LT. 0.0_wp ) THEN 988 LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i) 989 LE_veg(j,i) = 0.0_wp 990 ENDIF 991 ENDIF 992 993 tend = - LE_liq(j,i) * drho_l_lv 994 995 m_liq_p(j,i) = m_liq(j,i) + dt_3d * ( tsc(2) * tend & 996 + tsc(3) * tm_liq_m(j,i) ) 997 998 ! 999 !-- Check if reservoir is overfull -> reduce to maximum 1000 !-- (conservation of water is violated here) 1001 m_liq_p(j,i) = MIN(m_liq_p(j,i),m_liq_max) 1002 1003 ! 1004 !-- Check if reservoir is empty (avoid values < 0.0) 1005 !-- (conservation of water is violated here) 1006 m_liq_p(j,i) = MAX(m_liq_p(j,i),0.0_wp) 1007 1008 1009 ! 1010 !-- Calculate m_liq tendencies for the next Runge-Kutta step 1011 IF ( timestep_scheme(1:5) == 'runge' ) THEN 1012 IF ( intermediate_timestep_count == 1 ) THEN 1013 tm_liq_m(j,i) = tend 1014 ELSEIF ( intermediate_timestep_count < & 1015 intermediate_timestep_count_max ) THEN 1016 tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp & 1017 * tm_liq_m(j,i) 1018 ENDIF 1019 ENDIF 1020 1021 ! 1022 !-- Calculate moisture flux in the atmosphere 1023 qsws(j,i) = LE(j,i) / rho_lv 1024 1025 ENDIF 1026 967 1027 ENDDO 1028 ENDDO 968 1029 969 1030 … … 1221 1282 DO i = nxlg, nxrg 1222 1283 DO j = nysg, nyng 1223 1224 1284 k = nzb_s_inner(j,i) 1225 1285 !
Note: See TracChangeset
for help on using the changeset viewer.
