assim_nature_isba_enkf.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
SUBROUTINE assim_nature_isba_enkf(I, HPROGRAM, KI, PT2M, PHU2M, HTEST)

! -----------------------------------------------------------------------------
!
! Land Data Assimilation System based on an Extended Kalman Filter
!
! Revised version : JFM (15 September 2008)
!
! The control vector can be any element of (TG1,TG2,WG1,WG2) - Choice in namelist
!
! The observations can be any element of (T2M,HU2M,WG1) - Choice in namelist
!
! Possibility to evolve the B matrix in the cycling - otherwise SEKF
!
! First version including patches (15 October 2008)
! Trygve Aspelien, Separating IO  06/2013
! Alina Barbu: bug correction of B matrix, otherwise understimation of the gain matrix (11/2013) 
! Alina Barbu: equivalent analysis of B matrix to ensure its symetric and positiv definiteness properties (11/2013) 
  
! -----------------------------------------------------------------------------
!
USE modd_surfex_mpi,    ONLY : nrank, npio
USE modd_assim,         ONLY : nobstype, xerrobs, nvar, nprintlev, cvar, &
                               xf_patch, xf,cobs,cfile_format_obs,nens, &
                               nechgu, nboutput, nobs, xyo, lenkf, ldenkf, &
                               lpb_correlations, lperturbation_run, &
                               lbias_correction, xinfl
! 
USE modd_surf_par,      ONLY : xundef
!
USE modd_isba_n,        ONLY : isba_t
!
#ifdef SFX_ARO
USE yommp0,             ONLY : myproc 
#endif
!
USE yomhook,            ONLY : lhook,dr_hook
USE parkind1,           ONLY : jprb
!
USE modi_abor1_sfx
USE modi_add_forecast_to_date_surf
!
USE modi_outer_product
USE mode_ekf
USE mode_random
!
! -----------------------------------------------------------
!
IMPLICIT NONE
!
TYPE(isba_t), INTENT(INOUT) :: i
!
 CHARACTER(LEN=6),   INTENT(IN) :: hprogram     ! program calling surf. schemes
INTEGER,            INTENT(IN) :: ki
REAL, DIMENSION(:), INTENT(IN) :: pt2m
REAL, DIMENSION(:), INTENT(IN) :: phu2m
 CHARACTER(LEN=2),   INTENT(IN) :: htest        ! must be equal to 'OK'
!
!    Declarations of local variables
!
 CHARACTER(LEN=30)  :: ybgfile
 CHARACTER(LEN=17)  :: ylfname
 CHARACTER(LEN=9)   :: yfname
 CHARACTER(LEN=7)   :: ymyproc
 CHARACTER(LEN=1)   :: ychar
!
! Local Matrix for Analysis calculation
!
!  Allocation
!  Perturbed simulations
!
! Initial values (to be analysed)
! Observations
!
! Temporary vectors used by the EKF approach
REAL,DIMENSION(KI) :: zcofswi                     ! dynamic range (Wfc - Wwilt)
!
REAL,DIMENSION(NENS+1,NOBSTYPE*NBOUTPUT) :: zyf            ! Vector of model observations (averaged)
!
REAL,DIMENSION(NOBSTYPE*NBOUTPUT,NENS) :: zinnov
!
REAL,DIMENSION(NVAR,NOBSTYPE*NBOUTPUT) :: zbht
REAL,DIMENSION(NOBSTYPE*NBOUTPUT,NOBSTYPE*NBOUTPUT) :: zr        ! covariance matrix of observation errors
REAL,DIMENSION(NOBSTYPE*NBOUTPUT,NOBSTYPE*NBOUTPUT) :: zk1, zk2, zhbht
REAL,DIMENSION(NOBSTYPE*NBOUTPUT) :: zp, zx
!
REAL, DIMENSION(I%NPATCH,NVAR,NENS) :: zxincr
REAL, DIMENSION(I%NPATCH,NVAR,NENS) :: za, zf
REAL, DIMENSION(I%NPATCH,NOBSTYPE*NBOUTPUT,NENS) :: zf_patch
REAL, DIMENSION(I%NPATCH,NVAR) :: zf_mean, za_mean
REAL, DIMENSION(NVAR) :: zf_agg, za_agg
!
REAL,DIMENSION(:,:,:),ALLOCATABLE :: zf_mean0, zf_patch_mean
!
REAL :: ztime, zalpha                    ! current time since start of the run (s)

