init_obscovar_pdaf.F90 Source File

!-------------------------------------------------------------------------------------------
!Copyright (c) 2013-2016 by Wolfgang Kurtz and Guowei He (Forschungszentrum Juelich GmbH)
!
!This file is part of TSMP-PDAF
!
!TSMP-PDAF is free software: you can redistribute it and/or modify
!it under the terms of the GNU Lesser General Public License as published by
!the Free Software Foundation, either version 3 of the License, or
!(at your option) any later version.
!
!TSMP-PDAF is distributed in the hope that it will be useful,
!but WITHOUT ANY WARRANTY; without even the implied warranty of
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!GNU LesserGeneral Public License for more details.
!
!You should have received a copy of the GNU Lesser General Public License
!along with TSMP-PDAF.  If not, see <http://www.gnu.org/licenses/>.
!-------------------------------------------------------------------------------------------
!
!
!-------------------------------------------------------------------------------------------
!init_obscovar_pdaf.F90: TSMP-PDAF implementation of routine
!                        'init_obscovar_pdaf' (PDAF online coupling)
!-------------------------------------------------------------------------------------------

!$Id: init_obscovar_pdaf.F90 1441 2013-10-04 10:33:42Z lnerger $
!BOP
!
! !ROUTINE: init_obscovar_pdaf --- Initialize observation error covariance matrix
!
! !INTERFACE:
SUBROUTINE init_obscovar_pdaf(step, dim_obs, dim_obs_p, covar, m_state_p, &
    isdiag)

    ! !DESCRIPTION:
    ! User-supplied routine for PDAF.
    ! Used in the filters: EnKF
    !
    ! The routine is called during the analysis
    ! step when an ensemble of observations is
    ! generated by PDAF\_enkf\_obs\_ensemble.
    ! It has to initialize the global observation
    ! error covariance matrix.
    !
    ! !REVISION HISTORY:
    ! 2013-02 - Lars Nerger - Initial code
    ! Later revisions - see svn log
    !
    ! !USES:
    USE mod_assimilation, &
        ONLY: rms_obs, obs_pdaf2nc
    USE mod_parallel_pdaf, ONLY: mype_world
    USE mod_parallel_pdaf, ONLY: abort_parallel
    use mod_read_obs, only: multierr,clm_obserr, pressure_obserr
    USE mod_tsmp, ONLY: point_obs
    use netcdf

    IMPLICIT NONE

    ! !ARGUMENTS:
    INTEGER, INTENT(in) :: step                 ! Current time step
    INTEGER, INTENT(in) :: dim_obs              ! Dimension of observation vector
    INTEGER, INTENT(in) :: dim_obs_p            ! PE-local dimension of obs. vector
    REAL, INTENT(out) :: covar(dim_obs,dim_obs) ! Observation error covar. matrix
    REAL, INTENT(in) :: m_state_p(dim_obs_p)    ! Observation vector
    LOGICAL, INTENT(out) :: isdiag              ! Whether matrix R is diagonal

    ! !CALLING SEQUENCE:
    ! Called by: PDAF_enkf_obs_ensemble    (as U_init_obs_covar)
    !EOP

    ! *** local variables ***
    !   INTEGER :: i          ! Index of observation component
    !   REAL :: variance_obs  ! ariance of observations


    ! **********************
    ! *** INITIALIZATION ***
    ! **********************


    ! ******************************************************
    ! *** Initialize observation error covariance matrix ***
    ! ******************************************************

! *** local variables ***
  INTEGER :: i          ! Index of observation component
  REAL :: variance_obs  ! ariance of observations
  integer :: ncid,j,status,varid


! **********************
! *** INITIALIZATION ***
! **********************

  variance_obs = rms_obs ** 2

  covar(:, :) = 0.0


! *************************************
! ***   Initialize covariances      ***
! *************************************

! *** We assume uncorrelated Measurements here, thus R is diagonal

  !DO i = 1, dim_obs
  !   covar(i, i) = variance_obs
  !ENDDO

  if(multierr.ne.1) then
    DO i = 1, dim_obs
       covar(i, i) = variance_obs
    ENDDO
  endif

 
  if(multierr.eq.1) then

    ! Check that point observations are used
    if (.not. point_obs .eq. 1) then
      print *, "TSMP-PDAF mype(w)=", mype_world, ": ERROR(1) `point_obs.eq.1` needed for using obs_pdaf2nc."
      call abort_parallel()
    end if

    do i=1,dim_obs
#if defined CLMSA
      covar(i,i) = clm_obserr(obs_pdaf2nc(i))*clm_obserr(obs_pdaf2nc(i))
#else
      covar(i,i) = pressure_obserr(obs_pdaf2nc(i))*pressure_obserr(obs_pdaf2nc(i))
#endif
    enddo
  endif
  ! The matrix is diagonal
  ! This setting avoids the computation of the SVD of COVAR
  ! in PDAF_enkf_obs_ensemble
  isdiag = .TRUE.

!    WRITE (*,*) 'TEMPLATE init_obscovar_pdaf.F90: Set observation covariance matrix here!'
!
!    ! covar = ?
!    do i = 1, dim_obs
!        do j = 1, i
!            covar(i,j) = dot_product(wff(i,:), wff(j,:))
!            covar(j,i) = covar(i,j)
!        end do
!    end do
!    covar = covar / float(dim_obs)
!    ! Define whether the matrix is diagonal (.true./.false.)
!    isdiag = .false.

    ! kuw start
    !call check( nf90_open('corr.ncdf', nf90_nowrite, ncid) )
    !call check( nf90_inq_varid(ncid, 'corr', varid) )
    !call check( nf90_get_var(ncid,varid,covar))
    !call check(nf90_close(ncid))
    !do i=1,dim_obs
    !  do j=1,dim_obs
    !    covar(i,j) = covar(i,j)*variance_obs
    !  end do
    !end do
    !isdiag=.false.
    ! kuw end

END SUBROUTINE init_obscovar_pdaf

! subroutine check(status)

!   use netcdf
!   integer, intent ( in) :: status

!   if(status /= nf90_noerr) then
!      print *, trim(nf90_strerror(status))
!      stop "Stopped"
!   end if
! end subroutine check