enkf_drv.f90

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

      SUBROUTINE enkf_drv (lon ,lat ,istep ,lun_rst ,rdlai) 

!  ------------------------------------------------------------------------
! |  most procedures for EnKF are parasitized on the original CLM code.    |
! |  This clm_ensemble_drv procedure used to propogate ensemble state      |
! |  varialbes and their derived time-varying parameters(or varaibles)     |
! |  in time is modified from clm_drv,so if want to know physical meaning  |
! |  of each variables(or parameters), please read the procedure clm_drv.  |
! |                                                                        |
! |                                                                        |
! |                                                                        |
! |                                                                        |
! |    source file: enkf_drv.f90                                           |
! |  first version: Qin Jun, Dai yongjiu 5/27/2004                         |
! |  modified version: Qin Jun  8/19/2004                                  |
!  ------------------------------------------------------------------------

!=======================================================================
!      Source file: clm_drv.f90
! Original version: Yongjiu Dai, September 15, 1999
!
!=======================================================================
      USE clm2d_module
      USE phycon_module
      USE enkfpar_module
      use clmpar_module
      USE clmtc_module
      USE clmtv_module
      USE miscell_module
      USE enkftv_module

      IMPLICIT NONE
! ------------------- input variables -----------------------------
   integer, INTENT(in) :: &
      lon,                &! atm number of longitudes
      lat,                &! atm number of latitudes
      istep,              &! atm time step index
      lun_rst             ! logical unit number of restart file

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

! ---------------------- local variables --------------------------
   integer i,j,k,l,m       ! loop/array indices
   real pgcm(kpt), calday
   real solisb(kpt), solisd(kpt)

! ====================iteration cycle number for ensemble==============

   integer::ii,jj,kk,mm
     
! ======================================================================
!  [1] map atmospheric fields to force clm: [lon] x [lat] grid ->
!     [lpt] vector of land points -> [kpt] vector of subgrid points
! ======================================================================

      do k = 1, kpt                             !clm vector index
         l = klnd(k)                            !land point index
         i = ixy(l)                             !longitude index
         j = jxy(l)                             !latitude index
         pco2m(k) = pco2xy(i,j)
         po2m(k)  = po2xy(i,j)
         tm(k)    = tgcmxy(i,j)
         us(k)    = ugcmxy(i,j)
         vs(k)    = vgcmxy(i,j)
         qm(k)    = qgcmxy(i,j)
         pgcm(k)  = pgcmxy(i,j)
         psrf(k)  = psrfxy(i,j)

         sols(k)  = solsxy(i,j)
         soll(k)  = sollxy(i,j)
         solsd(k) = solsdxy(i,j)
         solld(k) = solldxy(i,j)
         frl(k)   = flwdsxy(i,j)

! Next sets upper limit of air temperature for snowfall at 275.65K.
! This cut-off was selected based on Fig. 1, Plate 3-1, of Snow
! Hydrology (1956).
         if(prcxy(i,j)+prlxy(i,j) > 0.) then
            if(tm(k) > tfrz + tcrit)then
              iptype(k) = 1
              rainrate(k) = prcxy(i,j)+prlxy(i,j)
              snowrate(k) = 0.
            else
              iptype(k) = 2
              rainrate(k) = 0.
              snowrate(k) = prcxy(i,j)+prlxy(i,j)
            endif
         else
           iptype(k) = 0
           rainrate(k) = 0.
           snowrate(k) = 0
         endif

         hu(k) = zgcmxy(i,j,1)
         ht(k) = zgcmxy(i,j,2)
         hq(k) = zgcmxy(i,j,3)

         rhoair(k) = (pgcm(k)-0.378*qm(k)*pgcm(k)/(0.622+0.378*qm(k))) &
                   / (rgas*tm(k))
      end do
!==============================================================================
!==============================================================================
!============================begin data assimilation code======================

