clmtvi.f90

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

  SUBROUTINE clmtvi (kpt   ,msl  ,maxsn  ,lun_ini,&
                     ivt   ,z0m  ,albsol ,albvgs ,albvgl ,chil  ,ref ,tran ,porsl ,&
                     glai  ,gsai ,dlon   ,dlat   ,rdlsf  ,rdlai ,&
                     mcsec ,jday ,xerr   ,zerr ,&
                     dz    ,z    ,zi     ,tss    ,wliq   ,wice ,&
                     rootr ,tlsun,tlsha  ,ldew   ,sag    ,scv  ,snowdp ,&
                     etrc  ,tg   ,albg   ,albv   ,alb    ,ssun ,ssha ,tranc,thermk,extkb,extkd, &
                     cosz  ,green,fveg   ,fsno   ,sigf   ,lai  ,sai ,&
                     snl)

!=======================================================================
!      Source file: clmtvi.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! 1, the observational initial values are the first choice
! 2, if there are no observational initial data, arbitrarily set initial fields,
!    and then spin up the first model year to equilibrium state, and take
!    the equilibrium state variables as the initial values.
!
!    the arbitrary initial data are created following the rule:
!   (1) foliage temperature is initially set to lowest atmospheric model air-temperature.
!   (2) canopy water storage set to zero.
!   (3) soil temperatures are initialized as in bucket type parameterizations
!       using the lowest atmospheric model air-temperature and a climatological
!       deep-ground temperature.
!   (4) soil moistures are initialized to a percentage of field capacity, the
!       percent of liquid water and ice lens are determined by the layer temperatures.
!   (5) if know the depth of snow, subdivide the snow pack up to five layers which
!       follows the rule: from top layer to bottom layer
!       minimum thickness:     0.010, 0.015, 0.025, 0.055, 0.115 (m),
!       and maximum thickness: 0.02, 0.05, 0.11, 0.23, and  >0.23 m,
!       the snow layer temperature set to surface air temperature,
!       if air temperature great than freezing point, set to 273.16 K.
!       if no any information on snow is available, snow mass for areas of
!       permanent land ice is initially set to 50000 kg m-2. Other areas,
!       all snow related variables are set to 0.
!=======================================================================

  IMPLICIT NONE 

! Dummy argument
  logical, INTENT(in) :: &
        rdlsf,           &! true if read in initial data set
        rdlai             ! true if read in LAI data set

  integer, INTENT(in) :: & 
        kpt,             &! total number of clm land points, including subgrid points
        msl,             &! soil layer number
        maxsn,           &! maximum snow layer number
        jday,            &! julian day
        mcsec             ! current seconds (0 - 86400)

  real, INTENT(in) :: &
        dlon(kpt),       &! longitude in radians: + = EH
        dlat(kpt)         ! latitude in radians: + = NH

  real, INTENT(in) :: &
        glai(18,18,12),  &!
        gsai(18,18,12)    !

  real, INTENT(inout) :: &
        dz(maxsn+1:msl, 1:kpt),  &! layer thickiness [m]
        z (maxsn+1:msl, 1:kpt),  &! node depth [m]
        zi(maxsn  :msl, 1:kpt)    ! interface depth [m]

  integer, INTENT(in) :: &
        ivt(kpt)          ! index for land cover type [-]
 
  real, INTENT(in) :: &
        albsol(kpt),     &! soil albedo for different coloured soils [-]
        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)
        z0m(kpt),        &! roughness
        porsl(1:msl,1:kpt)        ! porosity of soil

  integer, INTENT(out) :: &
        snl(kpt)          ! number of snow layers

  integer, INTENT(in) :: &
        lun_ini           ! logical unit number of initial data

  real, INTENT(out) :: &
        tss (maxsn+1:msl,1:kpt), &! soil temperature [K]
        wliq(maxsn+1:msl,1:kpt), &! liquid water in layers [kg/m2]
        wice(maxsn+1:msl,1:kpt), &! ice lens in layers [kg/m2]
        tlsun(kpt),      &! sunlit leaf temperature [K]
        tlsha(kpt),      &! shaded leaf temperature [K]
        ldew(kpt),       &! depth of water on foliage [mm]
        sag(kpt),        &! non dimensional snow age [-]
        scv(kpt),        &! snow cover, water equivalent [mm]
        snowdp(kpt),     &! snow depth [meter]
        rootr(1:msl,1:kpt),  &! root resistance of a layer, all layers add to 1.0
        etrc(kpt),       &! maximum possible transpiration rate [mm/s]
        tg(kpt)           ! ground surface temperature [K]

  real, INTENT(out) :: &
        xerr(1:kpt),     &! accumulation of the error of water banace
        zerr(1:kpt)       ! accumulation of the error of energy balance

  real, INTENT(out) :: &
        green(kpt),      &!
        fveg(kpt),       &! fraction of vegetated cover
        fsno(kpt),       &! fraction of snow cover on ground
        sigf(kpt),       &! fraction of veg cover, excluding snow-covered veg [-]
        lai(kpt),        &! leaf area index
        sai(kpt)          ! stem area index

  real, INTENT(out) :: &
        cosz(kpt),       &! cosine of solar zenith angle
        albg(2,2,kpt),   &! albedo, ground [-]
        albv(2,2,kpt),   &! albedo, vegetation [-]
        alb (2,2,kpt),   &! averaged albedo [-]
        ssun(2,2,kpt),   &!
        ssha(2,2,kpt),   &!
        tranc(2,2,kpt),  &!
        thermk(kpt),     &!
        extkb(kpt),      &!
        extkd(kpt)        !