INTEGER :: iobscount
INTEGER :: iyear                      ! current year (UTC)
INTEGER :: imonth                     ! current month (UTC)
INTEGER :: iday                       ! current day (UTC)
INTEGER :: ihour
INTEGER :: iresp                      ! return code
INTEGER :: istep                      ! 
INTEGER :: imyproc
INTEGER :: iobs, iens
INTEGER :: istat, icpt, iunit
!
INTEGER :: ii,j,k,jj,l,k1,l1
!
LOGICAL :: gbexists
!
REAL(KIND=JPRB)                            :: zhook_handle
!
!
IF (lhook) CALL dr_hook('ASSIM_NATURE_ISBA_ENKF',0,zhook_handle)

#ifdef USE_SODA

!
!############################# BEGINNING ###############################
!
IF (htest/='OK') THEN
  CALL abor1_sfx('ASSIM_NATURE_ISBA_ENKF: FATAL ERROR DURING ARGUMENT TRANSFER')
END IF
!
IF ( nprintlev>0 .AND. nrank==npio) THEN
  WRITE(*,*)
  WRITE(*,*) '   --------------------------'
  WRITE(*,*) '   |   ENTERING  VARASSIM   |'
  WRITE(*,*) '   --------------------------'
  WRITE(*,*)
ENDIF
!
#ifdef SFX_ARO
IF ( myproc > 0 ) THEN 
  imyproc = myproc
ELSE
  imyproc = 1
ENDIF
#else
imyproc = nrank+1
#endif
!
WRITE(ymyproc(1:7),'(I7.7)') imyproc
!
IF ( nprintlev > 0 .AND. nrank==npio ) WRITE(*,*) 'number of patches =',i%NPATCH
!
!############################# INITIALISATIONS ###############################
!
!   Read CLAY fraction to  compute the SWI range (Wfc - Wwilt)
!   (XSIGMA is defined in terms of SWI), need to convert to equivalent v/v
!   using same clay fraction in both layers
!   Read SAND fraction to compute the saturation for conversion of ERS SWI
!
 CALL cofswi(i%XCLAY(:,1),zcofswi)
!DO I=1,KI
  !ZSMSAT (I) = 0.001 * (-1.08*100.*XSAND(I,1) + 494.305)
  !ZWILT  (I) = 0.001 * 37.1342 * ((100.*XCLAY(I,1))**0.5) 
