water.f90

CLM集合卡曼滤波数据同化算法

  SUBROUTINE water ( ivt     ,lb      ,iub    ,dtime  ,&
             z      ,dz      ,zi      ,bsw    ,porsl  ,&
             phi0   ,hksati  ,rootfr  ,rootr  ,tss    ,&
             wliq   ,wice    ,pg_rain ,sm     ,etr    ,&
             qseva  ,qsdew   ,qsubl   ,qfros  ,etrc   ,&
             rsur   ,rnof    ,fevpg   ,fevpa  ,lfevpa )

!=======================================================================
!      Source file: water.f
! Original version: Yongjiu Dai, September 15, 1999
!
! this is the main subroutine to execute the calculation of water processes
!
! (1) water flow within snow (see clm_snowwater.f90)
!
! (2) water flow within soi (see clm_soilwater.f90)
!
! (3) runoff
!   the original code was provide by Robert E. Dickinson based on following clues:
!   exponential decrease of Ksat added, a water table level determination
!   level added including highland and lowland levels and fractional area
!   of wetland (water table above the surface. Runoff is parametrized
!   from the lowlands in terms of precip incident on wet areas and a base
!   flow, where these are estimated using ideas from TOPMODEL.
!
!   the original scheme was modified by Z.-L. Yang and G.-Y. Niu,
!   o  using a new method to determine water table depth and
!      the fractional wet area (fcov)
!   o  computing runoff (surface and subsurface) from this
!      fraction and the remaining fraction (i.e. 1-fcov)
!   o  for the 1-fcov part, using BATS1e method to compute
!     surface and subsurface runoff.
!
!   Revision histroy:
!   the original code on soil moisture and runoff were provided by
!   R. E. Dickinson in July 1996
!   15 September 1999: Yongjiu Dai; Initial code
!   12 November 1999:  Z.-L. Yang and G.-Y. Niu
!=======================================================================

  USE PHYCON_MODULE ! physical constant
  IMPLICIT NONE
!    
!-------------------------- Dummy argument ------------------------------
!   
  integer, INTENT(in) :: &
        ivt,             &! land cover type
        lb,              &! lower bound of array
        iub               ! upper bound of array

  real, INTENT(in) :: &
        dtime,           &! time step (s)
        z (lb : iub),    &! layer depth (m)
        dz(lb : iub),    &! layer thickness (m)
        zi(lb - 1 : iub)  ! interface level below a "z" level (m)

! Input soil physical parameters
  real, INTENT(in) :: &
        bsw(1 : iub),    &! Clapp-Hornberger "B"
        porsl(1 : iub),  &! saturated volumetric soil water content(porosity)
        phi0(1 : iub),   &! saturated soil suction (mm)
        hksati(1 : iub), &! hydraulic conductivity at saturation (mm h2o/s)
        rootfr(1 : iub)   ! fraction of roots in a layer, 

! variables
  real, INTENT(in) :: &
        tss(lb : iub),   &! soil/snow skin temperature (K)
        pg_rain,         &! rainfall after removal of interception (mm h2o/s)
        sm,              &! snow melt (mm h2o/s)
        etr,             &! actual transpiration (mm h2o/s)
        qseva,           &!ground surface evaporation rate (mm h2o/s)
        qsdew,           &!ground surface dew formation (mm h2o /s) [+]
        qsubl,           &!sublimation rate from snow pack (mm h2o /s) [+]
        qfros             !surface dew added to snow pack (mm h2o /s) [+]

  real, INTENT(inout) :: &
        wice(lb : iub),  &! ice lens (kg/m2)
        wliq(lb : iub),  &! liquid water (kg/m2)
        rootr(1 : iub),  &! root resistance of a layer, all layers add to 1.0
        fevpg,           &! evaporation heat flux from ground [mm/s]
        fevpa,           &! evapotranspiration from canopy height to atmosphere [mm/s]
        lfevpa            ! latent heat flux from canopy height to atmosphere [W/m2]

  real, INTENT(out) :: &
        rsur,            &! surface runoff (mm h2o/s)
        rnof,            &! total runoff (mm h2o/s)
        etrc              ! maximum possible transpiration rate (mm h2o/s)
!                    
!----------------------- local variables ----------------------------
!                   
  integer                &!
       i,                &! loop counter
       j,                &! do loop/array indices
       jwatp              ! index for high water table level

  real hk(1 : iub),      &! hydraulic conductivity (mm h2o/s)
       dhkdw(1 : iub),   &! d(hk)/d(vol_liq)
       dwat(1 : iub),    &! change in soil water
       fcov,             &! fractional area with water table at surface
       gwat,             &! net water input from top
       qinfl,            &! infiltration rate (mm h2o/s)
       qout_snowb,       &! rate of water out of snow bottom (mm/s)
       roota,            &! accumulates root resistance factors
       rresis(1 : iub),  &! soil water contribution to root resistance
       rsubst,           &! subsurface runoff (mm h2o/s)
       s(1 : iub),       &! wetness of soil (including ice)
       s_node,           &! vol_liq/porosity
       smpmax,           &! wilting point potential in mm
       smp_node,         &! matrix potential
       vol_liq(1 : iub), &! partitial volume of liquid water in layer
       vol_ice(1 : iub), &! partitial volume of ice lens in layer
  eff_porosity(1 : iub), &! effective porosity = porosity - vol_ice
       wat1,             &! water in top soil layer
       wtsub,            &! weight factor for subsurface runoff 
       xs,               &! excess soil water above saturation
       zwice,            &! the sum of ice mass of soil (kg/m2)
       zwt,              &! water table depth
       zmm (1 : iub),    &! layer depth (mm)
       dzmm(1 : iub),    &! layer thickness (mm)
       zimm(0 : iub)      ! interface level below a "z" level (mm)

! for Z.-L. Yang & G.-Y. Niu's modification
  real zmean              ! The surface soil layers contributing to runoff
  real wmean              ! The averaged soil wetness in surface soil layers
  real fz                 ! coefficient for water table depth
                          ! 
  real zsat               ! hydraulic conductivity weighted soil thickness
  real wsat               ! hydraulic conductivity weighted soil wetness
  real rsubstw            ! subsurface runoff from "wet" part (mm h2o/s)
  real rsubstd            ! subsurface runoff from "dry" part (mm h2o/s)
  real dzksum             ! hydraulic conductivity weighted soil thickness
!
!=======================================================================
! [1] the change of snow mass and the snow water onto soil
!=======================================================================
!
      if (lb >=1)then
         gwat = pg_rain + sm - qseva
      else
         call snowwater ( lb, dtime, &
                          pg_rain, qseva, qsdew, qsubl, qfros, &
                          dz(lb), wice(lb), wliq(lb), qout_snowb )
         gwat = qout_snowb
      endif
!
!=======================================================================
! [2] Surface runoff
!=======================================================================
      if(ivt/=11 .AND. ivt/=15)then   ! NOT wetland and land ice
! Porosity of soil, partitial volume of ice and liquid
         zwice = 0.
         do i = 1, iub
            zwice = zwice + wice(i)
            vol_ice(i) = min(porsl(i), wice(i)/(dz(i)*dice))
            eff_porosity(i) = porsl(i)-vol_ice(i)
            vol_liq(i) = min(eff_porosity(i), wliq(i)/(dz(i)*rhowat))
            s(i) = min(1.,(vol_ice(i)+vol_liq(i))/porsl(i))
         enddo

! Determine water table 
! ---------------------------------   
         wmean = 0.   
         fz    = 1.0                    
         do i  = 1, iub                
            wmean = wmean + s(i)*dz(i)   
         enddo                       
                                 
         zwt = fz * (zi(iub) - wmean)  
                                   
! saturation fraction
         fcov = wtfact*min(1.,exp(-zwt))
        
! modified surface runoff according to the concept of TOPMODEL    
         wmean = 0.                                                
         zmean = 0.                                               
                                                               
         do i = 1, 3                                            
            zmean = zmean + dz(i)                                 
            wmean = wmean + s(i) * dz(i)                         
         enddo                                               
         wmean = wmean / zmean                              
                                                         
         rsur =  max(0.,     fcov *                 gwat) &  
              +  max(0., (1.-fcov)*min(1.,wmean**4)*gwat)   
        
! Add in hillslope runoff
! Infiltration into surface soil layer (minus the evaporation)
! and the derivative of evaporation with respect to water in top soil layer

         qinfl = gwat - rsur 

      else
        if(ivt==11)then
           rsur=0.;   qinfl=gwat  ! wetland
        endif
        if(ivt==15)then
           rsur=gwat; qinfl=0.    ! land ice
        endif
      endif
!
!=======================================================================
! [3] Set up r, a, b, and c vectors for tridiagonal solution and renew bl
!=======================================================================
      if(ivt/=11 .AND. ivt/=15)then   ! NOT wetland and land ice
! following length units are all in millimeters
         zimm(0) = zi(0)*1.e3
         do i = 1, iub
            zmm(i) = z(i)*1.e3
            dzmm(i) = dz(i)*1.e3
            zimm(i) = zi(i)*1.e3
         enddo

         call soilwater ( iub,dtime, &
                          porsl,phi0,bsw,hksati, &
                          zmm,dzmm,zimm,tss(1),vol_liq,vol_ice,eff_porosity, &
                          qinfl,etr,rootr,dwat,hk,dhkdw  )
!
! Renew the mass of liquid water
         do i= 1, iub
            wliq(i) = max(0.,wliq(i) + dwat(i)*dzmm(i))
!*02/02/2001
            if(dwat(i) > eff_porosity(i)-vol_liq(i)) &
               dwat(i) = eff_porosity(i)-vol_liq(i)
         enddo
      endif

!=======================================================================
! [4] Streamflow and total runoff
!=======================================================================
      if(ivt/=11 .AND. ivt/=15)then      ! NOT wetland and land ice
! The amount of streamflow is assumed maintained by flow from the 
! lowland water table with different levels contributing according to 
! their thickness and saturated hydraulic conductivity, i.e. a given 
! layer below the water table interface loses water at a rate per unit 
! depth given by rsubst*hk/(sum over all layers below this water table 
! of hk*dz). Because this is a slow smooth process, and not strongly 
! coupled to water in any one layer, it should remain stable for 
! explicit time differencing. Hence, for simplicity it is removed
! explicitly prior to the main soil water calculation.
! Another assumption: no subsurface runoff for ice mixed soil 

! Subsurface runoff
      rsubst = 0.
      rsubstw = 0.
      rsubstd = 0.
                                                             
      if(zwice <= 0.)then 
        zsat = 0.
        wsat = 0.
        dzksum = 0.
        do i = 6,iub-1
           zsat = zsat + dz(i)*hk(i)
           wsat = wsat + s(i)*dz(i)*hk(i)
           dzksum  = dzksum   + hk(i)*dz(i)
        end do 
        wsat = wsat / zsat 
                                                  
!*      rsubstd =  (1.-fcov)*4.e-2* wsat ** (2.*bsw(iub)+3.)  ! mm/s 
!*      rsubstw =  fcov * 1.e-5 * exp(-zwt)              ! mm/s
                                                                
        rsubst = rsubstd + rsubstw

!*      do i = 6, iub-1
!*         wliq(i) = wliq(i) - dtime * rsubst * dz(i)*hk(i) /dzksum
!*      enddo
      endif

! Determine water in excess of saturation and
! Add excess water back to underlying soil layers 
! Any left over put into runoff

        xs = max(0., wliq(1)-(pondmx+eff_porosity(1)*dzmm(1)))
        wliq(1) = min(pondmx+eff_porosity(1)*dzmm(1), wliq(1))

        do i = 2, iub
           wliq(i) = wliq(i) + xs
           xs = max(0., wliq(i)-eff_porosity(i)*dzmm(i))     ! [mm]
           wliq(i) = min(eff_porosity(i)*dzmm(i), wliq(i))
        end do
 
! Sub-surface runoff and drainage 
        rsubst = rsubst + xs/dtime &
               + hk(iub) + dhkdw(iub)*dwat(iub) ! [mm/s]
 
! Collect total runoff
        rnof  = rsubst + rsur                                     ! [mm/s]

      else 
        if(ivt==11)then
           rsubst=0.; rnof=0.          ! wetland
        endif
        if(ivt==15)then
           rsubst=0.; rnof=rsur        ! land ice
        endif
      endif

!=======================================================================
! [5] Implicit evaporation for application to soil temperature equation
!=======================================================================

      if(ivt/=11 .AND. ivt/=15)then             ! NOT wetland and lanc ice
! Renew the ice and liquid mass due to condensation
        if(lb >= 1)then
           wliq(1) = wliq(1) + qsdew*dtime
           wice(1) = wice(1) + (qfros-qsubl)*dtime
        endif
      endif
!     write(6,100) ((wliq(i)+wice(i))/dzmm(i),i=1,10)
100   format(10f6.2)

!*      do i = 1, iub
!*          if(zi(i-1) > 2.0)then
!*           wliq(i) = dzmm(i)*porsl(i)
!*           wice(i) = 0.
!*          endif
!*      enddo

!=======================================================================
! [6] Root resistance factors and maximum possible transpiration rate
!=======================================================================
! New potential and root resistance factors
      if(ivt/=11 .AND. ivt/=15)then             ! NOT wetland and lanc ice
        roota = 1.e-10                    ! must be non-zero to begin
        do i = 1, iub
          if(tss(i) > tfrz)then
             smpmax = -1.5e5
         	   s_node = max(wliq(i)/(dzmm(i)*eff_porosity(i)),0.01)
             smp_node = max(smpmax, -phi0(i)*s_node**(-bsw(i))) 
             rresis(i)=(1.-smp_node/smpmax)/(1.+phi0(i)/smpmax) 
             rootr(i) = rootfr(i)*rresis(i)
	     roota    = roota + rootr(i)
          else
             rootr(i) = 0.
          endif
        end do
!
! Normalize root resistances to get layer contribution to ET
        do i = 1, iub
           rootr(i) = rootr(i)/roota
        end do
!	
! Determine maximum possible transpiration rate 
        etrc = trsmx0*roota

!       write(66,200) (rootr(i),i=1,10), roota, etrc
200     format(12e12.4)
      else
        if(ivt==11)then
           rootr(1:iub)=(/(rootfr(i), i=1,iub)/); etrc=trsmx0
        endif
        if(ivt==15)then
           rootr(1:iub)=(/(0., i=1,iub)/); etrc=0.
        endif
      endif

  END SUBROUTINE water