clmtci.f90

来自「CLM集合卡曼滤波数据同化算法」· F90 代码 · 共 358 行

  SUBROUTINE clmtci (kpt    ,msl    ,ivt    ,isc    ,sand   ,clay   ,&
                     csol   ,porsl  ,phi0   ,bsw    ,dkmg   ,dksatu ,&
                     dkdry  ,albsol ,hksati ,z0m    ,displa ,sqrtdi ,&
                     effcon ,vmax25 ,slti   ,hlti   ,shti   ,hhti   ,&                 
                     trda   ,trdm   ,trop   ,gradm  ,binter ,extkn  ,&
                     albvgs ,albvgl ,chil   ,ref    ,tran   ,rootfr ,&
                     zsoi   ,dzsoi  ,zsoih)
!=======================================================================
!      source file: clmtci.f90
! Original version: Yongjiu Dai, September 15, 1999
!
!=======================================================================

  USE CLM_TABLE_MODULE

  IMPLICIT NONE

!-----------------------------------------------------------------------
! Dummy argument

  integer, INTENT(in) :: &
        kpt,                 &!total number of clm land points, including subgrid points
        msl,                 &!number of model soil layers
        ivt(kpt),            &!index for vegetation type [-]
        isc(kpt)              !index for soil color type [-]

  real, INTENT(inout) :: &
        sand(kpt),           &! percent of snad
        clay(kpt)             ! percent of clay

! Soil parameters
  real, INTENT(out) :: &
        csol  (1:msl,1:kpt), &!heat capacity of soil solids [J/(m3 K)]
        porsl (1:msl,1:kpt), &!fraction of soil that is voids [-]
        phi0  (1:msl,1:kpt), &!minimum soil suction [mm]
        bsw   (1:msl,1:kpt), &!clapp and hornbereger "b" parameter [-]
        dkmg  (1:msl,1:kpt), &!thermal conductivity of soil minerals [W/m-K]
        dksatu(1:msl,1:kpt), &!thermal conductivity of saturated soil [W/m-K]
        dkdry (1:msl,1:kpt), &!thermal conductivity for dry soil  [W/(m-K)]
        hksati(1:msl,1:kpt), &!hydraulic conductivity at saturation [mm h2o/s] 
        albsol(kpt)           !soil albedo for different coloured soils [-]

! Vegetation static parameters
  real, INTENT(out) :: &
        z0m(kpt),            &!aerodynamic roughness length [m]
        displa(kpt),         &!displacement height [m]
        sqrtdi(kpt),         &!inverse sqrt of leaf dimension [m**-0.5]
        effcon(kpt),         &! quantum efficiency of RuBP regeneration (mol CO2 / mol quanta)
        vmax25(kpt),         &! maximum carboxylation rate at 25 C at canopy top
        slti(kpt),           &! slope of low temperature inhibition function      (0.2)
        hlti(kpt),           &! 1/2 point of low temperature inhibition function  (288.16)
        shti(kpt),           &! slope of high temperature inhibition function     (0.3)
        hhti(kpt),           &! 1/2 point of high temperature inhibition function (313.16)
        trda(kpt),           &! temperature coefficient in gs-a model             (1.3)
        trdm(kpt),           &! temperature coefficient in gs-a model             (328.16)
        trop(kpt),           &! temperature coefficient in gs-a model             (298.16)
        gradm(kpt),          &! conductance-photosynthesis slope parameter
        binter(kpt),         &! conductance-photosynthesis intercep
        extkn(kpt),          &!

        albvgs(kpt),         &!veg. albedo for wavelengths < 0.7 microns [-]
        albvgl(kpt),         &!veg. albedo for wavelengths > 0.7 microns [-]
        chil(kpt),           &! leaf angle distribution factor
        ref(2,2,kpt),        &! leaf reflectance (iw=iband, il=life and dead)
        tran(2,2,kpt),       &! leaf transmittance (iw=iband, il=life and dead)
        rootfr(1:msl,1:kpt)   !fraction of roots in each soil layer 

! Soil layer thickness, depths
  real, INTENT(out) :: &
        dzsoi(1:msl,1:kpt),  &!soil node thickness [m]
        zsoi(1:msl,1:kpt),   &!soil layer depth [m]
        zsoih(0:msl,1:kpt)    !interface level below a zsoi level [m]

