average.f

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

! ----------------------------------------------------------------------
! for CLM daily/monthly/annual average mean
! ----------------------------------------------------------------------
      IMPLICIT NONE
      integer msl,i,j,m,mstep, mday, jday, k, mi, intf,mn
      parameter (msl=10)
!     parameter (mstep=315552, mday=6570, intf=48, mi=1)    ! Valdai data
!     parameter (mstep=17520, mday=365, intf=48, mi=1)      ! Cabauw / Hapex
!     parameter (mstep=17520, mday=365, intf=24, mi=9)      ! Arme
      parameter (mstep=17544, mday=365, intf=24, mi=1)      ! Abracos
!     parameter (mstep=35088, mday=365, intf=48, mi=5)      ! fife data
!     parameter (mstep=26280, mday=365, intf=72, mi=6)      ! Tucson data
!     parameter (mstep=17520, mday=365, intf=48, mi=1)      ! Boreas

      real calday     ! current julian day including fraction

! Radiative fluxes
      real solar,     ! incident solar radiation (W/m2)
     $     sabv,      ! solar absorbed by vegetation (W/m2)
     $     sabg,      ! solar absorbed by ground (W/m2)
     $     frl,       ! downward infrared (longwave) radiation (W/m2)
     $     olrg,      ! outgoing longwave radiation from canopy + ground

! Sensible fluxes
     $     fsena,     ! sensible heat to atmosphere (total) (W/m2)
     $     fsenl,     ! sensible heat from leaves to canopy air (W/m2)
     $     fseng,     ! sensible heat flux from ground (W/m2)

! Latent heat fluxes
     $     fevpa,     ! evapotranspiration to atmosphere (total) (W/m2)
     $     fevpl,     ! evaporation+transpiration from leaves (mm/s)
     $     etr,       ! transpiration rate (mm/s)
     $     fevpg,     ! evaporation heat flux from ground (mm/s)

     $     fgrnd,     ! ground heat flux ( positive surface towards soil) (W/m2)
     $     taux,      ! stress
     $     tauy,      ! stress

! Water balance
     $     prcp,      ! precipitation (mm h2o/s)
     $     rsur,      ! surface runoff (mm h2o/s)
     $     rnof,      ! total runoff (mm h2o/s)
     $     ldew,      ! depth of water on foliage (mm)
     $     scv,       ! snow cover, water equivalent (mm)
     $     snowdp,    ! snow depth (meter)

! Photosynthesis
     $     rst,       !
     $     assim,     !
     $     respc,     !

! State variables
     $     tm,        ! air temperature (K)
     $     trad,      ! radiative temperature (K)
     $     tref,      ! reference height temperature (K)
     $     tlsun,     ! leaf temperature (K)
     $     tlsha,     ! leaf temperature (K)
     $     tg,        ! ground surface temperature (K)

     $     wind,      ! wind component (m/s)
     $     qm,        ! specific humidity (kg/kg)

     $     sigf,      ! fraction of veg cover, excluding snow-covered veg (-)
     $     fsno,      ! fraction of soil covered by snow (-)

     $     tss(msl),  ! soil temperature (K)
     $     wliq(msl), ! liquid water in layers (kg/m2)
     $     wice(msl)  ! ice lens in layers (water equivalent) (kg/m2)

      real X(65), Y(65)
      integer julia(12)
      data julia / 31,  59,  90, 120, 151, 181, 212,
     $            243, 273, 304, 334, 365 /

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

      m = 0
      mn = mi
      do i = 2, 65
         X(I) = 0.
         Y(I) = 0.
      enddo

!     open(11,file='VALDAI-output',form='unformatted')
!     open(12,file='CLM_valdai.daily',form='formatted')
!     open(13,file='CLM_valdai.monthly',form='formatted')

!     open(11,file='CAB-output',form='unformatted')
!     open(12,file='CLM_cabauw.daily',form='formatted')
!     open(13,file='CLM_cabauw.monthly',form='formatted')

!     open(11,file='HAP-output',form='unformatted')
!     open(12,file='CLM_hapex.daily',form='formatted')
!     open(13,file='CLM_hapex.monthly',form='formatted')

!     open(11,file='ARM-output',form='unformatted')
!     open(12,file='CLM_arme.daily',form='formatted')
!     open(13,file='CLM_arme.monthly',form='formatted')

      open(11,file='ABRJ-output_twoleaf',form='unformatted')
      open(12,file='CLM_abrj_twoleaf.daily',form='formatted')
      open(13,file='CLM_abrj_twoleaf.monthly',form='formatted')

!     open(11,file='ABNS-output_oneleaf',form='unformatted')
!     open(12,file='CLM_abns.daily_oneleaf',form='formatted')
!     open(13,file='CLM_abns.monthly_oneleaf',form='formatted')

!     open(11,file='TUC-output',form='unformatted')
!     open(12,file='CLM_tuc.daily',form='formatted')
!     open(13,file='CLM_tuc.monthly',form='formatted')

!     open(11,file='FIFE-output',form='unformatted')
!     open(12,file='CLM_fife.daily',form='formatted')
!     open(13,file='CLM_fife.monthly',form='formatted')

