albedo.f90

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

  SUBROUTINE 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) 

!=======================================================================
!      Source file: albedo.f90
!     Initilizated: Yongjiu Dai, September 15, 1999
!
! Calculates fragmented albedos (direct and diffuse) in
! wavelength regions split at 0.7um.
! 
! (1) soil albedos: as in BATS formulations, which are the function of
!     soil color and moisture in the surface soil layer
! (2) snow albedos: as in BATS formulations, which are inferred from
!     the calculations of Wiscombe and Warren (1980) and the snow model
!     and data of Anderson(1976), and the function of snow age, grain size,
!     solar zenith angle, pollution, the amount of the fresh snow
! (3) canopy albedo: has been developed to capture the essential features of
!     a two-stream approximation model while forgoing the complexity of the full
!     treatment. It combines soil and canopy albedo by simple rules that are
!     formulated to reduce to correct asymptotic limits for thick and thin canopies
!     and provide reasonable results for intermediate values of leaf area index
! (4) glacier albedos: as in BATS, which are set to constants (0.8 for visible beam,
!     0.55 for near-infrared)
! (5) lake and wetland albedos: as in BATS, which depend on cosine solar zenith angle,
!     based on data in Henderson-Sellers (1986). The frozen lake and wetland albedos
!     are set to constants (0.6 for visible beam, 0.4 for near-infrared)
! (6) over the snow covered tile, the surface albedo is estimated by a linear
!     combination of albedos for snow, canopy and bare soil (or lake, wetland, glacier).
!=======================================================================

  Use phycon_module

  IMPLICIT NONE

!------------------------- Dummy Arguments -----------------------------
! ground cover index
  integer, INTENT(in) :: &
        kpt,       &! number of clm land points, including subgrid points
        ivt(kpt)    ! index for land cover type [-]

! parameters
  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)
        fveg(kpt),      &! fractional vegetation cover [-]
        green(kpt),     &! green leaf fraction
        lai(kpt),       &! leaf area index (LAI+SAI) [m2/m2]
        sai(kpt)         ! stem area index (LAI+SAI) [m2/m2]


! variables
  real, dimension(kpt), INTENT(in) :: &
        cosz,      &! cosine of solar zenith angle [-]
        wt,        &! fraction of vegetation covered by snow [-]
        fsno,      &! fraction of soil covered by snow [-]
        ssw,       &! water volumetric content of soil surface layer [m3/m3]
        scv,       &! snow cover, water equivalent [mm]
        sag,       &! non dimensional snow age [-]
        tg          ! ground surface temperature [K]

! output
  real, INTENT(out) :: &
        alb(2,2,kpt),   &! averaged albedo [-]
        albg(2,2,kpt),   &! albedo, ground
        albv(2,2,kpt),   &! albedo, vegetation [-]
        ssun(2,2,kpt),   &! sunlit canopy absorption for solar radiation
        ssha(2,2,kpt),   &! shaded canopy absorption for solar radiation,
                          ! normalized by the incident flux
        tranc(2,2,kpt),  &! canopy transmittances for solar radiation
        thermk(kpt),     &! canopy gap fraction for tir radiation
        extkb(kpt),      &! (k, g(mu)/mu) direct solar extinction coefficient
        extkd(kpt)        ! diffuse and scattered diffuse PAR extinction coefficient

