clm_main.f90

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

! ======================================================================

      do k = 1, kpt
         scvold(k) = scv(k)          ! snow mass at previous time step
         newnode(k) = 0              ! signification when snow node will be initialized
         if(ivt(k) == 17) cycle
         if(snl(k) < 0)then
           do i = snl(k)+1, 0  ! ice fraction of snow at previous time step
              fiold(i,k) = wice(i,k)/(wliq(i,k)+wice(i,k)) 
           enddo
         endif
      enddo
!      print*,'initial set'
! ======================================================================
! [2] Canopy interception and precipitation onto ground surface
! ======================================================================
!
! 2.1 Add precipitation to leaf water 

      do k = 1, kpt
         if(ivt(k)/=15 .AND. ivt(k)/=17)then 
            call leafdew (dtime,rainrate(k),snowrate(k),tm(k),&
                           sigf(k),lai(k),sai(k),ldew(k),etrrun(k),tti(k))
         else
            ldew(k)=0.; etrrun(k)=0.; tti(k)=0.    ! land ice and lake
         endif
      enddo
!      print*,'leafdew'
! 2.2 Precipitation onto ground

      do k = 1, kpt
         if(ivt(k) == 17) cycle
         pg(k) = (1.-sigf(k))*(rainrate(k)+snowrate(k)) &
               + tti(k) + etrrun(k)                     ! kg/(m2 s)
         if(sigf(k) <= 0.001) pg(k) = rainrate(k) + snowrate(k)

! 2.3 The percentage of liquid water by mass, which is arbitrarily set to 
!     vary linearly with air temp, from 0% at 273.16 to 40% max at 275.16.

         if(iptype(k) <= 1)then

           flfall(k) = 1.        ! fraction of liquid water within falling precip.
           pg_snow(k) = 0.       ! ice onto ground (mm/s)
           pg_rain(k) = pg(k)    ! liquid water onto ground (mm/s)
           dz_snowf(k) = 0.      ! rate of snowfall, snow depth/s (m/s)

         else

           if(tm(k) <= tfrz)then
             flfall(k) = 0.
           else if(tm(k) <= tfrz+2.)then
             flfall(k) = -54.632 + 0.2*tm(k)
           else
             flfall(k) = 0.4
           endif
 
!        use Alta relationship, Anderson(1976); LaChapelle(1961), 
!        U.S.Department of Agriculture Forest Service, Project F, 
!        Progress Rep. 1, Alta Avalanche Study Center:Snow Layer Densification.

           if(tm(k) > tfrz + 2.)then
             bifall(k) =169.
           else if(tm(k) > tfrz - 15.)then
             bifall(k)=50. + 1.7*(tm(k) - tfrz + 15.)**1.5
           else
             bifall(k)=50.
           endif
 
           pg_snow(k) = pg(k)*(1.-flfall(k))                 
           pg_rain(k) = pg(k)*flfall(k)
           dz_snowf(k) = pg_snow(k)/bifall(k)                
           snowdp(k) = snowdp(k) + dz_snowf(k)*dtime         
           scv(k) = scv(k) + pg_snow(k)*dtime      ! snow water equivalent (mm)
           
           if(ivt(k)==11 .AND. tg(k)>tfrz)then
              scv(k)=0.; snowdp(k)=0.; sag(k)=0.
           endif

         endif
      enddo
!      print*,' Canopy interception and precipitation onto ground surface'
! ======================================================================
! [3] When the snow accumulation exceeds 10 mm, initialize a snow layer
! ======================================================================
      do k = 1, kpt
         if(ivt(k) == 17) cycle
         if(snl(k) == 0 .AND. pg_snow(k) > 0.0 &
                              .AND. snowdp(k) >= 0.01)then  
            snl(k) = -1
            newnode(k) = 1
            dz(0,k)  = snowdp(k)             ! meter
            z (0,k)  = -0.5*dz(0,k)
            zi(-1,k) = -dz(0,k)