! Local variables
  integer i, j, k             !loop indices
  real  scalez                !scale factor for vertical coordinate [m]
  real  bd                    !bulk density of dry soil material [kg/m^3]
  real  dkm                   !
  real  hkdepth               !length scale for decrease of hksati
  real  xksat                 !maximum hydraulic conductivity of soil [mm/s]

  real  dzlak(1:msl), zlak(1:msl)
  integer idlak

  integer iw, id 
!-----------------------------------------------------------------------
! soil layer thickness, depths
! Spatial discretization for soil, Unit is meter

      scalez = 0.025

! Define layer structure, vertical profiles of Ksat, and root fractions.
! here "0" refers to soil sourface, "msl" refers to the bottom of model soil
      do k = 1, kpt

         idlak = 1                                      !assumed all lakes are deep lake
         if(ivt(k) == 17)then                           !LAKE
            if(idlak == 1) then
               dzlak = (/1., 2., 3., 4., 5., 7., 7., 7., 7., 7./)
                zlak = (/0.5, 1.5, 4.5, 8.0, 12.5, 18.5, 25.5, 32.5, 39.5, 46.5/)
            else
               dzlak = (/.25, 0.5, 0.75, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)     ! shallow lake
                zlak = (/ 0.125,  0.5,  1.125,  2.,  3.,  4.,  5.,  6.,  7.,  8./)
            end if

            zsoih(0,k)   = 0.
            do j = 1, msl
               dzsoi(j,k) = dzlak(j)
                zsoi(j,k) =  zlak(j)
                zsoih(j,k) = zsoih(j-1,k)+dzsoi(j,k)
            enddo

         else                                           !NON LAKE LAND
            do j = 1, msl
              zsoi(j,k) = scalez*(exp(0.5*(j-0.5))-1.)     !node depths
            end do

            dzsoi(1,k)  = 0.5*(zsoi(1,k)+zsoi(2,k))        !=zsoih(1)
            dzsoi(msl,k)= zsoi(msl,k)-zsoi(msl-1,k)
            do j = 2,msl-1
              dzsoi(j,k)= 0.5*(zsoi(j+1,k)-zsoi(j-1,k))    !thickness b/n two interfaces
            end do

            zsoih(0,k)   = 0.
            zsoih(msl,k) = zsoi(msl,k) + 0.5*dzsoi(msl,k)
            do j = 1, msl-1
               zsoih(j,k)= 0.5*(zsoi(j,k)+zsoi(j+1,k))     !interface depths
            enddo
         endif
      enddo

!-----------------------------------------------------------------------
! soil thermal and hydraulic properties

      do k = 1, kpt
         if (ivt(k)/=15 .AND. ivt(k)/=17) then  ! not ice, lake
            do j = 1, msl

! (1) modification 1
! RJ forest
!             if(zsoih(j,k) <= 0.2)then
!                sand(k) = 88. + (82.-88.)/(.2-0.)*(zsoih(j,k)-0.)
!                clay(k) = 4.  + (10. - 4.)/(.2-0.)*(zsoih(j,k)-0.)
!             else if(zsoih(j,k) <= .4)then
!                sand(k) = 82. + (77.-82.)/(.4-.2)*(zsoih(j,k)-.2)
!                clay(k) = 10.  + (16. - 10.)/(.4-.2)*(zsoih(j,k)-.2)
!             else if(zsoih(j,k) <= .6)then
!                sand(k) = 77. + (63.-77.)/(.6-.4)*(zsoih(j,k)-.4)
!                clay(k) = 16.  + (18.- 16.)/(.6-.4)*(zsoih(j,k)-.4)
!             else if(zsoih(j,k) <= .8)then
!                sand(k) = 63. + (58.-63.)/(.8-.6)*(zsoih(j,k)-.6)
!                clay(k) = 18. + (35. -18.)/(.8-.6)*(zsoih(j,k)-.6)
!             else if(zsoih(j,k) <= 1.)then
!                sand(k) = 58. + (58.-58.)/(1.-.8)*(zsoih(j,k)-.8)
!                clay(k) = 35. + (35.-35.)/(1.-.8)*(zsoih(j,k)-.8)
!             else if(zsoih(j,k) <= 1.2)then
!                sand(k) = 58. + (53.-58.)/(1.2-1.)*(zsoih(j,k)-1.)
!                clay(k) = 35. + (36.-35.)/(1.2-1.)*(zsoih(j,k)-1.)
!             else 
!                sand(k) = 53.
!                clay(k) = 36.
!             endif
 
