thermalk.f90

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

  SUBROUTINE thermalk ( ivt  ,lb   ,iub   ,& 
                        tss  ,wice ,wliq  ,scv    ,&
                        csol ,dkmg ,dkdry ,dksatu ,porsl ,&
                        dz   ,z    ,zi    ,tk     ,cv )

! ======================================================================
!      Source file: thermalk.f90
! Original version: Yongjiu Dai, September 15, 1999
!
!calculation of thermal conductivities and heat capacities of snow / soil layers
!
!(1) the volumetric heat capacity is calculated as a linear combination
!    in terms of the volumetric fraction of the constituent phases.
!
!(2) The thermal conductivity of soil is computed from
!    the algorithm of Johansen (as reported by Farouki 1981), and of snow is from
!    the formulation used in SNTHERM (Jordan 1991).
!
!    The thermal conductivities at the interfaces between two neighbor layers
!    (j, j+1) are derived from an assumption that the flux across the interface
!    is equal to that from the node j to the interface and the flux from the
!    interface to the node j+1.
! ======================================================================

  USE PHYCON_MODULE ! physical constant
  IMPLICIT NONE

! input
  integer, INTENT(in) :: & 
        ivt,          &! land cover type
        lb,           &! lower bound of array
        iub            ! upper bound of array

  real, INTENT(in) :: &
        dz  (lb:iub), &! layer thickiness [m]
        z   (lb:iub), &! node depth [m]
        zi(lb-1:iub), &! interface depth [m]
        tss (lb:iub), &! Nodal temperature [K]
        wice(lb:iub), &! ice lens [kg/m2]
        wliq(lb:iub), &! liqui water [kg/m2]
        scv            ! snow water equivalent [mm]

  real, INTENT(in) :: &
        csol(1:iub),  &! heat capacity of soil soilds [J/(m3 K)]
        dkmg(1:iub),  &! dkm**dmvol, where dkm is the mineral conductivity and
        dkdry(1:iub), &! thermal conductivity for dry soil [W/m-K]
        dksatu(1:iub),&! Thermal conductivity of saturated soil [W/m-K]
        porsl(1:iub)   ! fractional volume between soil grains=1.-dmvol

! output
  real, INTENT(out) :: &
        cv(lb:iub),   &! heat capacity [J/(m2 K)]
        tk(lb:iub)     ! thermal conductivity [W/(m K)]

! local
  real  rhosnow,      &! partitial density of water (ice + liquid)
        dksat,        &! thermal conductivity for saturated soil (j/(k s m))
        dke,          &! kersten number
        fl,           &! fraction of liquid or unfrozen water to total water
        satw,         &! relative total water content of soil.
        thk(lb:iub)    ! thermal conductivity of layer

  integer i

!-----------------------------------------------------------------------
! Thermal conductivity of soil from Farouki (1981),
      do i = 1, iub

         if(ivt/=11 .AND. ivt/=15)then          ! NOT glacier and wetland
            satw = (wliq(i)/rhowat + wice(i)/dice) / (dz(i)*porsl(i))
            satw = min(1., satw)

            if(satw > .1e-6)then          

               fl = wliq(i)/(wice(i)+wliq(i))
   
               if(tss(i) >= tfrz) then       ! Unfrozen soil
                  dke = max(0., alog10(satw) + 1.0)
                  dksat = dksatu(i)
               else                          ! Frozen soil
                  dke = satw
                  dksat = dkmg(i)*0.249**(fl*porsl(i))*2.29**porsl(i)
               end if
               thk(i) = dke*dksat + (1.-dke)*dkdry(i)
            else    
               thk(i) = dkdry(i)
            end if
         else
            thk(i) = tkwat
            if(tss(i) < tfrz) thk(i) = tkice
         endif
      enddo

! Thermal conductivity of snow, which from Jordan (1991) pp. 18
      if(lb < 1)then
        do i = lb, 0
          rhosnow = (wice(i)+wliq(i))/dz(i)

!*        thk(i) = tkair+(7.75e-5*rhosnow+1.105e-6*rhosnow*rhosnow)*(tkice-tkair)

          thk(i) = 0.088+(7.75e-5*rhosnow+1.105e-6*rhosnow*rhosnow)*(tkice-0.088)
        enddo
      endif

! Thermal conductivity at the layer interface
      do i = lb, iub-1

! the following consideration is try to avoid the snow conductivity 
! to be dominant in the thermal conductivity of the interface. 
! Because when the distance of bottom snow node to the interfacee 
! is larger than that of interface to top soil node,
! the snow thermal conductivity will be dominant, and the result is that 
! lees heat tranfer between snow and soil 
         if((i==0) .AND. (z(i+1)-zi(i)<zi(i)-z(i)))then
            tk(i) = 2.*thk(i)*thk(i+1)/(thk(i)+thk(i+1))
            tk(i) = max(0.5*thk(i+1),tk(i))
         else
            tk(i) = thk(i)*thk(i+1)*(z(i+1)-z(i)) &
                  /(thk(i)*(z(i+1)-zi(i))+thk(i+1)*(zi(i)-z(i)))
         endif
      enddo
      tk(iub) = 0.
!
! ----------------------------------------------------------------------
! Heat capacity
!
! Soil heat capacity, which from de Vires (1963)

      do i = 1, iub
        if(ivt/=11 .AND. ivt/=15)then
           cv(i) = csol(i)*(1.-porsl(i))*dz(i) + (wice(i)*ci + wliq(i)*cl)
        else
           cv(i) = (wice(i)*ci + wliq(i)*cl)
        endif
      end do
      if(lb == 1 .AND. scv > 0.) cv(1) = cv(1) + ci*scv

! Snow heat capacity
      if(lb < 1)then
        do i = lb, 0
        cv(i) = cl*wliq(i) + ci*wice(i)
        enddo
      endif

  END SUBROUTINE thermalk