enkf_analysis.f90

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

   subroutine enkf_analysis(istep)

!  ----------------------------------------------------------------------------
! | This subroutine is used for updating the ensemble of state variables and   |
! | getting analysis after assimilation of remote sensing data. In the current |
! | version, only two bands microwave observations are considered and          |
! | snow-covered situation is not considered too.                              |
! |                                                                            |
! | source file: enkf_analysis.f90                                             |
! | author : QinJun, 6/3/2004                                                  |
!  ----------------------------------------------------------------------------

         use enkfpar_module ! many other parameters must be added into it

	 use clmtv_module

	 use clmtc_module

	 use enkftv_module

         use clmpar_module 

         use enkfobs_module

	 implicit none

	 real,dimension(1:Ne) ::sim_obs ! used to store simulated observations,&                                        ! MODIS-TEMP 
	                               ! from model state variables for each kpt
	 real    ::ave_obs ! the average of simulated observations for each kpt

	 real,dimension(1:Ne) ::per_obs ! sat_obs plus perturbations &
                                        ! with distribtuon N(0,Q)

	 real,dimension(1:10) ::PfHT    ! cross-covariance matrix between&
                                        ! the ensemble
                           ! of state variables and the ensemble of  
	                   ! simulated observations before
			   ! assimilation of remote sensing data for each kpt

	 real ::HPfHT  ! covariance matrix of the&
                       ! ensemble of 
	               ! simulated observations for each kpt

	 real,dimension(1:10) ::GAIN    ! kalman filter gain

	 real ::SQRTQ   ! square root of matrix Q

	 real ::RAND    ! multivariate random vector 1 by 1

	 real ::COVQ    ! multivariate random vector 1 by 1 with distribtion
	                ! N(0,Q)

	 real,dimension(1:10)  ::OFFSET  ! temporary use

         real::gasdev

	 integer::i,j,k,m,n

         integer::istep

!    call enkf_judge( )! judge whether or not observations exist 
     if(sat_obs.EQ.-999.0) then

       OBS = .FALSE.
     else

       OBS = .TRUE.
     end if

print*,'obs',obs




    if( OBS .EQV. .FALSE. ) then

       do i = 1, kpt

          tss (:,i) = 0.0


          do j = 1, Ne

                   do m = 1, 10

                      tss (m,i) = tss (m,i) + en_tss (m,i,j)

                   end do

          end do

          do m = 1, 10

                   tss (m,i) = tss (m,i)/Ne

          end do

       tg(i)=tss(1,i)         

       end do
!print*,'tss', tss
!print*,'wliq',wliq
!print*,'wice',wice
!print*,'tlsun',tlsun
!print*,'ldew',ldew
       

    else


!-------------------------------------------------------------------------
       do i = 1, kpt
	     
! ------------------------------------------------------------------------
! ** generate the ensemble of microwave observations at two frequencies **
! ------------------------------------------------------------------------
print*, ' start observation operator'
         do j = 1, Ne

           sim_obs(j)=reg_one+reg_two*en_tss(1,i,j)
         end do 
print*,'sim_obs',sim_obs

! --------------------------------------------------------------------
! ** evaluate the average of the ensemble of simulated modis temp 
! --------------------------------------------------------------------
if(Ne .GT. 0) then
        ave_obs = 0.0

        do j = 1, Ne

          ave_obs = ave_obs + sim_obs(j)

	end do

	  ave_obs = ave_obs/Ne

!print*,sim_obs(:,1)
!print*,sim_obs(:,2)
!print*,sim_obs(:,3)
!print*,sim_obs(:,4)
!print*,' ave_obs',ave_obs
! --------------------
! ** evaluate HPfHT **
! --------------------

        HPfHT = 0.0

	do j = 1, Ne


           HPfHT = HPfHT + (sim_obs(j)&
                 - ave_obs)*(sim_obs(j)- ave_obs) 
				 

        end do