!-------------------------- Local variables ----------------------------

  integer          &!
        iw,        &! wavelength (1=visible, 2=near-infrared)
        id,        &! 1=direct, 2=diffuse
        k           ! looping indx

   real age,       &! factor to reduce visible snow alb due to snow age [-]
        albg0,     &! temporary varaiable [-]
        albsno(2,2),&! snow albedo [-]
        albv0(2),  &! vegetation albedo [-]
        alwet,     &! decrease in soil albedo due to wetness [-]
        beta0,     &! upscattering parameter for direct beam [-]
        cff,       &! snow alb correction factor for zenith angle > 60 [-]
        conn,      &! constant (=0.5) for visible snow alb calculation [-]
        cons,      &! constant (=0.2) for nir snow albedo calculation [-]
        czen,      &! cosine of solar zenith angle > 0 [-]
        czf,       &! solar zenith correction for new snow albedo [-]
        dfalbl,    &! snow albedo for diffuse nir radiation [-]
        dfalbs,    &! snow albedo for diffuse visible solar radiation [-]
        dralbl,    &! snow albedo for visible radiation [-]
        dralbs,    &! snow albedo for near infrared radiation [-]
        fsol1,     &! solar flux fraction for wavelength < 0.7 micron [-]
        fsol2,     &! solar flux fraction for wavelength > 0.7 micron [-]
        lsai,      &! leaf and stem area index (LAI+SAI) [m2/m2]
        scat(2),   &! single scattering albedo for vir/nir beam [-]
        sl,        &! factor that helps control alb zenith dependence [-]
        snal0,     &! alb for visible,incident on new snow (zen ang<60) [-]
        snal1,     &! alb for NIR, incident on new snow (zen angle<60) [-]
        tdiffs,    &! difference of air temperature and freezing temp [K]
        tff,       &! exp(-LSAI)
        tffd,      &! exp(-0.5*LSAI/czen)
        ti,        &! correction due to scattering
        upscat,    &! upward scattered fraction for direct beam [-]
        zkat(2),   &! temporary
        zkatd(2)    ! temporary
     
! ----------------------------------------------------------------------
! 1. Initial set
! ----------------------------------------------------------------------
! Division of solar flux for wavelength less or greater than 0.7 micron
      fsol1 = 0.5      ! shortwave
      fsol2 = 0.5      ! longwave

! Short and long wave albedo for new snow
      snal0 = 0.85     ! shortwave
      snal1 = 0.65     ! long wave

! Set initial leaf scattering reflectance. Note: "scat" may use different
! value for different vegetation latter
      beta0 = 0.5
      scat(1) = 0.15
      scat(2) = 0.85

! ----------------------------------------------------------------------
      do k = 1, kpt

! Set default soil and vegetation albedos and solar absorption
          do id=1,2
             do iw=1,2
                alb (iw,id,k) = 0. ! averaged
                albg(iw,id,k) = 0. ! ground
                albv(iw,id,k) = 0. ! vegetation
                ssun(iw,id,k) = 0.
                ssha(iw,id,k) = 0.
                tranc(iw,id,k)= 0.
             enddo
          enddo

          thermk(k) = 1.e-3
          extkb(k) = 1.e-6
          extkd(k) = 0.718

          lsai = lai(k) + sai(k)
          if(cosz(k)<=0.) CYCLE

! Get albedo over land if solar angle is poitive
! Cosine of zenith angle

          czen = max(cosz(k),0.001)

! Set initial snow albedo
          albsno(1,1) = 0.
          albsno(2,1) = 0.
          albsno(1,2) = 0.
          albsno(2,2) = 0.

! ----------------------------------------------------------------------
! 2. Albedo for snow cover.
!    Snow albedo depends on snow-age, zenith angle, and thickness
!    of snow age gives reduction of visible radiation
! ----------------------------------------------------------------------
          if(scv(k) > 0.) then

             cons = 0.2
             conn = 0.5

! Sl helps control albedo zenith dependence
             sl  = 2.0

! Correction for snow age
             age = 1.-1./(1.+sag(k))
             dfalbs = snal0*(1.-cons*age)

! Czf corrects albedo of new snow for solar zenith

             cff    = ((1.+1./sl)/(1.+czen*2.*sl )- 1./sl)
             cff    = max(cff,0.)
             czf    = 0.4*cff*(1.-dfalbs)
             dralbs = dfalbs+czf
             dfalbl = snal1*(1.-conn*age)
             czf    = 0.4*cff*(1.-dfalbl)
             dralbl = dfalbl+czf
       
             albsno(1,1) = dralbs
             albsno(2,1) = dralbl
             albsno(1,2) = dfalbs
             albsno(2,2) = dfalbl

          endif

! ----------------------------------------------------------------------
! 3. Get albedo over land
! ----------------------------------------------------------------------
! 3.1 Bare soil albedos, depends on moisture
          if(ivt(k)/=11 .AND. ivt(k)/=15 .AND. ivt(k)/=17) then    
                                ! not wetland, permanent ice and water
             alwet = max((11.-40.0*ssw(k)),0.) *0.01
             alwet = min(alwet,albsol(k))
             albg0 = albsol(k)+alwet
             albg(1,1,k) = albg0
             albg(2,1,k) = 2.*albg0
     
             do iw=1,2                        ! Diffused albedos for bare soil
                albg(iw,2,k) = albg(iw,1,k)
             enddo