! Currently, the water temperature for the precipitation is simply set 
! as the surface air temperature
            sag(k) = 0.                      ! snow age
            tss (0,k) = min(tfrz, tm(k))     ! K
            wice(0,k) = scv(k)               ! kg/m2
            wliq(0,k) = 0.                   ! kg/m2
            fiold(0,k) = 1.
!**         write(6,*) 'snow layer is built'
         endif

! The change of ice partial density of surface node due to precipitation
! only ice part of snowfall is added here, the liquid part will be added latter
         if(snl(k) < 0 .AND. newnode(k) == 0)then
            wice(snl(k)+1,k) = wice(snl(k)+1,k)+dtime*pg_snow(k)
            dz(snl(k)+1,k) = dz(snl(k)+1,k)+dz_snowf(k)*dtime
         endif
      enddo
!    print*,' When the snow accumulation exceeds 10 mm, initialize a snow layer'

! ======================================================================
! [4] Energy AND Water balance 
!    (please take care the dummy arrays as assumed-shape arrays, 
!     for expample, z(2:5,1) is a 1D array with 4 elements;
!     z(2:5,1:1) is a 4x1 2D array)
! ======================================================================

      do k = 1, kpt
         if(ivt(k) /= 17)then                ! NOT lake
            lb  = snl(k) + 1           ! lower bound of array 
            iub = msl                  ! upper bound of array

         call thermal ( ivt(k)      ,lb         ,iub       ,dtime      ,&
              csol(1,k) ,porsl(1,k) ,phi0(1,k)  ,bsw(1,k)  ,dkmg(1,k)  ,&  
              dkdry(1,k),dksatu(1,k),lai(k)     ,sai(k)    , z0m(k)    ,&
              displa(k) ,sqrtdi(k)  ,rootfr(1,k),effcon(k) ,vmax25(k)  ,&
              slti(k)   ,hlti(k)    ,shti(k)    ,hhti(k)   ,trda(k)    ,&
              trdm(k)   ,trop(k)    ,gradm(k)   ,binter(k) ,extkn(k)   ,&
              hu(k)     ,ht(k)      ,hq(k)      ,us(k)     ,vs(k)      ,&
              tm(k)     ,qm(k)      ,rhoair(k)  ,psrf(k)   ,pco2m(k)   ,&
              po2m(k)   ,cosz(k)    ,parsun(k)  ,parsha(k) ,sabvsun(k) ,&
              sabvsha(k),sabg(k)    ,frl(k)     ,extkb(k)  ,extkd(k)   ,&
              thermk(k) ,fsno(k)    ,sigf(k)    ,dz(lb,k)  ,z(lb,k)    ,&
              zi(lb-1,k),tlsun(k)   ,tlsha(k)   ,tss(lb,k) ,wice(lb,k) ,&
              wliq(lb,k),ldew(k)    ,scv(k)     ,snowdp(k) ,imelt(lb,k),&
              taux(k)   ,tauy(k)    ,fsena(k)   ,fevpa(k)  ,lfevpa(k)  ,&
              fsenl(k)  ,fevpl(k)   ,etr(k)     ,fseng(k)  ,fevpg(k)   ,&
              olrg(k)   ,fgrnd(k)   ,etrc(k)    ,qseva(k)  ,qsdew(k)   ,&
              qsubl(k)  ,qfros(k)   ,sm(k)      ,tref(k)   ,trad(k)    ,&
              rst(k)    ,assim(k)   ,respc(k)   ,errore(k) ,solisb(k) ,solisd(k))
