Changeset 1500


Ignore:
Timestamp:
Dec 3, 2014 5:42:41 PM (7 years ago)
Author:
maronga
Message:

bugfixes and adjustments in land surface model

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/Makefile

    r1497 r1500  
    2020# Current revisions:
    2121# ------------------
    22 #
     22# Bugfix: missing adjustments for land surface model and radiation model
    2323#
    2424# Former revisions:
     
    298298check_open.o: modules.o mod_kinds.o mod_particle_attributes.o
    299299check_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
    301301close_file.o: modules.o mod_kinds.o
    302302compute_vpt.o: modules.o mod_kinds.o
     
    398398poisfft.o: modules.o cpulog.o fft_xy.o mod_kinds.o tridia_solver.o
    399399poismg.o: modules.o cpulog.o mod_kinds.o
    400 prandtl_fluxes.o: modules.o mod_kinds.o
     400prandtl_fluxes.o: modules.o mod_kinds.o land_surface_model.o
    401401pres.o: modules.o cpulog.o mod_kinds.o poisfft.o
    402402print_1d.o: modules.o cpulog.o mod_kinds.o
     
    422422sum_up_3d_data.o: modules.o cpulog.o mod_kinds.o
    423423surface_coupler.o: modules.o cpulog.o mod_kinds.o
    424 swap_timelevel.o: modules.o cpulog.o mod_kinds.o
     424swap_timelevel.o: modules.o cpulog.o mod_kinds.o land_surface_model.o
    425425temperton_fft.o: modules.o mod_kinds.o
    426426time_integration.o: modules.o advec_ws.o buoyancy.o calc_mean_profile.o \
  • palm/trunk/SOURCE/check_parameters.f90

    r1497 r1500  
    958958       ENDIF   
    959959
    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
    968970
    969971       IF ( .NOT. prandtl_layer )  THEN
  • palm/trunk/SOURCE/land_surface_model.f90

    r1497 r1500  
    2020! Current revisions:
    2121! -----------------
    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.
    2325!
    2426! Former revisions:
     
    3739! H-TESSEL. The implementation is based on the formulation implemented in the
    3840! 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
    3952!------------------------------------------------------------------------------!
    4053     USE arrays_3d,                                                            &
     
    4255
    4356     USE cloud_parameters,                                                     &
    44          ONLY: cp, l_d_r, l_v, rho_l, r_d, r_v
     57         ONLY: cp, l_d_r, l_v, precipitation_rate, rho_l, r_d, r_v
    4558
    4659     USE control_parameters,                                                   &
    4760         ONLY: dt_3d, humidity, intermediate_timestep_count,                   &
    48                intermediate_timestep_count_max, pt_surface, rho_surface,       &
    49                surface_pressure, timestep_scheme, tsc
     61               intermediate_timestep_count_max, precipitation, pt_surface,     &
     62               rho_surface, surface_pressure, timestep_scheme, tsc
    5063
    5164     USE indices,                                                              &
     
    119132              root_fraction = (/0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp/), & !: distribution of root surface area to the individual soil layers
    120133              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)
    122135
    123136    REAL(wp), DIMENSION(0:soil_layers) ::   &
     
    532545       G       = 0.0_wp
    533546       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
    535554       LE_veg  = 0.0_wp
    536555       LE_soil = LE
     
    658677          DO  j = nysg, nyng
    659678             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
    660683             us(j,i) = 0.1_wp
    661684             ts(j,i) = (pt(k+1,j,i) - pt(k,j,i)) / r_a(j,i)
     
    717740       DO i = nxlg, nxrg
    718741          DO j = nysg, nyng
    719 
     742             k = nzb_s_inner(j,i)
    720743
    721744!
     
    726749                lambda_skin = lambda_skin_u(j,i)
    727750             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)
    736758
    737759!
     
    793815
    794816             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
    799826             ENDIF
    800827
     
    808835
    809836
    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
    815843
    816844             f_H  = rho_cp / r_a(j,i)
     
    831859             Rn(j,i) = Rn(j,i) + sigma_SB * T_0(j,i) ** 4
    832860
    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
    843887
    844888             tend = 0.0_wp
     
    877921             G(j,i)         = lambda_skin * (T_0_p(j,i) - T_soil(0,j,i))
    878922             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
    889935
    890936!              IF ( i == 1 .AND. j == 1 )  THEN
     
    898944!              ENDIF
    899945
    900 
     946!
     947!--          Calculate the true surface resistance
    901948             IF ( LE(j,i) .EQ. 0.0 )  THEN
    902 !               PRINT*, "+++ Evapotranspiration -> 0"
    903949                r_s(j,i) = 1.0E10
    904950             ELSE
     
    908954
    909955!
    910 !--          Calculate change in liquid water reservoir due to dew fall or
    911 !--          evaporation of liquid water (to do: add interception from rainfall)
    912              IF ( q_s .LE. q_p(k+1,j,i))  THEN
    913 !
    914 !--             Check if reservoir is full (avoid values > m_liq_max)
    915 !--             In that case, LE_liq goes to LE_soil. In this case
    916 !--             LE_veg is zero anyway (because c_liq = 1), so that tend is
    917 !--             zero and no further check is needed
    918                 IF ( m_liq(j,i) .EQ. m_liq_max )  THEN
    919                    LE_soil(j,i) = LE_soil(j,i) + LE_liq(j,i)
    920                    LE_liq(j,i) = 0.0_wp
    921                 ENDIF
    922 
    923 !
    924 !--             In case LE_veg becomes negative (unphysical behavior), let
    925 !--             the water enter the liquid water reservoir as dew on the
    926 !--             plant
    927                 IF ( LE_veg(j,i) .LT. 0.0_wp )  THEN
    928                    LE_liq(j,i) = LE_liq(j,i) + LE_veg(j,i)
    929                    LE_veg(j,i) = 0.0_wp
    930                 ENDIF
    931              ENDIF                   
    932  
    933              tend = - LE_liq(j,i) * drho_l_lv
    934  
    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 maximum
    941 !--          (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 step
    952              IF ( timestep_scheme(1:5) == 'runge' )  THEN
    953                 IF ( intermediate_timestep_count == 1 )  THEN
    954                    tm_liq_m(j,i) = tend
    955                 ELSEIF ( intermediate_timestep_count <                         &
    956                          intermediate_timestep_count_max )  THEN
    957                    tm_liq_m(j,i) = -9.5625_wp * tend + 5.3125_wp * tm_liq_m(j,i)
    958                 ENDIF
    959              ENDIF
    960 
    961 !
    962956!--          Calculate fluxes in the atmosphere
    963957             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
    9671027          ENDDO
     1028       ENDDO
    9681029
    9691030
     
    12211282       DO i = nxlg, nxrg   
    12221283          DO j = nysg, nyng
    1223 
    12241284             k = nzb_s_inner(j,i)
    12251285!
Note: See TracChangeset for help on using the changeset viewer.