assim_nature_isba_ekf.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_ekf (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_isba_n, ONLY : isba_t
!
USE modd_surfex_mpi,    ONLY : nrank, npio
!
USE modd_assim,         ONLY : lbev, lbfixed, nobstype, xerrobs, xqcobs, nnco, nvar, nncv, &
                               xscale_q, nprintlev, cvar, xsigma, cbio, xi,        &
                               xf_patch, xf, cobs, xscale_qlai,cfile_format_obs,   &
                               xalph,nechgu, nboutput, xtprt, xlai_pass, xbio_pass,&
                               nobs, xyo
! 
USE modd_surf_par,      ONLY : xundef
!
!
#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 mode_ekf
!
! -----------------------------------------------------------
!
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=19)  :: 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(KI) :: ZSMSAT                      ! saturation  
!REAL,DIMENSION(KI) :: ZWILT
!
REAL,DIMENSION(KI,I%NPATCH,NVAR) :: zcoef
REAL,DIMENSION(KI,I%NPATCH,NVAR) :: zeps            ! The perturbation amplitude
!
REAL,DIMENSION(NVAR+1,NOBSTYPE) :: zyf            ! Vector of model observations (averaged) 
!
REAL,DIMENSION(KI*I%NPATCH*NVAR*I%NPATCH*NVAR) :: zblong
REAL,DIMENSION(KI,I%NPATCH*NVAR,I%NPATCH*NVAR) :: zb           ! background error covariance matrix
!
REAL,DIMENSION(I%NPATCH*NVAR,I%NPATCH*NVAR) :: zltm         ! linear tangent matrix for the f'ward model
REAL,DIMENSION(I%NPATCH*NVAR,I%NPATCH*NVAR) :: zq           ! model error matrix
!
REAL,DIMENSION(NOBSTYPE*NBOUTPUT,I%NPATCH*NVAR) :: zho, zhowr             ! Jacobian of observation operator
REAL,DIMENSION(I%NPATCH*NVAR,NOBSTYPE*NBOUTPUT) :: zhot            ! Transpose of HO
REAL,DIMENSION(I%NPATCH*NVAR,NOBSTYPE*NBOUTPUT) :: zgain           ! Kalman gain (used explicitly for Ba) 
!
REAL,DIMENSION(NOBSTYPE*NBOUTPUT,NOBSTYPE*NBOUTPUT) :: zr        ! covariance matrix of observation errors
REAL,DIMENSION(NOBSTYPE*NBOUTPUT,NOBSTYPE*NBOUTPUT) :: zk1
REAL,DIMENSION(NOBSTYPE*NBOUTPUT) :: zx,zb2,zp
!
REAL,DIMENSION(I%NPATCH*NVAR,I%NPATCH*NVAR) :: zkrk
REAL,DIMENSION(I%NPATCH*NVAR,I%NPATCH*NVAR) :: zidkh
REAL,DIMENSION(I%NPATCH*NVAR,I%NPATCH*NVAR) :: zident          ! identitiy matrix, used for Ba
!
REAL,DIMENSION(I%NPATCH) :: zvlaimin
!
REAL :: ztime                      ! 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
INTEGER :: istat, icpt, iunit
!
INTEGER :: ji,jj,jk,jjj,jl,k1,l1
!
LOGICAL :: gbexists
!
REAL(KIND=JPRB)                            :: zhook_handle
!
!
IF (lhook) CALL dr_hook('ASSIM_NATURE_ISBA_EKF',0,zhook_handle)

#ifdef USE_SODA