! local
  real  ssw(kpt),        &! water volumetric content of soil surface layer
        wt(kpt),         &! fraction of vegetation covered by snow [-]
        calday            ! current julian day including fraction

  integer i, k

!-----------------------------------------------------------------------
! initialize water and temperature based on:
!    o rdlsf = true : read initial data set
!    o rdlsf = false: arbitrary initialization
!-----------------------------------------------------------------------
! (1) read initial data set

      if (rdlsf) then

         read(lun_ini) (snl      (k),k=1,kpt)
         read(lun_ini) ((dz    (i,k),i=-4,10),k=1,kpt)
         read(lun_ini) ((z     (i,k),i=-4,10),k=1,kpt)
         read(lun_ini) ((zi    (i,k),i=-5,10),k=1,kpt)
         read(lun_ini) ((tss   (i,k),i=-4,10),k=1,kpt)
         read(lun_ini) ((wliq  (i,k),i=-4,10),k=1,kpt)
         read(lun_ini) ((wice  (i,k),i=-4,10),k=1,kpt)

         read(lun_ini) (ldew     (k),k=1,kpt)
         read(lun_ini) (tlsun    (k),k=1,kpt)
         read(lun_ini) (tlsha    (k),k=1,kpt)
         read(lun_ini) (scv      (k),k=1,kpt)
         read(lun_ini) (snowdp   (k),k=1,kpt)
         read(lun_ini) (tg       (k),k=1,kpt)

         sag(1:kpt) = 0.

         close (lun_ini)

!-----------------------------------------------------------------------
! (2) arbitrary initialization, set water and temperatures to constant values
!-----------------------------------------------------------------------
      else
         ldew  (1:kpt) = 1.  
         scv   (1:kpt) = 0.
         sag   (1:kpt) = 0.
         snowdp(1:kpt) = 0.
         tlsun (1:kpt) =283.0 
         tlsha (1:kpt) =283.0 
         tg    (1:kpt) =283.0 

         snl(1:kpt) = 0 

         do k = 1, kpt
            do i = 1, msl
               tss(i,k) = 283.
               wice(i,k) = 0.
               if(ivt(k)==15 .OR. ivt(k)==17)then
                  wliq(i,k)=dz(i,k)*1000.
               else
                  wliq(i,k) = dz(i,k)*1.*porsl(i,k)*1000.
               endif
            enddo
         enddo

         tss(-4:0,1:kpt) = -999.
         wice(-4:0,1:kpt) = -999.
         wliq(-4:0,1:kpt) = -999.
         dz(-4:0,1:kpt) = -999.
         z(-4:0,1:kpt) = -999.
         zi(-5:-1,1:kpt) = -999.

      endif

!-----------------------------------------------------------------------
! (3) the remaining variables are not part of the initial data set
!-----------------------------------------------------------------------

! vegetation cover fraction and LAI, SAI 
      calday = float(jday) + float(mcsec)/86400.
      call laiini(kpt,ivt,dlat,calday,glai,gsai,green,fveg,lai,sai)

! fraction of snow cover
      call fractsnow (kpt,fveg,z0m,snowdp,wt,sigf,fsno)

! cosine of solar zenith angle
      call clmzen (kpt, calday, dlon, dlat, cosz)

! surface albedo
      do k = 1, kpt
         ssw(k) = 1.
         if(ivt(k) /= 17) ssw(k) = min(1.,1.e-3*wliq(1,k)/dz(1,k))
      enddo

      call albedo (kpt,ivt,albsol,albvgs,albvgl,&
                   chil,ref,tran,fveg,green,lai,sai,cosz,wt,fsno,ssw,tg,&
                   scv,sag,alb,albg,albv,ssun,ssha,tranc,thermk,extkb,extkd)

! initial set for rooting fraction and miximum of transpiration

      do k = 1, kpt
         if(ivt(k) == 17) cycle
         do i = 1, msl
           rootr(i,k) = 0.
         end do

         etrc(k) = 1.e-12
      end do

! initial set for the balance errors
      do k = 1, kpt
         xerr(k) = 0.
         zerr(k) = 0.
      enddo
    
  END SUBROUTINE clmtvi