!ENDDO
!
! ====================================================================
!
! Analysis
!
! ====================================================================
!
!   Time reinitialization 
iyear  = i%TTIME%TDATE%YEAR
imonth = i%TTIME%TDATE%MONTH
iday   = i%TTIME%TDATE%DAY
!
ihour = 0
ztime = float(nechgu) * 3600.
!
!############################# READS OBSERVATIONS ###############################
!
! Map the variables in case we read them from CANARI inline/offline FA files
! At the moment only T2M and HU2M can be used. If other variables should be used 
! they must be added to the interface or be read from file.
IF ( trim(cfile_format_obs) == "FA" ) THEN
  !
  DO iobs = 1,nobstype
    SELECT CASE (trim(cobs(iobs)))   
      CASE("T2M") 
        xyo(:,iobs) = pt2m(:)
      CASE("HU2M")   
        xyo(:,iobs) = phu2m(:)
      CASE("WG1","LAI")  
        CALL abor1_sfx("Mapping of "//cobs(iobs)//" is not defined in ASSIM_NATURE_ISBA_EKF!")
    END SELECT                 
  ENDDO
  !  
ENDIF
!
!//////////////////////TO WRITE OBS/////////////////////////////////////
IF ( nprintlev > 0 ) OPEN (unit=111,file='OBSout.'//ymyproc,status='unknown',iostat=istat)
DO ii = 1,ki
  IF ( minval(i%XWGI(ii,1,:))>0. ) THEN
    xyo(ii,:) = xundef
    IF ( nprintlev > 1 ) WRITE(*,*) 'OBSERVATION FOR POINT ',ii,' REMOVED'
  ENDIF
  IF ( nprintlev > 0 ) WRITE (111,*) xyo(ii,:)
ENDDO
IF ( nprintlev > 0 ) CLOSE(111)
!//////////////////////TO WRITE OBS/////////////////////////////////////
!
! Recentering THE FORECAST ENSEMBLE MEMBERS
IF ( lbias_correction ) THEN
  !
  ALLOCATE(zf_mean0(ki,i%NPATCH,nvar))
  ALLOCATE(zf_patch_mean(ki,i%NPATCH,nobs))
  !
  DO ii = 1,ki
    DO j=1,i%NPATCH
      DO l = 1,nvar
        zf_mean0(ii,j,l) = sum(xf(ii,j,1:nens,l))/REAL(nens)
      ENDDO
      DO k = 1,nobs
        zf_patch_mean(ii,j,k) = sum(xf_patch(ii,j,1:nens,k))/REAL(nens)
      ENDDO
    ENDDO
  ENDDO
  !
  DO ii = 1,ki
    DO j = 1,i%NPATCH
      DO iens = 1,nens
        !
        DO l = 1,nvar
          IF ( xf(ii,j,iens,l) - zf_mean0(ii,j,l) + xf(ii,j,nens+1,l)>0.0 ) THEN
            xf(ii,j,iens,l) = xf(ii,j,iens,l) - zf_mean0(ii,j,l) + xf(ii,j,nens+1,l)
          ENDIF
        ENDDO
        !
        DO k = 1,nobs
          IF ( xf_patch(ii,j,iens,k) - zf_patch_mean(ii,j,k) + xf_patch(ii,j,nens+1,k)>0.0 ) THEN
            xf_patch(ii,j,iens,k) = xf_patch(ii,j,iens,k) - zf_patch_mean(ii,j,k) + xf_patch(ii,j,nens+1,k)
          ENDIF
        ENDDO
        !
      ENDDO
    ENDDO
  ENDDO
  !
  DEALLOCATE(zf_mean0,zf_patch_mean)
  !
ENDIF
!
!############################# ANALYSIS ###############################
!
IF ( nprintlev > 0 ) THEN
  IF (nrank==npio) THEN
    WRITE(*,*) 'calculating jacobians',nobs
    WRITE(*,*) ' and then PERFORMING ANALYSIS'
  ENDIF
  !
  !//////////////////////TO WRITE ANALYSIS ARRAYS/////////////////////////////////////
  ! WRITE OUT OBS AND YERROR FOR DIAGNOSTIC PURPOSES
  OPEN (unit=111,file='OBSERRORout.'//ymyproc,status='unknown',iostat=istat)
  ! *** Write innovations in ASCII file ***
  OPEN (unit=112,file='INNOV.'//ymyproc,status='unknown',iostat=istat)
  ! Write analysis results and increments in ASCII file
  OPEN (unit=113,file='ANAL_INCR.'//ymyproc,status='unknown',iostat=istat)
  ! **** Write out the observation operator + Gain matrix ****
  OPEN (unit=114,file='K1.'//ymyproc,status='unknown',iostat=istat)
  OPEN (unit=115,file='F_AGG.'//ymyproc,status='unknown',iostat=istat)
  OPEN (unit=116,file='A_AGG.'//ymyproc,status='unknown',iostat=istat)
ENDIF
!//////////////////////TO WRITE ANALYSIS ARRAYS/////////////////////////////////////
!
iobscount = 0
!
 CALL init_random_seed()
!
DO ii=1,ki
  !
  !---------------- MEAN SIMULATED OBS AVERAGED OVER TILES-----------------------
  zyf(:,:) = 0. 
  DO j=1,i%NPATCH
    IF (i%XPATCH(ii,j) > 0.0) THEN
      WHERE ( xf_patch(ii,j,:,:)/=xundef ) 
        zyf(:,:) = zyf(:,:) + i%XPATCH(ii,j)*xf_patch(ii,j,:,:)
      ENDWHERE
    ENDIF
  ENDDO
  IF ( nprintlev > 1 ) WRITE(*,*) 'read in sim obs yf', zyf(:,1)
  !
  !
  zr(:,:) = 0. ! Observation error matrix
  zinnov(:,:) = 0.
  !
  zxincr(:,:,:) = 0.0
  zk1(:,:) = 0.0
  zk2(:,:) = 0.0
  zp(:) = 0.0
  zx(:) = 0.0
  zbht(:,:) = 0.0
  zhbht(:,:) = 0.0
  
  DO istep=1,nboutput
    !
    DO k = 1,nobstype
      !
      k1 = (istep-1)*nobstype + k
      !
!--------------------- SET OBSERVATION ERROR ------------------      
      zr(k1,k1) = xerrobs(k)*xerrobs(k)
      IF ( cobs(k) .EQ. "LAI" ) THEN
        zr(k1,k1) = zr(k1,k1) * xyo(ii,k1)*xyo(ii,k1)
      ELSEIF (cobs(k) .EQ. "WG1") THEN
        ! convert R for wg1 from SWI  to abs value
        zr(k1,k1) = zr(k1,k1) * zcofswi(ii)*zcofswi(ii)
      ENDIF
      !
      ! Apply quality control
      IF ( abs(xyo(ii,k1)-zyf(nens+1,k1)) > 0.5 .OR. zr(k1,k1) < 0. ) xyo(ii,k1) = xundef
      !
      IF( xyo(ii,k1).NE.xundef .AND. xyo(ii,k1).NE.999.0 ) THEN         !if obs available
        !
        DO iens = 1,nens
          !
          zinnov(k1,iens) = xyo(ii,k1) - zyf(iens,k1)
          IF (lenkf) THEN
            zinnov(k1,iens) = zinnov(k1,iens) + random_normal() * (xerrobs(k)*zcofswi(ii))
          ENDIF
          !
        ENDDO
        !
        iobscount = iobscount + 1
        !
      ENDIF
      !
    ENDDO
    !
    IF ( nprintlev > 0 ) THEN
      WRITE(111,*) zr(:,:)
      DO k = 1,nobs
        WRITE(112,*) (sum(zinnov(k,:))/(nens*1.0)), xyo(ii,k)
      ENDDO
    ENDIF
    !
  ENDDO
  !
  IF (nprintlev > 0) WRITE(*,*) 'PERFORMING ANALYSIS'
  !
  !---------------******  SOIL ANALYSIS *******--------------------------
  !
  DO iens = 1,nens
    zf(:,:,iens) = xf(ii,:,iens,:)
    zf_patch(:,:,iens) = xf_patch(ii,:,iens,:)
  ENDDO

  !
  zf_mean(:,:) = sum(zf(:,:,:),dim=3)/REAL(nens)
  !
  IF (.NOT.lperturbation_run) THEN
    !
    DO j = 1,i%NPATCH
      !
      CALL outer_product(nens,nvar,nobs,zf(j,:,:),i%XPATCH(ii,j)*zf_patch(j,:,:),&
                         zbht(:,:),zhbht(:,:),lpb_correlations,cvar,cobs)
      !
      zk1(:,:) =  zhbht(:,:) + zr(:,:)
      zk2(:,:) = zk2(:,:) + zhbht(:,:) 
      CALL choldc(nobs,zk1(:,:),zp(:))
      !
      DO iens = 1,nens           ! Cholesky decomposition (1)
        !
        CALL cholsl(nobs,zk1(:,:),zp(:),zinnov(:,iens),zx(:))   ! Cholesky decomposition (2)       
        zxincr(j,:,iens) = matmul(zbht(:,:),zx(:))
        !
        DO l = 1,nvar
          za(j,l,iens) = zf(j,l,iens)
          IF (cvar(l)/="WG3" .AND. cvar(l)/="TG3") THEN
            za(j,l,iens) = za(j,l,iens) + zxincr(j,l,iens)
          ENDIF
        ENDDO
        !
      ENDDO
      !
    ENDDO
    !
    za_mean(:,:) = sum(za(:,:,:),dim=3)/REAL(nens)
    !
    IF (nprintlev>1) WRITE(114,*) (zk2(:,:) + zr(:,:))
    !
    DO iens = 1,nens
      !
      DO l = 1,nvar
        !
        DO j = 1,i%NPATCH       
          !
          IF ( ldenkf .AND. cvar(l)/="WG3" .AND. cvar(l)/="TG3" ) THEN
            !
            DO k = 1,nobs
              zalpha = 1.0 / ( 1.0 + sqrt( zr(k,k)/(zk2(k,k) + zr(k,k)) ) )
              za(j,l,iens) = za_mean(j,l) + (1.0-zalpha) * ( zf(j,l,iens)-zf_mean(j,l) ) &
                                          + zalpha * ( za(j,l,iens)-za_mean(j,l))
            ENDDO
            !
          ENDIF                 
          !apply inflation factor to the ensemble spread
          za(j,l,iens) = za_mean(j,l) + xinfl(l) * (za(j,l,iens) - za_mean(j,l))
          !
        ENDDO
        !
      ENDDO
      !
    ENDDO
    !
    WHERE (za(:,:,:)<0.0)
      za(:,:,:) = zf(:,:,:)
    END WHERE
    !
    !
  ELSE
    !
    za(:,:,:) = zf(:,:,:)
    !
  ENDIF
  !
  za_mean(:,:) = sum(za(:,:,:),dim=3)/REAL(nens)
  !
  zf_agg(:) = 0.
  za_agg(:) = 0.
  DO l = 1,nvar
    DO j = 1,i%NPATCH
      IF (za_mean(j,l)/=xundef .AND. zf_mean(j,l)/=xundef) THEN
        zf_agg(l) = zf_agg(l) + i%XPATCH(ii,j) * zf_mean(j,l)
        za_agg(l) = za_agg(l) + i%XPATCH(ii,j) * za_mean(j,l)
      ENDIF
    ENDDO
  ENDDO
  !
  IF (nprintlev>1) THEN
    WRITE(115,*) (zf_agg(l), l=1,nvar)
    WRITE(116,*) (za_agg(l), l=1,nvar)
  ENDIF
  !
  !############################# GET VARIABLES FOR OUTPUT WRITING ###############################
  !
  DO iens = 1,nens
    xf(ii,:,iens,:) = za(:,:,iens)
  ENDDO
  !
  IF (lbias_correction) xf(ii,:,nens+1,:) = za_mean(:,:)
  !
ENDDO
!
IF ( nprintlev > 0 ) THEN
  CLOSE(111)
  CLOSE(112)
  CLOSE(113)
  CLOSE(114)
  CLOSE(115)
  CLOSE(116)
ENDIF
!
IF (lhook) CALL dr_hook('ASSIM_NATURE_ISBA_ENKF',1,zhook_handle)

#endif

!
END SUBROUTINE assim_nature_isba_enkf