!
!############################# BEGINNING ###############################
!
IF (htest/='OK') THEN
  CALL abor1_sfx('ASSIM_NATURE_ISBA_EKF: 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 JI=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
!
! Set control variables
zident(:,:) = 0.                   ! identity matrix
DO jl = 1,nvar
  !
  DO jj = 1,i%NPATCH
    zident(jj+i%NPATCH*(jl-1),jj+i%NPATCH*(jl-1)) = 1.0
  ENDDO
  !
  WHERE ( xi(:,:,jl)/=xundef )
    zeps(:,:,jl) = xtprt(jl) * xi(:,:,jl)
  ELSEWHERE
    zeps(:,:,jl) = 1.
  ENDWHERE
  !
  IF ( trim(cvar(jl))=='WG2' .OR. trim(cvar(jl))=='WG1' ) THEN
    !
    DO ji = 1,ki
      zcoef(ji,:,jl) = zcofswi(ji)*zcofswi(ji)
    ENDDO
    !
  ELSEIF ( trim(cvar(jl))=='LAI' .AND. lbfixed ) THEN
    !
    DO ji = 1,ki
      DO jj = 1,i%NPATCH
        IF ( xlai_pass(ji,jj)/=xundef .AND. xlai_pass(ji,jj)>=2. ) THEN
          zcoef(ji,jj,jl) = xlai_pass(ji,jj)*xlai_pass(ji,jj)
        ELSE 
          zcoef(ji,jj,jl) = 0.4*0.4/(xsigma(jl)*xsigma(jl))
        ENDIF
      ENDDO
    ENDDO
    !
  ELSE
    !
    zcoef(:,:,jl) = 1.
    !
  ENDIF
  !
ENDDO
!
!############################# B CALCULATION ###############################
!
! ----------------------
! VARASSIM OPTION : LBEV
! ----------------------
!   Calculate the LTM, and evolve B. 
!
! Set the B input file depending of an existing B was found or not
ybgfile = "BGROUNDin."//trim(ymyproc)
INQUIRE (file=trim(ybgfile),exist=gbexists)
!
IF ( lbev .AND. gbexists ) THEN
  !
  zb(:,:,:) = 0.
  CALL b_big_loop(i, &
                  "READ",ybgfile,zb)
  IF ( nprintlev > 0 ) WRITE(*,*) 'read previous B matrix  ==>',zb(1,1,1),nvar
  !
ELSEIF ( lbev .OR. lbfixed ) THEN
  !
  ! Initialization of B 
  zb(:,:,:) = 0.
  DO ji = 1,ki
    DO jl = 1,nvar
      DO jj = 1,i%NPATCH   
        !
        l1 = jj + i%NPATCH *(jl-1)
        zb(ji,l1,l1) = xsigma(jl)*xsigma(jl) * zcoef(ji,jj,jl)
        !
      ENDDO
    ENDDO
    !
  ENDDO
  !
  IF ( lbev ) THEN
    !
    zblong(:) = 0.
    icpt = 0
    !
    DO ji = 1,ki
      DO jj = 1,i%NPATCH
        DO jjj = 1,i%NPATCH
          DO jl = 1,nvar
            DO jk = 1,nvar
              !
              l1 = jj + i%NPATCH * (jl-1)
              !
              icpt = icpt + 1
              zblong(icpt) = zb(ji,l1,l1)
              !
            ENDDO
          ENDDO
        ENDDO
      ENDDO
    ENDDO
    !
    zb(:,:,:) = 0.
    CALL b_big_loop(i, &
                  "BUIL","",zb,zblong)
    IF ( nprintlev > 0 ) WRITE(*,*) 'Initialized B'
    !
  ENDIF
  !
ELSE
  !
  CALL abor1_sfx("LBEV or LBFIXED should be .TRUE.!")
  !
ENDIF
!
IF ( lbev ) THEN
  !
!//////////////////////TO WRITE LTM/////////////////////////////////////
  IF (nprintlev>0) THEN
    iunit = 120
    DO jl=1,nvar
      DO jk=1,nvar
        iunit = iunit + 1
        WRITE(ychar,'(I1)') jk
        ylfname='LTM_del'//trim(cvar(jk))//'_del'//trim(cvar(jl))//"."//trim(ymyproc)
        OPEN(unit=iunit,file=ylfname,form='FORMATTED',status='UNKNOWN',position='APPEND')
      ENDDO
    ENDDO
  ENDIF
!/////////////////////TO WRITE LTM////////////////////////////////////////
  DO ji = 1,ki
    !
    ! calculate LTM
    zltm(:,:) = 0.0
    iunit = 120
    DO jl = 1,nvar    ! control variable (x at previous time step)
      DO jk = 1,nvar 
        iunit = iunit + 1
        DO jj = 1,i%NPATCH 
          !
          l1 = jj + i%NPATCH*(jl-1)
          k1 = jj + i%NPATCH*(jk-1)
          !
          IF ( i%XPATCH(ji,jj)>0.0 .AND. xf(ji,jj,jl+1,jk).NE.xundef .AND. xf(ji,jj,1,jk).NE.xundef ) THEN
            !
            ! Jacobian of fwd model
            zltm(l1,k1) = ( xf(ji,jj,jl+1,jk) - xf(ji,jj,1,jk) ) / zeps(ji,jj,jl)
            ! impose upper/lower limits 
            zltm(l1,k1) = max(-0.1, zltm(l1,k1))
            zltm(l1,k1) = min( 1.0, zltm(l1,k1))
            !
          ENDIF
          !
          IF (nprintlev>0) WRITE (iunit,*) zltm(l1,k1)
          !
        ENDDO
      ENDDO
    ENDDO
    !
!//////////////////////TO WRITE LTM/////////////////////////////////////   
    IF (nprintlev>0) THEN
      iunit = 120
      DO jl=1,nvar
        DO jk=1,nvar
          iunit = iunit + 1
          CLOSE(iunit)
        ENDDO
      ENDDO
      !//////////////////////TO WRITE LTM/////////////////////////////////////
      WRITE(*,*) 'LTM d(wg2)/d(wg2)', zltm(1,1)
    ENDIF
    !
    ! evolve B 
    zb(ji,:,:) = matmul(zltm(:,:),matmul(zb(ji,:,:),transpose(zltm(:,:))))
    !
    !
    !   Adding model error to background error matrix 
    zq(:,:) = 0.0
    DO jl=1,nvar
      DO jj=1,i%NPATCH
        !
        l1 = jj+i%NPATCH*(jl-1)
        !
        zq(l1,l1) = xsigma(jl)*xsigma(jl)
        !
        IF (trim(cvar(jl)) == 'LAI') THEN
          zq(l1,l1) = xscale_qlai*xscale_qlai * zq(l1,l1)
        ELSE
          zq(l1,l1) = xscale_q*xscale_q * zq(l1,l1) * zcoef(ji,jj,jl)
        ENDIF
        !
      ENDDO
    ENDDO
    !
    ! B is the forecast matrix - need to add Q
    IF ( nprintlev > 0 ) THEN
      WRITE(*,*) 'B before wg2 wg2 ==> ',zb(ji,1,1)/zcofswi(ji),zb(ji,1,1)
      WRITE(*,*) 'Q value wg2 wg2 ==> ',zq(1,1)/zcofswi(ji),zq(1,1)
    ENDIF
    !
    zb(ji,:,:) = zb(ji,:,:) + zq(:,:)
    !
    IF ( nprintlev > 0 ) WRITE(*,*) 'B after wg2 wg2 ==>',zb(ji,1,1)/zcofswi(1),zb(ji,1,1)
    !
  ENDDO
  !
  ! write out the LTM for the forward model
  ! Write out current B
  IF (nprintlev>0) THEN
    ybgfile="BGROUNDout_LBEV."//trim(ymyproc)
    CALL b_big_loop(i, &
                    "WRIT",ybgfile,zb)
    WRITE(*,*) 'store B matrix after TL evolution ==>',zb(1,1,1)
    WRITE(*,*) 'writing out B'
  ENDIF
  !
ENDIF
!
! ====================================================================
!
! 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.'//trim(ymyproc),status='unknown',iostat=istat)
DO ji = 1,ki
  IF ( minval(i%XWGI(ji,1,:))>0. ) THEN
    xyo(ji,:) = xundef
    IF ( nprintlev > 0 ) WRITE(*,*) 'OBSERVATION FOR POINT ',ji,' REMOVED'
  ENDIF
  IF ( nprintlev > 0 ) WRITE (111,*) xyo(ji,:)
ENDDO
IF ( nprintlev > 0 ) CLOSE(111)
!//////////////////////TO WRITE OBS/////////////////////////////////////
!
!############################# 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.'//trim(ymyproc),status='unknown',iostat=istat)
  ! *** Write innovations in ASCII file ***
  OPEN (unit=112,file='INNOV.'//trim(ymyproc),status='unknown',iostat=istat)
  ! Write analysis results and increments in ASCII file
  OPEN (unit=113,file='ANAL_INCR.'//trim(ymyproc),status='unknown',iostat=istat)
  ! **** Write out the observation operator + Gain matrix ****
  iunit = 150
  DO jl = 1,nvar
    DO jk=1,nobstype
      iunit = iunit + 1
      WRITE(ychar,'(I1)') jk
      yfname='HO_'//trim(cvar(jl))//'_v'//ychar
      OPEN(unit=iunit,file=yfname,form='FORMATTED',status='UNKNOWN',iostat=istat)
    ENDDO
  ENDDO
ENDIF
!//////////////////////TO WRITE ANALYSIS ARRAYS/////////////////////////////////////
!
IF (i%NPATCH==12) THEN
  zvlaimin = (/0.3,0.3,0.3,0.3,1.0,1.0,0.3,0.3,0.3,0.3,0.3,0.3/)
ELSE
  zvlaimin = (/0.3/)
ENDIF
!
ALLOCATE(i%XINCR(ki,i%NPATCH*nvar))
i%XINCR(:,:) = 0.
!
iobscount = 0
DO ji=1,ki
  !
!---------------- MEAN SIMULATED OBS AVERAGED OVER TILES-----------------------
  zyf(:,:) = 0. 
  DO jj=1,i%NPATCH
    IF (i%XPATCH(ji,jj) > 0.0) THEN
      WHERE ( xf_patch(ji,jj,:,:)/=xundef ) 
        zyf(:,:) = zyf(:,:) + i%XPATCH(ji,jj)*xf_patch(ji,jj,:,:)
      ENDWHERE
    ENDIF
  ENDDO
  IF ( nprintlev > 0 ) WRITE(*,*) 'read in sim obs yf', zyf(:,1)
  !
  !
  zr(:,:) = 0. ! Observation error matrix
  !
  zho(:,:) = xundef  ! Linearized observation matrix
  zhowr(:,:) = xundef
  zb2(:)   = xundef  ! Innovation vector
  
  DO istep=1,nboutput
    !
    DO jk = 1,nobstype
      !
      k1 = (istep-1)*nobstype + jk
      !
!--------------------- SET OBSERVATION ERROR ------------------      
      zr(k1,k1) = xerrobs(jk)*xerrobs(jk)
      IF ( cobs(jk) .EQ. "LAI" ) THEN
        zr(k1,k1) = zr(k1,k1) * xyo(ji,k1)*xyo(ji,k1)
      ELSEIF (cobs(jk) .EQ. "WG1") THEN
        ! convert R for wg1 from SWI  to abs value
        zr(k1,k1) = zr(k1,k1) * zcofswi(ji)*zcofswi(ji)
      ENDIF
      !
      ! Apply quality control
      IF ( ( abs( xyo(ji,k1)-zyf(1,jk) ) > xqcobs(jk) ) .OR. (zr(k1,k1) .LT. 0 ) ) THEN 
        xyo(ji,k1) = 999.0
      ENDIF
      !      
!--------------------- CALCULATE JACOBIANS ------------------         
      DO jl=1,nvar
        DO jj=1,i%NPATCH
          !
          l1 = jj + i%NPATCH*(jl-1)
          !
          IF ( i%XPATCH(ji,jj)>0.0 .AND. xf_patch(ji,jj,jl+1,jk).NE.xundef .AND. xf_patch(ji,jj,1,jk).NE.xundef ) THEN 
            zhowr(k1,l1) = i%XPATCH(ji,jj)*(xf_patch(ji,jj,jl+1,jk) - xf_patch(ji,jj,1,jk))/zeps(ji,jj,jl)
          ENDIF
          !
          IF( (xyo(ji,k1).NE.xundef) .AND. (xyo(ji,k1).NE.999.0) ) THEN         !if obs available
            ! Jacobian of obs operator
            zho(k1,l1) = i%XPATCH(ji,jj)*(xf_patch(ji,jj,jl+1,jk) - xf_patch(ji,jj,1,jk))/zeps(ji,jj,jl)
            ! impose limits  
            zho(k1,l1) = max(-0.1, zho(k1,l1))
            zho(k1,l1) = min( 1.0, zho(k1,l1))
            ! innovation vector
            zb2(k1) = xyo(ji,k1) - zyf(1,jk)
            iobscount = iobscount + 1
          ELSE  !if no obs available
            ! set obs operator and innovation to zero if no obs available
            zho(k1,l1) = 0.0
            zb2(k1) = 0.0 
          ENDIF
          !
        ENDDO
      ENDDO
      !
    ENDDO
    !
  ENDDO
  !
  IF ( nprintlev > 0 ) THEN
    WRITE(111,*) zr(:,:)
    WRITE(112,*) zb2(:)
  ENDIF
  
!---------------******  SOIL ANALYSIS *******--------------------------
  zhot(:,:) = 0.
  zk1(:,:) = 0.
  zp(:) = 0.
  zx(:) = 0.
  !
  zhot(:,:) = transpose(zho(:,:))
  zk1(:,:) = matmul(zho(:,:),matmul(zb(ji,:,:),zhot(:,:))) + zr(:,:)
  CALL choldc(nobstype,zk1(:,:),zp(:))                         ! Cholesky decomposition (1)
  CALL cholsl(nobstype,zk1(:,:),zp(:),zb2(:),zx(:))            ! Cholesky decomposition (2)
  i%XINCR(ji,:) = matmul(zb(ji,:,:),matmul(zhot(:,:),zx(:)))
  DO jl=1,nvar
    DO jj=1,i%NPATCH
      !
      l1 = jj+i%NPATCH*(jl-1)
      !
      ! Update the modified values
      IF ( trim(cvar(jl))=="LAI" ) THEN
        i%XINCR(ji,l1) = max( i%XINCR(ji,l1), zvlaimin(jj)-xf(ji,jj,1,jl) )
        xbio_pass(ji,jj) = xbio_pass(ji,jj) + i%XINCR(ji,l1)*xalph(jj)
      ELSEIF ( xf(ji,jj,1,jl)+i%XINCR(ji,l1)<0. ) THEN
        i%XINCR(ji,l1) = 0.
      ENDIF
      !
      xf(ji,jj,1,jl) = xf(ji,jj,1,jl) + i%XINCR(ji,l1)
      !
      ! For no only warn if we have negative values.
      IF ( nprintlev > 0 ) THEN
        IF ( xf(ji,jj,1,jl) < 0. ) WRITE(*,*) "WARNING X<0. for ",ji,jj," for variable ",trim(cvar(jl))
      ENDIF
      !
    ENDDO
  ENDDO
  !
  IF ( nprintlev > 0 ) THEN
    DO jj=1,i%NPATCH
      WRITE(113,*) (xf(ji,jj,1,jl),jl=1,nvar), (i%XINCR(ji,jj+i%NPATCH*(jl-1)),jl=1,nvar)
    ENDDO
  ENDIF
  
!--------------------ANALYSIS OF B (FOR USE IN NEXT CYCLE)-------------------
  ! Ba = (I-KH)Bf(I-KH)t+KRKt
  ! K = BfHt{K1}**-1
  ! K1 needs PATCH dim added
  
  zgain(:,:) = 0.
  zidkh(:,:) = 0.
  zkrk(:,:) = 0.
  !
  ! K1 = (R+H.B.HT) (calculate inverse -> output goes to K1)
  CALL inverse_matrix(nobs,zk1(:,:),zp(:))
  zgain(:,:) = matmul(zb(ji,:,:),matmul(zhot(:,:),zk1(:,:)))
  zidkh(:,:) = zident(:,:) - matmul(zgain(:,:),zho(:,:))
  zkrk(:,:) = matmul(zgain(:,:),matmul(zr(:,:),transpose(zgain(:,:))))
  IF (.NOT.lbfixed)  zb(ji,:,:) = matmul(zidkh(:,:),matmul(zb(ji,:,:),transpose(zidkh(:,:)))) + zkrk(:,:)
  
  IF ( nprintlev > 0 ) THEN
    iunit = 150
    DO jl = 1,nvar
      DO jk = 1,nobstype
        iunit = iunit + 1
        DO jj=1,i%NPATCH
          WRITE(iunit,*) zhowr(jk,jj+i%NPATCH*(jl-1)),zgain(jj+i%NPATCH*(jl-1),jk)
        ENDDO
      ENDDO
    ENDDO
  ENDIF
  !  
ENDDO
!
!
!//////////////////////TO WRITE ANALYSIS ARRAYS/////////////////////////////////////
IF ( nprintlev > 0 ) THEN
  CLOSE(111)
  CLOSE(112)
  CLOSE(113)
  iunit = 150
  DO jl = 1,nvar
    DO jk = 1,nobstype
      iunit = iunit + 1
      CLOSE(iunit)
    ENDDO
  ENDDO
ENDIF
!//////////////////////TO WRITE ANALYSIS ARRAYS/////////////////////////////////////
!
IF (lbev .OR. nprintlev>0) THEN
  ! Write out analysed B (for use in next cycle)
  ybgfile = "BGROUNDout_ASSIM."//trim(ymyproc)
  CALL b_big_loop(i, &
                    "WRIT",ybgfile,zb)
ENDIF
!
IF ( nprintlev > 0 ) THEN
  iobscount = iobscount / i%NPATCH / nvar
  IF (nrank==npio) THEN
    WRITE(*,*)
    WRITE(*,*) '   ---------------------------------------'
    WRITE(*,*) '   |   EXITING VARASSIM AFTER ANALYSIS   |'
    WRITE(*,*) '   ---------------------------------------'
    WRITE(*,*)
  ENDIF
  WRITE(*,*) 'Number of assimilated observations =',iobscount
  WRITE(*,*)
ENDIF
!
!############################# GET VARIABLES FOR OUTPUT WRITING ###############################
DO jl=1,nvar
  !
  ! Update the modified values
  SELECT CASE (trim(cvar(jl)))
    CASE("TG1")
      i%XTG(:,1,:) = xf(:,:,1,jl)
    CASE("TG2")
      i%XTG(:,2,:) = xf(:,:,1,jl)
    CASE("WG1")
      i%XWG(:,1,:) = xf(:,:,1,jl)
    CASE("WG2")
      i%XWG(:,2,:) = xf(:,:,1,jl)
    CASE("LAI") 
      i%XLAI(:,:) = xf(:,:,1,jl)
      SELECT CASE (trim(cbio))
        CASE("BIOMA1","BIOMASS1")
          i%XBIOMASS(:,1,:) = xbio_pass(:,:)
        CASE("BIOMA2","BIOMASS2")
          i%XBIOMASS(:,2,:) = xbio_pass(:,:)
        CASE("RESPI1","RESP_BIOM1")
          i%XRESP_BIOMASS(:,1,:) = xbio_pass(:,:)
        CASE("RESPI2","RESP_BIOM2")
          i%XRESP_BIOMASS(:,2,:) = xbio_pass(:,:)
        CASE("LAI")
          i%XLAI(:,:) = xbio_pass(:,:)
        CASE default
          CALL abor1_sfx("Mapping of "//cbio//" is not defined in EKF!")
      END SELECT
    CASE default
      CALL abor1_sfx("Mapping of "//trim(cvar(jl))//" is not defined in EKF!")
  END SELECT
ENDDO
!
#endif
!
IF (lhook) CALL dr_hook('ASSIM_NATURE_ISBA_EKF',1,zhook_handle)
!
END SUBROUTINE assim_nature_isba_ekf