laiini.f90

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

      subroutine laiini(kpt,ivt,dlat,calday,glai,gsai,green,fveg,lai,sai)

!=======================================================================
      implicit none

      integer, INTENT(in) :: &
            kpt               ! total number of clm land points, including subgrid points

      integer, INTENT(in) :: &
            ivt(kpt)          ! index for land cover type [-]

      real, INTENT(in) :: &
            dlat(kpt),       &! latitude in radians: + = NH
            calday,          &! current julian day including fraction
            glai(18,18,12),  &! LAI
            gsai(18,18,12)    ! SAI

      real, INTENT(out) :: &
            green(kpt),      &!
            fveg(kpt),       &! fraction of vegetation cover
            lai(kpt),        &! LAI
            sai(kpt)          ! SAI
      
      integer &
            i,               &!
            k,               &!
            ireg              ! regional index

      real  pie               ! pie
      real  frad              ! converts degrees to radians.

      real  ldoy              ! intermediate variable
      real  ndoy              ! intermediate variable
      real  cday              ! intermediate variable
      real  fact              ! intermediate variable

      real  rmon(12)          !
      data rmon/31.,59.,90.,120.,151.,181., &
                212.,243.,273.,304.,334.,365./

      logical observation
      real glai_obs(12), gsai_obs(12),fveg_obs(12)
!-----------------------------------------------------------------------
! (1) Valdai
!     data glai_obs/.00,0.00,0.60,0.70,2.11,4.18,4.07,2.54,0.80,0.60,0.50,0.00/
!     data gsai_obs /12*0.8/
!     data fveg_obs /12*1.0/
! (2) Cabauw
!     data glai_obs/0.456, 0.62, 0.726, 0.91, 1.296, 1.638, 1.458, 0.624, 0.66, 0.594, 0.55, 0.45/
!     data gsai_obs /12*4./
!     data fveg_obs /0.92, 0.93, 0.94, 0.95, 0.97, 0.98, 0.99, 0.98, 0.97, 0.96, 0.95, 0.93/
! (3) Hapex
!     data glai_obs/4*0.0, 1.0, 4*3.0, 3*0.0/
!     data gsai_obs /4*0.0, 5*0.5, 3*0.0/
!     data fveg_obs /4*0.0, 0.5, 4*0.9, 3*0.0/
! (4) Abracos forest
!     data glai_obs/4.48, 4.09, 3.74, 6*3.68, 3.85, 4.36, 4.60/
!     data gsai_obs /12*1./
!     data fveg_obs /12*1.0/
! (5) Abracos pasture
!     data glai_obs/1.42, 2.42, 3.00, 3.90, 3.22, 2.55, 1.86, 1.75, 1.66, 0.88, 0.72, 0.57/
!     data gsai_obs /12*0.5/
!     data fveg_obs /12*1.0/
! (6) Fife
!     data glai_obs/4*0.008, 1.105, 2.6, 2.7, 2.36, 1.3, 0.7, 0.1, 0.01/
!     data gsai_obs /12*0.8/
!     data fveg_obs /12*1.0/
! (7) Tucson
!     data glai_obs/1.55, 1.31, 1.20, 0.89, 0.808, 0.777, 0.806, 0.947, 1.085, 1.437, 1.679, 1.739/
!     data gsai_obs /12*2.0/
!     data fveg_obs /12*1.0/
! (8) Fort Peck
     data glai_obs/12*0.0/
     data gsai_obs /12*0.0/
     data fveg_obs /12*0.0/

!-----------------------------------------------------------------------
      observation = .true. 
!=======================================================================

      pie = 4.*atan(1.)
      frad = pie/180.     ! converts degrees to radians.

