assimilate_pdaf.F90 Source File


Source Code

!-------------------------------------------------------------------------------------------
!Copyright (c) 2013-2016 by Wolfgang Kurtz, Guowei He and Mukund Pondkule (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/>.
!-------------------------------------------------------------------------------------------
!
!
!-------------------------------------------------------------------------------------------
!assimilate_pdaf.F90: TSMP-PDAF implementation of routine
!                     'assimilate_pdaf' (PDAF online coupling)
!-------------------------------------------------------------------------------------------

!$Id$
!BOP
!
! !ROUTINE: assimilate_pdaf - Routine to control perform analysis step
!
! !INTERFACE:
SUBROUTINE assimilate_pdaf()

! !DESCRIPTION:
! This routine is called during the model integrations at each time 
! step. It check whether the forecast phase is completed. If so, 
! PDAF_put_state_X is called to perform the analysis step.
!
! !REVISION HISTORY:
! 2013-08 - Lars Nerger - Initial code
! Later revisions - see svn log
!
! !USES:
  USE mod_parallel_pdaf, &     ! Parallelization variables
       ONLY: abort_parallel, mype_world
  USE mod_assimilation, &      ! Variables for assimilation
       ONLY: filtertype
  ! USE PDAF_interfaces_module   ! Check consistency of PDAF calls

  IMPLICIT NONE

! !CALLING SEQUENCE:
! Called by: step
! Calls: PDAF_assimilate_X
! Calls: PDAF_set_debug_flag
!EOP

! Local variables
  INTEGER :: status_pdaf       ! PDAF status flag


! ! External subroutines
! !  (subroutine names are passed over to PDAF in the calls to 
! !  PDAF_get_state and PDAF_assimilate_X. This allows the user 
! !  to specify the actual name of a routine. However, the 
! !  PDAF-internal name of a subroutine might be different from
! !  the external name!)
!
! ! Subroutines used with all filters
  EXTERNAL :: collect_state_pdaf, &    ! Routine to collect a state vector from model fields
       init_dim_obs_pdaf, &            ! Initialize dimension of observation vector
       obs_op_pdaf, &                  ! Implementation of the Observation operator
       init_obs_pdaf, &                ! Routine to provide vector of measurements
       prepoststep_ens_pdaf, &         ! User supplied pre/poststep routine
       prodRinvA_pdaf, &               ! Provide product R^-1 A for some matrix A
       init_obsvar_pdaf, &             ! Initialize mean observation error variance
       next_observation_pdaf, &        ! Provide time step, model time, and dimension of next observation
       distribute_state_pdaf           ! Routine to distribute a state vector to model fields
! ! Subroutines for local filters
  EXTERNAL :: init_n_domains_pdaf, &   ! Provide number of local analysis domains
       init_dim_l_pdaf, &              ! Initialize state dimension for local ana. domain
       init_dim_obs_l_pdaf,&           ! Initialize dim. of obs. vector for local ana. domain
       g2l_state_pdaf, &               ! Get state on local ana. domain from global state
       l2g_state_pdaf, &               ! Init global state from state on local analysis domain
       g2l_obs_pdaf, &                 ! Restrict a global obs. vector to local analysis domain
       init_obs_l_pdaf, &              ! Provide vector of measurements for local ana. domain
       prodRinvA_l_pdaf, &             ! Provide product R^-1 A for some local matrix A
       init_obsvar_l_pdaf, &           ! Initialize local mean observation error variance
       init_obs_f_pdaf, &              ! Provide full vector of measurements for PE-local domain
       obs_op_f_pdaf, &                ! Obs. operator for full obs. vector for PE-local domain
       init_dim_obs_f_pdaf             ! Get dimension of full obs. vector for PE-local domain
! ! Subroutines used in EnKF
  EXTERNAL :: add_obs_error_pdaf, &    ! Add obs. error covariance R to HPH in EnKF
       init_obscovar_pdaf              ! Initialize obs error covar R in EnKF
! ! Subroutines used for localization in LEnKF
  EXTERNAL :: localize_covar_pdaf       ! Apply localization to HP and HPH^T
! ! ! Subroutines used in NETF
!   EXTERNAL :: likelihood_pdaf          ! Compute observation likelihood for an ensemble member
! ! ! Subroutines used in LNETF
!   EXTERNAL :: likelihood_l_pdaf        ! Compute local observation likelihood for an ensemble member
! ! ! Subroutines used in LKNETF
!   EXTERNAL :: likelihood_hyb_l_pdaf, & ! Compute local likelihood awith hybrid weight for an ensemble member
!        prodRinvA_hyb_l_pdaf            ! Provide product R^-1 A for some matrix A including hybrid weight

! *** Switch on debug output ***
! *** for main process        ***
#ifdef PDAF_DEBUG
  IF (mype_world == 0) CALL PDAF_set_debug_flag(1)
#endif

! *********************************
! *** Call assimilation routine ***
! *********************************

  ! IF (filtertype == 1) THEN
  !    CALL PDAF_assimilate_seik(collect_state_pdaf, distribute_state_pdaf, &
  !         init_dim_obs_pdaf, obs_op_pdaf, init_obs_pdaf, prepoststep_ens_pdaf, &
  !         prodRinvA_pdaf, next_observation_pdaf, status_pdaf)
  IF (filtertype == 2) THEN
    CALL PDAF_assimilate_enkf(collect_state_pdaf, distribute_state_pdaf, &
      init_dim_obs_pdaf, obs_op_pdaf, &
      init_obs_pdaf, prepoststep_ens_pdaf, add_obs_error_pdaf, init_obscovar_pdaf, &
      next_observation_pdaf, status_pdaf)
  ! ELSE IF (filtertype == 3) THEN
  !    CALL PDAF_assimilate_lseik(collect_state_pdaf, distribute_state_pdaf, &
  !         init_dim_obs_f_pdaf, obs_op_f_pdaf, &
  !         init_obs_f_pdaf, init_obs_l_pdaf, prepoststep_ens_pdaf, &
  !         prodRinvA_l_pdaf, init_n_domains_pdaf, init_dim_l_pdaf, &
  !         init_dim_obs_l_pdaf, g2l_state_pdaf, l2g_state_pdaf, &
  !         g2l_obs_pdaf, init_obsvar_pdaf, init_obsvar_l_pdaf, next_observation_pdaf, status_pdaf)
  ELSE IF (filtertype == 4) THEN
     CALL PDAF_assimilate_etkf(collect_state_pdaf, distribute_state_pdaf, &
          init_dim_obs_pdaf, obs_op_pdaf, &
          init_obs_pdaf, prepoststep_ens_pdaf, prodRinvA_pdaf, init_obsvar_pdaf, next_observation_pdaf, status_pdaf)
  ELSE IF (filtertype == 5) THEN
     CALL PDAF_assimilate_letkf(collect_state_pdaf, distribute_state_pdaf, &
          init_dim_obs_f_pdaf, obs_op_f_pdaf, &
          init_obs_f_pdaf, init_obs_l_pdaf, prepoststep_ens_pdaf, &
          prodRinvA_l_pdaf, init_n_domains_pdaf, init_dim_l_pdaf, &
          init_dim_obs_l_pdaf, g2l_state_pdaf, l2g_state_pdaf, &
          g2l_obs_pdaf, init_obsvar_pdaf, init_obsvar_l_pdaf, next_observation_pdaf, status_pdaf)
  ELSE IF (filtertype == 6) THEN
     CALL PDAF_assimilate_estkf(collect_state_pdaf, distribute_state_pdaf, &
          init_dim_obs_pdaf, obs_op_pdaf, &
          init_obs_pdaf, prepoststep_ens_pdaf, prodRinvA_pdaf, init_obsvar_pdaf, next_observation_pdaf, status_pdaf)
  ELSE IF (filtertype == 7) THEN
     CALL PDAF_assimilate_lestkf(collect_state_pdaf, distribute_state_pdaf, &
          init_dim_obs_f_pdaf, obs_op_f_pdaf, &
          init_obs_f_pdaf, init_obs_l_pdaf, prepoststep_ens_pdaf, &
          prodRinvA_l_pdaf, init_n_domains_pdaf, init_dim_l_pdaf, &
          init_dim_obs_l_pdaf, g2l_state_pdaf, l2g_state_pdaf, &
          g2l_obs_pdaf, init_obsvar_pdaf, init_obsvar_l_pdaf, next_observation_pdaf, status_pdaf)
  ELSE IF (filtertype == 8) THEN
     CALL PDAF_assimilate_lenkf(collect_state_pdaf, distribute_state_pdaf, &
          init_dim_obs_pdaf, obs_op_pdaf, &
          init_obs_pdaf, prepoststep_ens_pdaf, localize_covar_pdaf, add_obs_error_pdaf, &
          init_obscovar_pdaf, next_observation_pdaf, status_pdaf)
  ! ELSE IF (filtertype == 9) THEN
  !    CALL PDAF_assimilate_netf(collect_state_pdaf, distribute_state_pdaf, &
  !         init_dim_obs_pdaf, &
  !         obs_op_pdaf, init_obs_pdaf, prepoststep_ens_pdaf, &
  !         likelihood_pdaf, next_observation_pdaf, status_pdaf)
  ! ELSE IF (filtertype == 10) THEN
  !    CALL PDAF_assimilate_lnetf(collect_state_pdaf, distribute_state_pdaf, &
  !         init_dim_obs_f_pdaf, &
  !         obs_op_f_pdaf, init_obs_l_pdaf, prepoststep_ens_pdaf, &
  !         likelihood_l_pdaf, init_n_domains_pdaf, init_dim_l_pdaf, &
  !         init_dim_obs_l_pdaf, g2l_state_pdaf, l2g_state_pdaf, &
  !         g2l_obs_pdaf, next_observation_pdaf, status_pdaf)
  ! ELSE IF (filtertype == 11) THEN
  !    CALL PDAF_assimilate_lknetf(collect_state_pdaf, distribute_state_pdaf, &
  !         init_dim_obs_f_pdaf, obs_op_f_pdaf, &
  !         init_obs_f_pdaf, init_obs_l_pdaf, prepoststep_ens_pdaf, &
  !         prodRinvA_l_pdaf, prodRinvA_hyb_l_pdaf, init_n_domains_pdaf, &
  !         init_dim_l_pdaf, init_dim_obs_l_pdaf, g2l_state_pdaf, l2g_state_pdaf, &
  !         g2l_obs_pdaf, init_obsvar_pdaf, init_obsvar_l_pdaf, &
  !         likelihood_l_pdaf, likelihood_hyb_l_pdaf, next_observation_pdaf, status_pdaf)
  ! ELSE IF (filtertype == 12) THEN
  !    CALL PDAF_assimilate_pf(collect_state_pdaf, distribute_state_pdaf, &
  !         init_dim_obs_pdaf, &
  !         obs_op_pdaf, init_obs_pdaf, prepoststep_ens_pdaf, &
  !         likelihood_pdaf, next_observation_pdaf, status_pdaf)
  END IF

  ! Check for errors during execution of PDAF

  IF (status_pdaf /= 0) THEN
     WRITE (*,'(/1x,a6,i3,a43,i4,a1/)') &
          'ERROR ', status_pdaf, &
          ' in PDAF_assimilate - stopping! (PE ', mype_world,')'
     CALL  abort_parallel()
  END IF

! *** Switch off debug output ***
#ifdef PDAF_DEBUG
  CALL PDAF_set_debug_flag(0)
#endif

END SUBROUTINE assimilate_pdaf