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

bugfix in land surface model and new option in spinup

File:
1 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
Note: See TracChangeset for help on using the changeset viewer.