!     open(11,file='NSAGRASS-output_oneleaf',form='unformatted')
!     open(12,file='CLM_nsagrass.daily_oneleaf',form='formatted')
!     open(13,file='CLM_nsagrass.monthly_oneleaf',form='formatted')
! ----------------------------------------------------------------------
      do 100 j = 1, mstep
      m = m+ 1

      read(11)                                 !
     $   calday,                               ! [1]
! ----------------------------------------------------------------------
! Fluxes
     $   solar,                                ! [2]  W/m2
     $   sabv,                                 ! [3]  W/m2
     $   sabg,                                 ! [4]  W/m2
     $   frl,                                  ! [5]  W/m2
     $   olrg,                                 ! [6]  W/m2

     $   fsena,                                ! [7]
     $   fsenl,                                ! [8]  W/m2
     $   fseng,                                ! [9]  W/m2

     $   fevpa,                                ! [10] W/m2  
     $   fevpl,                                ! [11] [mm/s]
     $   etr,                                  ! [12] [mm/s]
     $   fevpg,                                ! [13] [mm/s]

     $   fgrnd,                                ! [14] W/m2
     $   taux,                                 ! [15] (kg/m/s2)
     $   tauy,                                 ! [16] (kg/m/s2)

! ----------------------------------------------------------------------
! Water
     $   prcp,                                 ! [17] mm / s
     $   rsur,                                 ! [18] mm / s
     $   rnof,                                 ! [19] mm / s
     $   ldew,                                 ! [20] millimeter
     $   scv,                                  ! [21] millimeter
     $   snowdp,                               ! [22] meter

! ----------------------------------------------------------------------
! assimilation of CO2
     $   rst,                                  ! [23]
     $   assim,                                ! [24]
     $   respc,                                ! [25]

! ----------------------------------------------------------------------
! temperatures

     $   tm,                                   ! [26] K
     $   trad,                                 ! [27] K
     $   tref,                                 ! [28] K
     $   tlsun,                                ! [29] K
     $   tlsha,                                ! [30] K
     $   tg,                                   ! [31] K

     $   wind,                                 ! [32] m/s
     $   qm,                                   ! [33] kg/kg

     $   sigf,                                 ! [34] -
     $   fsno,                                 ! [35] -

     $  (tss (i),i=1, 10),                     ! [36 - 45] K
     $  (wliq(i),i=1, 10),                     ! [46 - 55] millimeter
     $  (wice(i),i=1, 10)                      ! [56 - 65] millimeter

         X(1)  = calday
         X(2)  = X(2)  + solar       
         X(3)  = X(3)  + sabv        
         X(4)  = X(4)  + sabg        
         X(5)  = X(5)  + frl         
         X(6)  = X(6)  + olrg        

         X(7)  = X(7)  + fsena       
         X(8)  = X(8)  + fsenl      
         X(9)  = X(9)  + fseng      

         X(10) = X(10) + fevpa      
         X(11) = X(11) + fevpl      
         X(12) = X(12) + etr        
         X(13) = X(13) + fevpg      

         X(14) = X(14) + fgrnd      
         X(15) = X(15) + taux
         X(16) = X(16) + tauy

         X(17) = X(17) + prcp       
         X(18) = X(18) + rsur       
         X(19) = X(19) + rnof       
         X(20) = X(20) + ldew      
         X(21) = X(21) + scv        
         X(22) = X(22) + snowdp    

         X(23) = X(23) + rst
         X(24) = X(24) + assim
         X(25) = X(25) + respc

         X(26) = X(26) + tm        
         X(27) = X(27) + trad       
         X(28) = X(28) + tref
         X(29) = X(29) + tlsun
         X(30) = X(30) + tlsha
         X(31) = X(31) + tg

         X(32) = X(32) + wind      
         X(33) = X(33) + qm        

         X(34) = X(34) + sigf      
         X(35) = X(35) + fsno      

         do i = 1,10
            X(35+i) = X(35+i) + tss(i) 
            X(45+i) = X(45+i) + wliq(i)
            X(55+i) = X(55+i) + wice(i)
         enddo

! daily average 
         if(m .eq. intf)then
            do i = 2, 65
               X(I) = X(I)/float(intf)
            enddo

            write(12,*) int(X(1)), (X(I),i=2,65)

C-------------------------
! monthly average
            do i = 2, 65
               Y(i) = Y(i) + X(i)
            enddo 
            if(int(X(1)) .eq. julia(mn))then
               do i = 2, 65
                  if(mn .le. 1)then
                     Y(i) = Y(i)/real(julia(mn))
                  else
                     Y(i) = Y(i)/real(julia(mn) - julia(mn-1))
                  endif
               enddo

               write(13,*) int(X(1)), (Y(i),i=2,65)
               mn= mn + 1
               if(int(X(1)).eq.julia(12)) mn = 1
               do i = 2, 65
                  Y(i) = 0.
               enddo
            endif

            m = 0
            do i = 2, 65
               X(I) = 0.
            enddo
C-------------------------
         endif

100   continue 
      close (11)
      close (12)
      close (13)

      end