! RJ pasture
!            if(zsoih(j,k) <= 0.1)then
!               sand(k) = 85.
!               clay(k) = 7.
!            else if(zsoih(j,k) <= 1.0)then
!               sand(k) = 85. + (53.-85.)/(1.-0.1)*(zsoih(j,k)-0.1)
!               clay(k) = 7. + (35.-7.)/(1.-0.1)*(zsoih(j,k)-0.1)
!            else if(zsoih(j,k) <= 1.5)then
!               sand(k) = 53.
!               clay(k) = 35.
!            else
!               sand(k) = 33.
!               clay(k) = 17.
!            endif
!
! -----------------------------------------------------
             porsl(j,k) = 0.489 - 0.00126*sand(k)
              phi0(j,k) = 10. * ( 10.**(1.88-0.0131*sand(k)) )
               bsw(j,k) = 2.91 + 0.159*clay(k)

                 dkm    = (8.80*sand(k)+2.92*clay(k)) / (sand(k)+clay(k))       ! W/(m K)
              csol(j,k) = (2.128*sand(k)+2.385*clay(k))/(sand(k)+clay(k))*1.e6  ! J/(m3 K)

                     bd = (1.- porsl(j,k))*2.7e3
              dkmg(j,k) = dkm**(1.-porsl(j,k))
            dksatu(j,k) = dkmg(j,k)*0.57**porsl(j,k)
             dkdry(j,k) = (.135*bd + 64.7) / (2.7e3 - 0.947*bd)

                  xksat = 0.0070556 * ( 10.**(-0.884+0.0153*sand(k)) ) ! mm/s

! ABRACOS forest
!             if(zsoih(j,k) <= 0.2)then
!               porsl(j,k) = 0.483
!                xksat = 63./3600.
!             else if(zsoih(j,k) <= .4)then
!                porsl(j,k) = 0.483 + (0.305-0.483)/(.4-.2)*(zsoih(j,k)-.2)
!                xksat = 63./3600.  + (66./3600. - 63./3600.)/(.4-.2)*(zsoih(j,k)-.2)
!             else if(zsoih(j,k) <= .6)then
!                porsl(j,k) = 0.305 + (0.343-0.305)/(.6-.4)*(zsoih(j,k)-.4)
!                xksat = 66./3600.  + (10./3600.- 66./3600.)/(.6-.4)*(zsoih(j,k)-.4)
!             else if(zsoih(j,k) <= .8)then
!                porsl(j,k) = 0.343 +(0.397-0.343)/(.8-.6)*(zsoih(j,k)-.6)
!                xksat = 10./3600. + (5./3600. -10./3600.)/(.8-.6)*(zsoih(j,k)-.6)
!             else if(zsoih(j,k) <= 1.)then
!                porsl(j,k) = 0.397 + (0.410-0.397)/(1.-.8)*(zsoih(j,k)-.8)
!                xksat = 5./3600.   
!             else if(zsoih(j,k) <= 1.2)then
!                porsl(j,k) = 0.410 + (0.408-0.410)/(1.2-1.)*(zsoih(j,k)-1.)
!                xksat = 5./3600.
!             else
!                porsl(j,k) = 0.418
!                xksat = 5./3600.
!             endif
 
! Define the vertical profile of saturated hydraulic conductivity
! from 1 (model soil bottom) to msl (soil surface layer)
            hkdepth = 0.5        !m length scale for Ksat decrease
            hksati(j,k) = xksat  !*exp(-zsoih(j,k)/hkdepth)
            end do

! only for ABRACOS forest
            hksati(10,k) = 0.

! Saturated soil albedo for visible beam
            albsol(k) = solour(isc(k))

         else                                   ! ice/glacier, lake
            do j = 1, msl
             porsl(j,k) = 1.
              phi0(j,k) = -999.
               bsw(j,k) = -999.

              dkmg(j,k) = -999.
              csol(j,k) = -999.

            dksatu(j,k) = -999.
             dkdry(j,k) = -999.
            hksati(j,k) = -999.

              albsol(k) = -999.
            end do

         endif
      end do

!-----------------------------------------------------------------------
! vegetation static parameters

      do k = 1, kpt
                 i = ivt(k)
            z0m(k) = z0m_table   (i)  
         displa(k) = displa_table(i)
         sqrtdi(k) = sqrtdi_table(i)  