! In the following, three subroutines are used to realize data assimilation
! 1) EnKF_ENSEMBLE_ASSIGN(), 2)EnKF_ANALYSIS(), 3)EnKF_ADD_NOISE(). AS for their
! detailed meanings and functionalities, read notes after their entities. Here,
! the process for data assimilation is outlined schematically and all details are
! omitted.
! 
!
!do kk=1,Ne
!    call EnKF_ENSEMBLE_ASSIGN()---> parasitize an ensemble member on CLM stem for 
!                                      propogating an ensemble state in time
!    call surrad()              ---> calculate solar radiation aborded by vegetation
!                                      and ground
!    call clm_main()            ---> updata an ensemble member in time 
!
!end do    
!
!    call EnKF_ANALYSIS()       ---> assimilate remote sensing data into CLM  
!
!    call EnKF_ADD_NOISE()      ---> add model noise into ensemble state variables
!                                      and these noisy ensemble state variables will be
!                                      used to produce derived parameters(diagnostic parameters)
!                                      for the next computation.
!
!do kk=1,Ne
!   call EnKF_ENSEMBLE_ASSIGN() ---> parasitize an ensemble member on CLM stem in order to
!                                      generate derived parameters(diagnostic parameters)
!                                      for the next calculation.
!   call ticktime()             ---> calendar information for the next time step for each
!                                      ensemble member
!   call laiini()               ---> ecosystem dynamics: phenology, vegetation, soil carbon
!                                      for each ensemble member
!   call fractsnow()            ---> fraction of snow cover for each ensemble member
!
!   call clmzen()               ---> compute the cosine of solar zenith angle for next-time step
!                                      for each ensemble member
!   call albedo()               ---> compute albedo (used in the next step radiation computations)
!                                      for each ensemble member
!end do    
!=========================================================================================

     
!    ----------------------------------------------------------------------
!   |  start updating ensemble members                                     |
!    ----------------------------------------------------------------------
	  
       do kk=1, Ne


! Initial set for water and Energy balance check

         do k = 1, kpt
            if(ivt(k)==17)then
               totwb(k) = 0.
               cycle            
	    endif

            totwb(k) = en_ldew(k,kk) + en_scv(k,kk)
            do i = 1, msl
               totwb(k) = totwb(k) + en_wice(i,k,kk) + en_wliq(i,k,kk)
            enddo
         enddo

! ======================================================================
!  [2] Solar absorbed by vegetation and ground
! ======================================================================

     call surrad(kpt ,en_sigf(:,kk) ,en_albg(:,:,:,kk) ,en_albv(:,:,:,kk)  ,&
        en_alb(:,:,:,kk) ,en_ssun(:,:,:,kk), en_ssha(:,:,:,kk) ,sols ,soll ,&
                solsd ,solld ,parsun, parsha, sabvsun, sabvsha, sabg ,sabvg,&
                 solisb, solisd )


! ======================================================================
!  [3] Main driver for CLM
!      Note: the modification should do here if pursue the multitask processes 
!      "process the "big" vector of [kpt] points as [numlv] "little"
!      vectors of [numpnt] points. [beg] is the starting location of the
!      [numpnt] points for the "little" vector in the "big" [kpt] vector.
!      multitask [numlv] processes"
! ======================================================================

print*, 'before clm_main',kk

         call clm_main ( kpt, ivt, msl, maxsnl, istep, dtime, dlat,  &

                  ! soil information
                    csol, porsl, phi0, bsw,                          &
                    dkmg, dksatu, dkdry,                             &
                    hksati,                                          &

                  ! vegetation information
                    z0m, displa, sqrtdi,                             &
                    effcon, vmax25, slti, hlti, shti, hhti,          &
                    trda, trdm, trop, gradm, binter, extkn,          &
                    en_lai(:,kk), en_sai(:,kk), rootfr,              &

                  ! atmospheric forcing
                    pco2m, po2m, us, vs, tm, qm, psrf, rhoair,       &
                    iptype, rainrate, snowrate,                      &
                    hu, ht, hq,                                      &

                  ! model variables needed for restart
                    en_snl(:,kk), en_dz(:,:,kk), en_z(:,:,kk), en_zi(:,:,kk),  &
                    en_tss(:,:,kk), en_tlsun(:,kk), en_tlsha(:,kk),            &
		    en_ldew(:,kk), en_wliq(:,:,kk),  en_wice(:,:,kk),          &
                    en_sag(:,kk), en_scv(:,kk),  en_snowdp(:,kk),              &
         	    en_rootr(:,:,kk), en_etrc(:,kk),   en_tg(:,kk),            &
                                                                  

                  ! fraction
                    en_fveg(:,kk),wt,en_fsno(:,kk),en_sigf(:,kk),ssw,          &

                  ! fluxes
                    taux, tauy,                                      &
                    fsena, fevpa, lfevpa, fsenl, fevpl, etr,         &
                    fseng, fevpg, olrg, fgrnd, trad, tref,           &
                    rsur, rnof, rst, assim, respc,                   &

                  ! net solar radiation
                   en_cosz(:,kk),frl,parsun,parsha,sabvsun,sabvsha,sabg,sabvg, &
                   en_extkb(:,kk), en_extkd(:,kk), en_thermk(:,kk), solisb, &
	           solisd,                                                  &

                  ! energy balance error
                    errore )
