Changeset 2881


Ignore:
Timestamp:
Mar 13, 2018 4:24:40 PM (7 years ago)
Author:
maronga
Message:

bugfix in land surface model and new option in spinup

Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r2805 r2881  
    2525! -----------------
    2626! $Id$
     27! Bugfix: wrong loop structure for soil moisture calculation
     28!
     29! 2805 2018-02-14 17:00:09Z suehring
    2730! Bugfix in initialization of roughness over water surfaces
    2831!
     
    47794782             ENDDO
    47804783
    4781           ENDIF
    4782 
    4783        ENDDO
    4784 
    4785 
    4786        DO  m = 1, surf%ns
    4787 
    4788           IF (  .NOT.  surf%water_surface(m)  .AND.  calc_soil_moisture )  THEN
    4789 
    4790 
    4791 !
    4792 !--          Prognostic equation for soil moisture content. Only performed,
    4793 !--          when humidity is enabled in the atmosphere.
    4794              IF ( humidity )  THEN
    4795 !
    4796 !--             Calculate soil diffusivity (lambda_w) at the _layer level
    4797 !--             using linear interpolation. To do: replace this with
    4798 !--             ECMWF-IFS Eq. 8.81
    4799                 DO  k = nzb_soil, nzt_soil-1
    4800                    surf%lambda_w(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
    4801                                         * 0.5_wp
    4802                    surf%gamma_w(k,m)  = ( gamma_temp(k+1)  +  gamma_temp(k) )  &
    4803                                         * 0.5_wp
    4804                 ENDDO
    4805 !
    4806 !
    4807 !--             In case of a closed bottom (= water content is conserved),
    4808 !--             set hydraulic conductivity to zero to that no water will be
    4809 !--             lost in the bottom layer. As gamma_w is always a positive value,
    4810 !--             it cannot be set to zero in case of purely dry soil since this
    4811 !--             would cause accumulation of (non-existing) water in the lowest
    4812 !--             soil layer
    4813                 IF ( conserve_water_content .AND.                              &
    4814                      surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp )  THEN
    4815 
    4816                    surf%gamma_w(nzt_soil,m) = 0.0_wp
    4817                 ELSE
    4818                    surf%gamma_w(nzt_soil,m) = gamma_temp(nzt_soil)
    4819                 ENDIF     
    4820 
    4821 !--             The root extraction (= root_extr * qsws_veg / (rho_l     
    4822 !--             * l_v)) ensures the mass conservation for water. The         
    4823 !--             transpiration of plants equals the cumulative withdrawals by
    4824 !--             the roots in the soil. The scheme takes into account the
    4825 !--             availability of water in the soil layers as well as the root
    4826 !--             fraction in the respective layer. Layer with moisture below
    4827 !--             wilting point will not contribute, which reflects the
    4828 !--             preference of plants to take water from moister layers.
    4829 !
    4830 !--             Calculate the root extraction (ECMWF 7.69, the sum of
    4831 !--             root_extr = 1). The energy balance solver guarantees a
    4832 !--             positive transpiration, so that there is no need for an
    4833 !--             additional check.
    4834                 m_total = 0.0_wp
    4835                 DO  k = nzb_soil, nzt_soil
    4836                     IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
    4837                        m_total = m_total + surf%root_fr(k,m)                   &
    4838                                          * surf_m_soil%var_2d(k,m)
    4839                     ENDIF
    4840                 ENDDO 
    4841                  IF ( m_total > 0.0_wp )  THEN
     4784
     4785             IF (  calc_soil_moisture )  THEN
     4786
     4787!
     4788!--             Prognostic equation for soil moisture content. Only performed,
     4789!--             when humidity is enabled in the atmosphere.
     4790                IF ( humidity )  THEN
     4791!
     4792!--                Calculate soil diffusivity (lambda_w) at the _layer level
     4793!--                using linear interpolation. To do: replace this with
     4794!--                ECMWF-IFS Eq. 8.81
     4795                   DO  k = nzb_soil, nzt_soil-1
     4796               
     4797                      surf%lambda_w(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
     4798                                           * 0.5_wp
     4799                      surf%gamma_w(k,m)  = ( gamma_temp(k+1)  +  gamma_temp(k) )  &
     4800                                           * 0.5_wp
     4801                                           
     4802                   ENDDO
     4803!
     4804!
     4805!--                In case of a closed bottom (= water content is conserved),
     4806!--                set hydraulic conductivity to zero to that no water will be
     4807!--                lost in the bottom layer. As gamma_w is always a positive value,
     4808!--                it cannot be set to zero in case of purely dry soil since this
     4809!--                would cause accumulation of (non-existing) water in the lowest
     4810!--                soil layer
     4811                   IF ( conserve_water_content .AND.                           &
     4812                        surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp )  THEN
     4813
     4814                      surf%gamma_w(nzt_soil,m) = 0.0_wp
     4815                   ELSE
     4816                      surf%gamma_w(nzt_soil,m) = gamma_temp(nzt_soil)
     4817                   ENDIF     
     4818
     4819!--                The root extraction (= root_extr * qsws_veg / (rho_l     
     4820!--                * l_v)) ensures the mass conservation for water. The         
     4821!--                transpiration of plants equals the cumulative withdrawals by
     4822!--                the roots in the soil. The scheme takes into account the
     4823!--                availability of water in the soil layers as well as the root
     4824!--                fraction in the respective layer. Layer with moisture below
     4825!--                wilting point will not contribute, which reflects the
     4826!--                preference of plants to take water from moister layers.
     4827!
     4828!--                Calculate the root extraction (ECMWF 7.69, the sum of
     4829!--                root_extr = 1). The energy balance solver guarantees a
     4830!--                positive transpiration, so that there is no need for an
     4831!--                additional check.
     4832                   m_total = 0.0_wp
    48424833                   DO  k = nzb_soil, nzt_soil
    48434834                      IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
    4844                          root_extr(k) = surf%root_fr(k,m)                      &
    4845                                       * surf_m_soil%var_2d(k,m) / m_total
    4846                       ELSE
    4847                          root_extr(k) = 0.0_wp
     4835                         m_total = m_total + surf%root_fr(k,m)                 &
     4836                                * surf_m_soil%var_2d(k,m)
    48484837                      ENDIF
    4849                    ENDDO
    4850                 ENDIF
    4851 !
    4852 !--             Prognostic equation for soil water content m_soil_h.
    4853                 tend(:) = 0.0_wp
    4854 
    4855                 tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (             &
     4838                   ENDDO 
     4839                   IF ( m_total > 0.0_wp )  THEN
     4840                      DO  k = nzb_soil, nzt_soil
     4841                         IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
     4842                            root_extr(k) = surf%root_fr(k,m)                   &
     4843                                           * surf_m_soil%var_2d(k,m) / m_total
     4844                         ELSE
     4845                            root_extr(k) = 0.0_wp
     4846                         ENDIF
     4847                      ENDDO
     4848                   ENDIF
     4849!
     4850!--                Prognostic equation for soil water content m_soil_h.
     4851                   tend(:) = 0.0_wp
     4852
     4853                   tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (          &
    48564854                         surf_m_soil%var_2d(nzb_soil+1,m)                      &
    48574855                         - surf_m_soil%var_2d(nzb_soil,m) )                    &
     
    48614859                            * ddz_soil(nzb_soil)
    48624860
    4863 
    4864                 DO  k = nzb_soil+1, nzt_soil-1
    4865                    tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)  &
     4861                   DO  k = nzb_soil+1, nzt_soil-1
     4862                      tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)  &
    48664863                             - surf_m_soil%var_2d(k,m) ) * ddz_soil_center(k)    &
    48674864                             - surf%gamma_w(k,m)                                 &
     
    48714868                             * surf%qsws_veg(m) * drho_l_lv)                     &
    48724869                             ) * ddz_soil(k)
    4873                 ENDDO
    4874                 tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)                  &
     4870                   ENDDO
     4871                   tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)               &
    48754872                                   - surf%lambda_w(nzt_soil-1,m)               &
    48764873                                   * ( surf_m_soil%var_2d(nzt_soil,m)          &
     
    48824879                                  ) * ddz_soil(nzt_soil)             
    48834880
    4884                 surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) =                    &
     4881                   surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) =                 &
    48854882                                       surf_m_soil%var_2d(nzb_soil:nzt_soil,m) &
    48864883                                         + dt_3d * ( tsc(2) * tend(:)          &
     
    48884885   
    48894886!
    4890 !--             Account for dry soils (find a better solution here!)
    4891                 DO  k = nzb_soil, nzt_soil
    4892                    IF ( surf_m_soil_p%var_2d(k,m) < 0.0_wp )  surf_m_soil_p%var_2d(k,m) = 0.0_wp
    4893                 ENDDO
    4894 
    4895 !
    4896 !--             Calculate m_soil tendencies for the next Runge-Kutta step
    4897                 IF ( timestep_scheme(1:5) == 'runge' )  THEN
    4898                    IF ( intermediate_timestep_count == 1 )  THEN
    4899                       DO  k = nzb_soil, nzt_soil
    4900                          surf_tm_soil_m%var_2d(k,m) = tend(k)
    4901                       ENDDO
    4902                    ELSEIF ( intermediate_timestep_count <                      &
    4903                             intermediate_timestep_count_max )  THEN
    4904                       DO  k = nzb_soil, nzt_soil
    4905                          surf_tm_soil_m%var_2d(k,m) = -9.5625_wp * tend(k)     &
     4887!--                Account for dry soils (find a better solution here!)
     4888                   DO  k = nzb_soil, nzt_soil
     4889                      IF ( surf_m_soil_p%var_2d(k,m) < 0.0_wp )  surf_m_soil_p%var_2d(k,m) = 0.0_wp
     4890                   ENDDO
     4891 
     4892!
     4893!--                Calculate m_soil tendencies for the next Runge-Kutta step
     4894                   IF ( timestep_scheme(1:5) == 'runge' )  THEN
     4895                      IF ( intermediate_timestep_count == 1 )  THEN
     4896                         DO  k = nzb_soil, nzt_soil
     4897                            surf_tm_soil_m%var_2d(k,m) = tend(k)
     4898                         ENDDO
     4899                      ELSEIF ( intermediate_timestep_count <                   &
     4900                               intermediate_timestep_count_max )  THEN
     4901                         DO  k = nzb_soil, nzt_soil
     4902                            surf_tm_soil_m%var_2d(k,m) = -9.5625_wp * tend(k)  &
    49064903                                                    + 5.3125_wp                &
    49074904                                                    * surf_tm_soil_m%var_2d(k,m)
    4908                       ENDDO
    4909 
     4905                         ENDDO
     4906
     4907                      ENDIF
     4908                     
    49104909                   ENDIF
     4910                   
    49114911                ENDIF
     4912
    49124913             ENDIF
    49134914
  • palm/trunk/SOURCE/modules.f90

    r2797 r2881  
    2525! -----------------
    2626! $Id$
     27! Added flag for switching on/off calculation of soil moisture
     28!
     29! 2797 2018-02-08 13:24:35Z suehring
    2730! +ghf_av
    2831!
     
    12141217    LOGICAL ::  bc_ns_dirrad = .FALSE.                           !< north-south boundary condition dirichlet/radiation?
    12151218    LOGICAL ::  bc_ns_raddir = .FALSE.                           !< north-south boundary condition radiation/dirichlet?
     1219    LOGICAL ::  calc_soil_moisture_during_spinup = .FALSE.       !< namelist parameter
    12161220    LOGICAL ::  call_microphysics_at_all_substeps = .FALSE.      !< namelist parameter
    12171221    LOGICAL ::  call_psolver_at_all_substeps = .TRUE.            !< namelist parameter
  • palm/trunk/SOURCE/parin.f90

    r2849 r2881  
    2525! -----------------
    2626! $Id$
     27! Added flag for switching on/off calculation of soil moisture
     28!
     29! 2849 2018-03-05 10:49:33Z Giersch
    2730! Position of d3par namelist in parameter file is unimportant now
    2831!
     
    443446             bottom_salinityflux, building_height, building_length_x,          &
    444447             building_length_y, building_wall_left, building_wall_south,       &
     448             calc_soil_moisture_during_spinup,                                 &
    445449             call_psolver_at_all_substeps, call_microphysics_at_all_substeps,  &
    446450             canyon_height,                                                    &
  • palm/trunk/SOURCE/time_integration_spinup.f90

    r2818 r2881  
    2525! -----------------
    2626! $Id$
     27! Added flag for switching on/off calculation of soil moisture
     28!
     29! 2818 2018-02-19 16:42:36Z maronga
    2730! Velocity components near walls/ground are now set to the profiles stored in
    2831! u_init and v_init. Activated soil moisture calculation during spinup.
     
    8386
    8487    USE control_parameters,                                                    &
    85         ONLY:  averaging_interval_pr, constant_diffusion, constant_flux_layer, &
     88        ONLY:  averaging_interval_pr, calc_soil_moisture_during_spinup,        &
     89               constant_diffusion, constant_flux_layer,                        &
    8690               coupling_start_time, current_timestep_number,                   &
    8791               data_output_during_spinup, disturbance_created, dopr_n, do_sum, &
     
    348352!--             Call for horizontal upward-facing surfaces
    349353                CALL lsm_energy_balance( .TRUE., -1 )
    350                 CALL lsm_soil_model( .TRUE., -1, .TRUE. )
     354                CALL lsm_soil_model( .TRUE., -1, calc_soil_moisture_during_spinup )
    351355!
    352356!--             Call for northward-facing surfaces
    353357                CALL lsm_energy_balance( .FALSE., 0 )
    354                 CALL lsm_soil_model( .FALSE., 0, .TRUE. )
     358                CALL lsm_soil_model( .FALSE., 0, calc_soil_moisture_during_spinup )
    355359!
    356360!--             Call for southward-facing surfaces
    357361                CALL lsm_energy_balance( .FALSE., 1 )
    358                 CALL lsm_soil_model( .FALSE., 1, .TRUE. )
     362                CALL lsm_soil_model( .FALSE., 1, calc_soil_moisture_during_spinup )
    359363!
    360364!--             Call for eastward-facing surfaces
    361365                CALL lsm_energy_balance( .FALSE., 2 )
    362                 CALL lsm_soil_model( .FALSE., 2, .TRUE. )
     366                CALL lsm_soil_model( .FALSE., 2, calc_soil_moisture_during_spinup )
    363367!
    364368!--             Call for westward-facing surfaces
    365369                CALL lsm_energy_balance( .FALSE., 3 )
    366                 CALL lsm_soil_model( .FALSE., 3, .TRUE. )
     370                CALL lsm_soil_model( .FALSE., 3, calc_soil_moisture_during_spinup )
    367371
    368372                CALL cpu_log( log_point(54), 'land_surface', 'stop' )
Note: See TracChangeset for help on using the changeset viewer.