! ABRJ -----------------------------
         effcon(k) = 0.08
         vmax25(k) = 68.e-6    !100.e-6
         slti(k)   = 0.2
         hlti(k)   = 288.16
         shti(k)   = 0.3
         hhti(k)   = 313.16
         trda(k)   = 1.3
         trdm(k)   = 328.16
         trop(k)   = 298.16
         gradm(k)  = 9.0
         binter(k) = 0.01
         extkn(k)  = 0.5
 
         chil(k) = 0.1
         ref(1,1,k)=0.1; ref(1,2,k)=0.16; ref(2,1,k)=0.45;  ref(2,2,k)=0.39
         tran(1,1,k)=0.05;  tran(1,2,k)=0.001; tran(2,1,k)=0.25; tran(2,2,k)=0.001
 
! ABNS -----------------------------
!        effcon(k) = 0.05
!        vmax25(k) = 3.e-5
!        slti(k)   = 0.2
!        hlti(k)   = 288.16
!        shti(k)   = 0.3
!        hhti(k)   = 313.16
!        trda(k)   = 1.3
!        trdm(k)   = 328.16
!        trop(k)   = 298.16
!        gradm(k)  = 4.0
!        binter(k) = 0.04
!        extkn(k)  = 0.78
! 
!        chil(k) = -0.3
!        ref(1,1,k)=0.105; ref(1,2,k)=0.36; ref(2,1,k)=0.58;  ref(2,2,k)=0.58
!        tran(1,1,k)=0.07;  tran(1,2,k)=0.22; tran(2,1,k)=0.25; tran(2,2,k)=0.38
!
! Boreas forest -----------------------------
!        effcon(k) = 0.08
!        vmax25(k) = 1.e-4
!        slti(k)   = 0.2
!        hlti(k)   = 278.16
!        shti(k)   = 0.3
!        hhti(k)   = 303.16
!        trda(k)   = 1.3
!        trdm(k)   = 328.16
!        trop(k)   = 298.16
!        gradm(k)  = 9.0
!        binter(k) = 0.01
!        extkn(k)  = 0.78
!
!        chil(k) = -0.3
!        ref(1,1,k)=0.07; ref(1,2,k)=0.16; ref(2,1,k)=0.35;  ref(2,2,k)=0.39
!        tran(1,1,k)=0.05;  tran(1,2,k)=0.001; tran(2,1,k)=0.1; tran(2,2,k)=0.001

! Boreas grass -----------------------------
!        effcon(k) = 0.05
!        vmax25(k) = 3.e-5
!        slti(k)   = 0.2
!        hlti(k)   = 288.16
!        shti(k)   = 0.3
!        hhti(k)   = 313.16
!        trda(k)   = 1.3
!        trdm(k)   = 328.16
!        trop(k)   = 298.16
!        gradm(k)  = 4.0
!        binter(k) = 0.04
!        extkn(k)  = 0.78
! 
!        chil(k) = -0.3
!        ref(1,1,k)=0.105; ref(1,2,k)=0.36; ref(2,1,k)=0.58;  ref(2,2,k)=0.58
!        tran(1,1,k)=0.07;  tran(1,2,k)=0.22; tran(2,1,k)=0.25; tran(2,2,k)=0.38

! ----------------------------------

         albvgs(k) = albvgs_table(i)
         albvgl(k) = albvgl_table(i)


! root fraction (computing from surface, d is depth in meter):
! Y = 1 -1/2 (exp(-ad)+exp(-bd) under the constraint that
! Y(d =0.1m) = 1-beta^(10 cm) and Y(d=d_obs)=0.99 with beta & d_obs
! given in Zeng et al. (1998).

         do j = 1, msl-1
            rootfr(j,k) = .5*( exp(-roota(i)*zsoih(j-1,k))  &
                             + exp(-rootb(i)*zsoih(j-1,k))  &
                             - exp(-roota(i)*zsoih(j  ,k))  &
                             - exp(-rootb(i)*zsoih(j  ,k)) )
         end do
         rootfr(msl,k) =  .5*( exp(-roota(i)*zsoih(msl-1,k))&
                             + exp(-rootb(i)*zsoih(msl-1,k)))

         do j=1,10
         rootfr(j,k) = dzsoi(j,k)/zsoih(10,k)
         enddo

      enddo

  END SUBROUTINE clmtci