if(.not.observation)then
      do k = 1, kpt
         do i = 1, 18
            if((dlat(k)> -90.*frad + 10.*float(i-1)*frad) .AND. &
               (dlat(k)<=-90.*frad + 10.*float(i)*frad))then
               ireg = i
               exit
            endif
         enddo

         green(k) = 1.

         if((calday < 31./2.) .OR. (calday >= 365 - 31./2.) ) then
             ndoy = 365. + 31./2.
             ldoy = 365. - 31./2.
             if(calday < 31./2.) then
                cday = 365. + calday
             else
                cday = calday
             endif
             fact = (cday - ldoy)/(ndoy - ldoy)
             lai(k) = glai(ivt(k),ireg,12) &
                     + fact*(glai(ivt(k),ireg,1) - glai(ivt(k),ireg,12))
             if(ivt(k)/=15 .and. ivt(k)/=17 .and. ivt(k)/=18)then
                lai(k) = max(lai(k), 1.e-6) 
             endif
             sai(k) = gsai(ivt(k),ireg,12) &
                    + fact*(gsai(ivt(k),ireg,1) - gsai(ivt(k),ireg,12))
         else

             do i = 1, 11
                if(i == 1)then
                   ldoy = 31./2
                else
                   ldoy = rmon(i) - (rmon(i)-rmon(i-1))/2.
                endif

                if(i == 11)then
                   ndoy = 365. - 31./2.
                else
                   ndoy = rmon(i) + (rmon(i+1)-rmon(i))/2.
                endif

                if((calday >= ldoy) .AND. (calday < ndoy))then
                    fact = (calday - ldoy)/(ndoy - ldoy)
                    lai(k) = glai(ivt(k),ireg,i) &
                            + fact*(glai(ivt(k),ireg,i+1) - glai(ivt(k),ireg,i))
                    if(ivt(k)/=15 .and. ivt(k)/=17 .and. ivt(k)/=18)then
                       lai(k) = max(lai(k), 1.e-6) 
                    endif
                    sai(k) = gsai(ivt(k),ireg,i) &
                            + fact*(gsai(ivt(k),ireg,i+1) - gsai(ivt(k),ireg,i))
                    exit
                endif
             enddo
         endif

         fveg(k)=1.           ! mosaic scheme use only
         green(k) = 1.
          print*,'fveg',fveg
         if(ivt(k)==15 .or. ivt(k)>=17)then
            fveg(k)=0.; green(k)=0.
            print*,'fveg',fveg
         endif

      enddo
else
      do k = 1,kpt 
         green(k) = 1.

      if(k == 1)then

         if((calday < 31./2.) .OR. (calday >= 365 - 31./2.) ) then
             ndoy = 365. + 31./2.
             ldoy = 365. - 31./2.
             if(calday < 31./2.) then
                cday = 365. + calday
             else
                cday = calday
             endif
             fact = (cday - ldoy)/(ndoy - ldoy)
             lai(k) = glai_obs(12) + fact*(glai_obs(1) - glai_obs(12))
             if(ivt(k)/=15 .and. ivt(k)/=17 .and. ivt(k)/=18)then
                lai(k) = max(lai(k), 1.e-6) 
             endif
             sai(k) = gsai_obs(12) + fact*(gsai_obs(1) - gsai_obs(12))
             fveg(k) = fveg_obs(12) + fact*(fveg_obs(1) - fveg_obs(12))
         else
             do i = 1, 11
                if(i == 1)then
                   ldoy = 31./2
                else
                   ldoy = rmon(i) - (rmon(i)-rmon(i-1))/2.
                endif

                if(i == 11)then
                   ndoy = 365. - 31./2.
                else
                   ndoy = rmon(i) + (rmon(i+1)-rmon(i))/2.
                endif

                if((calday >= ldoy) .AND. (calday < ndoy))then
                    fact = (calday - ldoy)/(ndoy - ldoy)
                    lai(k) = glai_obs(i) + fact*(glai_obs(i+1) - glai_obs(i))
                    if(ivt(k)/=15 .and. ivt(k)/=17 .and. ivt(k)/=18)then
                       lai(k) = max(lai(k), 1.e-6)
                    endif
                    sai(k) = gsai_obs(i) + fact*(gsai_obs(i+1) - gsai_obs(i))
                    fveg(k) = fveg_obs(i) + fact*(fveg_obs(i+1) - fveg_obs(i))
                    exit
                endif
             enddo
         endif

      else
         lai(k)=0.;sai(k)=0.;fveg(k)=0.;green(k)=0.
      endif
      enddo
endif

      end subroutine laiini