! 3.2 Albedos for permanent ice sheet. 
          else if (ivt(k)==15) then           ! permanent ice sheet

             albg(1,1,k) = 0.8
             albg(1,2,k) = 0.8
             albg(2,1,k) = 0.55
             albg(2,2,k) = 0.55

! 3.3 Albedo for wetland (swamps, rice paddies etc) and inland water
          else if (ivt(k)==11 .OR. ivt(k)==17) then             
             albg0 = 0.05/(czen+0.15)
             albg(1,1,k) = albg0
             albg(1,2,k) = albg0
             albg(2,1,k) = albg0
             albg(2,2,k) = albg0

             if(tg(k) < tfrz) then          ! frozen lake and wetland
                albg(1,1,k) = 0.6
                albg(1,2,k) = 0.6
                albg(2,1,k) = 0.4
                albg(2,2,k) = 0.4
             endif

          end if

! 3.4 Correction due to snow cover

          do id=1,2
             do iw=1,2
                albg(iw,id,k) = (1.-fsno(k))*albg(iw,id,k) &
                              + albsno(iw,id)*fsno(k)
             enddo
          enddo

! ----------------------------------------------------------------------
! 4. Canopy albedos  [R. E. Dickinson, 1998, personal communication]
! ----------------------------------------------------------------------
          if(fveg(k) > 0.001) then
 
!*Direct albedos for vegetation
!*           albv0(1) = albvgs(k)
!*           albv0(2) = albvgl(k)
!*
!*Upward scattered fraction for direct beam
!*           do iw = 1, 2
!*           upscat = scat(iw)*beta0
!*           zkat(iw) = upscat*lsai/albv0(iw)
!*           zkatd(iw) = zkat(iw)/czen
!*           enddo
!*
!*           tff = exp(-lsai)
!*           tffd = exp(-0.5*lsai/czen)
!*
!*Diffused albedos for vegetation 
!*           do iw = 1,2
!*              ti = 1.-exp(-2.*zkat(iw))
!*              albv(iw,2,k) = albv0(iw)*ti + albg(iw,2,k)*tff*tff
!*              absv(iw,2,k) = 1.- albv(iw,2,k) - (1.-albg(iw,2,k))*tff
!*           enddo
!*
!*Direct albedos for vegetation
!*           do iw = 1,2
!*              ti = 1. - exp(-zkatd(iw))
!*              albv(iw,1,k) = albv0(iw)*ti + albg(iw,1,k)*tffd*tff
!*              absv(iw,1,k) = 1.-albv(iw,1,k) - (1.-albg(iw,1,k))*tffd
!*           enddo
!*
!*Effective albedo over vegetation with snow
!*           do id = 1, 2
!*              do iw = 1, 2
!*                 albv(iw,id,k) = (1.-wt(k))*albv(iw,id,k) &
!*                               +     wt(k)*albsno(iw,id)
!*              enddo
!*           enddo
 
            ! two stream approximation
            !-----------------------------------------------------------
             call twostr ( chil(k), ref(1,1,k), tran(1,1,k), green(k), &
                           lai(k), sai(k), czen, albg(1,1,k), albv(1,1,k), &
                  tranc(1,1,k), thermk(k), extkb(k), extkd(k), ssun(1,1,k), ssha(1,1,k) ) 

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

          end if

! ----------------------------------------------------------------------
! 5. Weighted albedo over a pixel for shortwave and longwave
! ----------------------------------------------------------------------

          do id = 1, 2
             do iw = 1, 2
                if(fveg(k) <= 0.001)then
                   alb(iw,id,k) = albg(iw,id,k)
                else
                   alb(iw,id,k) = (1.-fveg(k))*albg(iw,id,k) &
                                +  fveg(k)*albv(iw,id,k)
                endif
             enddo
          enddo

! ----------------------------------------------------------------------
      enddo 

  END SUBROUTINE albedo