!             print*,'clm_main'
! ======================================================================
! [4] Water and Energy balance check
! ======================================================================

         call balance_a (kpt ,ivt, msl ,dtime ,en_ldew(:,kk) ,en_scv(:,kk) ,&
                      en_wice(1:msl,1:kpt,kk) ,en_wliq(1:msl,1:kpt,kk) ,    &
                       rainrate ,snowrate ,fevpa  ,rnof, totwb ,errore ,    &
                        en_xerr(:,kk) ,en_zerr(:,kk))

      end do

!    ----------------------------------------------------------------------
!   |  end updating ensemble members                                       |
!    ----------------------------------------------------------------------

!============================================================================
     print*,'terminate clm_main'
!    ----------------------------------------------------------------------
!   |  start the analysis part and adding-model-noise part                 |
!    ----------------------------------------------------------------------
print*,'istep',istep 
       
 
        call enkf_analysis(istep)
       print*,'enkf_analysis()'

        call enkf_add_noise(istep)
       print*,'enkf_add_noise()'



!    ----------------------------------------------------------------------
!   |  end the analysis part and adding-model-noise part                   |
!    ----------------------------------------------------------------------

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

!    ----------------------------------------------------------------------
!   |  start generating the time-varying parameters needed in the next     |
!   |  propogation of the ensemble members                                 |
!   |                                                                      |
!   |                                                                      |
!    ----------------------------------------------------------------------

      do kk=1, Ne


! =============================================================================
! [5] Calendar information for the next time step
! =============================================================================

         call ticktime (dtime ,en_jday(kk) ,en_iyear(kk) ,en_mcsec(kk))

! ==============================================================================
! [6] Ecosystem dynamics: phenology, vegetation, soil carbon
!     or call satellite data to get LAI,.., ect.;
!     Fraction of snow cover.
! =============================================================================
    
         calday = float(en_jday(kk)) + float(en_mcsec(kk))/86400.
         call laiini(kpt,ivt,dlat,calday,glai,gsai,en_green(:,kk),&
		     en_fveg(:,kk),en_lai(:,kk),en_sai(:,kk))
         call fractsnow (kpt,en_fveg(:,kk),z0m,en_snowdp(:,kk),wt,&
		         en_sigf(:,kk),en_fsno(:,kk))

! ==============================================================================
! [7] Compute the cosine of solar zenith angle for next-time step
! =============================================================================
     
         call clmzen (kpt, calday, dlon, dlat, en_cosz(:,kk))

! ==============================================================================
! [8] Compute albedo (used in the next step radiation computations)
! =============================================================================

         call albedo (kpt,ivt,albsol,albvgs,albvgl,&
                  chil,ref,tran,en_fveg(:,kk),en_green(:,kk),en_lai(:,kk), &
	      en_sai(:,kk),en_cosz(:,kk),wt,en_fsno(:,kk),ssw,en_tg(:,kk), &
             en_scv(:,kk),en_sag(:,kk),en_alb(:,:,:,kk),en_albg(:,:,:,kk), &
	            en_albv(:,:,:,kk),en_ssun(:,:,:,kk),en_ssha(:,:,:,kk), &
        en_tranc(:,:,:,kk),en_thermk(:,kk),en_extkb(:,kk),en_extkd(:,kk))

! ======================================================================
! [9] Return required surface fields to atmospheric model:
!     [kpt] vector of subgrid points -> [lpt] vector of land points ->
!     [lon] x [lat] grid. this mapping is for land points only.
!     non-land points are undefined.
! ======================================================================

      call clm2d_return (lon ,lat)

! ======================================================================
! [10] Write out the model variables for restart run 
! ======================================================================

!      call restartsave (lun_rst)



!    ----------------------------------------------------------------------
!   |  end generating the time-varying parameters needed in the next       |
!   |  propogation of the ensemble members                                 |
!   |                                                                      |
!   |                                                                      |
!    ----------------------------------------------------------------------

      end do

      END SUBROUTINE enkf_drv