!         print*,'thermal'
         call water (    ivt(k)     ,lb         ,iub       ,dtime      ,&
              z(lb ,k)  ,dz(lb,k)   ,zi(lb-1,k) ,bsw(1,k)  ,porsl(1,k) ,&
              phi0(1,k) ,hksati(1,k),rootfr(1,k),rootr(1,k),tss(lb,k)  ,&
              wliq(lb,k),wice(lb,k) ,pg_rain(k) ,sm(k)     ,etr(k)     ,&
              qseva(k)  ,qsdew(k)   ,qsubl(k)   ,qfros(k)  ,etrc(k)    ,&
              rsur(k)   ,rnof(k)    ,fevpg(k)   ,fevpa(k)  ,lfevpa(k)  )

         else if(ivt(k) == 17)then      ! lake 
         call lake ( 10 ,istep      ,1          ,dlat(k)   ,dtime      ,&
              z(1 ,k)   ,dz(1,k)    ,zi(0,k)    ,hu(k)     ,ht(k)      ,&
              hq(k)     ,us(k)      ,vs(k)      ,tm(k)     ,qm(k)      ,&
              snowrate(k),rhoair(k) ,psrf(k)    ,sabg(k)   ,frl(k)     ,&
              tg(k)     ,tss(1,k)   ,wliq(1,k)  ,wice(1,k) ,scv(k)     ,&
              snowdp(k) ,trad(k)    ,tref(k)    ,taux(k)   ,tauy(k)    ,&
              fsena(k)  ,fevpa(k)   ,lfevpa(k)  ,fseng(k)  ,fevpg(k)   ,&
              olrg(k)   ,fgrnd(k)  )
   
              lb = maxsnl+1
              iub = msl
              snl(k) = 0
              dz  (lb:0,k) = (/(0.,i=lb,0)/)
              z   (lb:0,k) = (/(0.,i=lb,0)/)
              zi(lb-1:0,k) = (/(0.,i=lb-1,0)/)
              tss (lb:0,k) = (/(0.,i=lb,0)/)
              wliq(lb:0,k) = (/(0.,i=lb,0)/)
              wice(lb:0,k) = (/(0.,i=lb,0)/) 

              tlsun(k)  = 0.
              tlsha(k)  = 0.
              ssw(k)   = 0.
              ldew(k)  = 0.
              etrc(k)  = 0.
              fsenl(k) = 0.
              fevpl(k) = 0.
              etr(k)   = 0.
              rsur(k)  = 0.
              rnof(k)  = 0.
              rst(k)   = 0.
              assim(k) = 0.
              respc(k) = 0.
              rootr(1:iub,k) = (/(0.,i= 1,iub)/)
              errore(k) = 0.

         endif

      enddo
!      print*,'Energy AND Water balance'
! ======================================================================
! [5] Compaction rate for snow 
!     Combine or divide thin or thick snow elements
! ======================================================================

      do k = 1, kpt
         if(ivt(k) == 17) cycle

         if(snl(k) < 0)then

! Natural compaction and metamorphosis. The compaction rate
! is recalculated for every new timestep

            lb  = snl(k) + 1   ! lower bound of array 
            call compact ( lb, dtime, &
                           imelt(lb:0,k),fiold(lb:0,k),tss(lb:0,k), &
                           wliq(lb:0,k) ,wice(lb:0,k) ,dz(lb:0,k) )

! Combine thin snow elements

            lb = maxsnl + 1
            call combin ( lb, snl(k),&
                          z(lb:1,k)   ,dz(lb:1,k)  ,zi(lb-1:1,k),&
                          wliq(lb:1,k),wice(lb:1,k),tss(lb:1,k),&
                          scv(k)      ,snowdp(k) )

! Divide thick snow elements

            if(snl(k) < 0) &
            call subdiv ( lb, snl(k), &
                          z(lb:0,k), dz(lb:0,k), zi(lb-1:0,k),  &
                          wliq(lb:0,k), wice(lb:0,k), tss(lb:0,k) )

! Set zero to the empty node

            if(snl(k) > maxsnl)then
               sag(k) = 0.
               do i = maxsnl+1, snl(k)
                  wice(i,k) = 0.
                  wliq(i,k) = 0.
                  tss(i,k) = 0.
                  dz(i,k) = 0.
               enddo
            endif
 
         endif
      enddo

! ======================================================================
! [6] Update the snow age 
! ======================================================================
!
      do k = 1, kpt
         if(ivt(k)/=17)then
            lb = snl(k) + 1
            tg(k) = tss(lb,k)
            ssw(k) = min(1.,1.e-3*wliq(1,k)/dz(1,k))
         endif
         call snowage (dtime, tg(k), scv(k), scvold(k), sag(k))
      enddo

      END SUBROUTINE clm_main