print*,HPfHT

                      HPfHT = HPfHT/(Ne-1)

print*,HPfHT
!  -----------------------------------------------------------------------------
! | evaluate the average of the ensemble of state variables.                    |
! | firstly , calculate the averaged temperature of soil (tss )                 |
! |                                                                             |
! | the averages of the ensemble of state variables are stored in the original  |
! | clmtv_module.                                                               |
!  -----------------------------------------------------------------------------
!print*, tss
                tss (:,i) = 0.0

               do j = 1, Ne

                  do m=1,10
		      tss (m,i) = tss (m,i) + en_tss (m,i,j)
                  end do

	       end do

	       do m = 1, 10
		   tss (m,i) = tss (m,i)/Ne
	       end do

!  ---------------------
!  **  evaluate PfHT  **
!  ---------------------

        PfHT(:) = 0.0
	 
        do j = 1, Ne

	   do m = 1, 10

                  PfHT(m) = PfHT(m) + (en_tss(m,i,j)&
                          - tss(m,i))*(sim_obs(j)- ave_obs)


	   end do 


	end do 
print*,PfHT
!----------------------------------------------------------------------------

        do m = 1, 10 

		      PfHT(m) = PfHT(m)/(Ne-1)


	end do
print*,'PfHT',PfHT
! ------------------------------------
! ** compute the kalman filter gain **
! ------------------------------------

        GAIN(:) = 0.0

        call MATADD(HPfHT,Q,1,1)

        call MATINV(HPfHT,1)
print*,'HPfHT',HPfHT
        call MATMULT(PfHT,HPfHT,GAIN,10,1,1)

!GAIN(:,:) = 0.0        
!stop
print*,'GAIN',GAIN
!if(obs.EQV..true.) stop
! -------------------------------------------------------------
! ** generate the ensemble of perturbing remote sensing data **
! -------------------------------------------------------------

                per_obs(:) = 0.0 
print*,'Q',Q,'SQRTQ',SQRTQ
!                call PDSQRT(Q,SQRTQ,1)
print*,'Q',Q,'SQRTQ',SQRTQ
                 SQRTQ=sqrt(Q)
		do j = 1, Ne

		    RAND = gasdev(idum)
!print*,'RAND',RAND
		    call MATMULT(SQRTQ,RAND,COVQ,1,1,1)
!print*,'COVQ',COVQ

!sat_obs(4,i) = 215.20249283636920267782



		      per_obs(j) = sat_obs + COVQ

	           end do
!print*,'per_obs',per_obs
print*,'perturbing'
print*, ' sat_obs',sat_obs
!  --------------------------------------------
!  ** update the ensemble of state varialbes **
!  --------------------------------------------

        do j = 1, Ne

           OFFSET(:) = 0.0 
!print*,'update the ensemble of state varialbes'
           call MATSUB(per_obs(j),sim_obs(j),1,1)
!print*,'MATSUB', per_obs
	   call MATMULT(GAIN,per_obs(j),OFFSET,10,1,1)

                
print*,'OFFSET',OFFSET
	   do m = 1, 10

	      en_tss(m,i,j) = en_tss(m,i,j) + OFFSET(m)

	   end do

		   en_tg(i,j) = en_tss(1,i,j)    ! in CLM, tg(?) = tss(1,?)
		   
	end do

! -----------------------------
! **  calculate the analysis ** 
! -----------------------------
!print*,'calculate the analysis'
        OFFSET(:) = 0.0 

        call MATSUB(sat_obs,ave_obs,1,1)

	call MATMULT(GAIN,sat_obs,OFFSET,10,1,1)

	do m = 1, 10

	    tss(m,i) = tss(m,i) + OFFSET(m)

	end do


		   
        tg(i) = tss(1,i)
  
    end if


      end do
print*,'end enkf_analysis'

end if

   write(8,*) tg

   end subroutine enkf_analysis