PDAF_gen_obs.F90 Source File


Source Code

! Copyright (c) 2004-2024 Lars Nerger
!
! This file is part of PDAF.
!
! 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.
!
! 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 Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with PDAF.  If not, see <http://www.gnu.org/licenses/>.
!
!$Id$
!BOP
!
! !ROUTINE: PDAF_gen_obs --- Generate synthetic observations
!
! !INTERFACE:
SUBROUTINE PDAF_gen_obs(step, dim_p, dim_obs_f, dim_ens, &
     state_p, Ainv, ens_p, &
     U_init_dim_obs_f, U_obs_op_f, U_get_obs_f, U_init_obserr_f, &
     U_prepoststep, screen, flag)
  
! !DESCRIPTION:
! This routine generates observations from a model state.
! Used are the general functionality of PDAF provided with
! the call-back routines init_dim_obs, obs_op.
!
! !  This is a core routine of PDAF and
!    should not be changed by the user   !
!
! !REVISION HISTORY:
! 2019-01 - Lars Nerger - Initial code
! Later revisions - see svn log
!
! !USES:
! Include definitions for real type of different precision
! (Defines BLAS/LAPACK routines and MPI_REALTYPE)
#include "typedefs.h"

  USE PDAF_timer, &
       ONLY: PDAF_timeit, PDAF_time_temp
  USE PDAF_memcounting, &
       ONLY: PDAF_memcount
  USE PDAF_mod_filter, &
       ONLY: obs_member
  USE PDAF_mod_filtermpi, &
       ONLY: mype, dim_ens_l

  IMPLICIT NONE

! !ARGUMENTS:
! ! Variable naming scheme:
! !   suffix _p: Denotes a full variable on the PE-local domain
! !   suffix _l: Denotes a local variable on the current analysis domain
! !   suffix _f: Denotes a full variable of all observations required for the
! !              analysis loop on the PE-local domain
  INTEGER, INTENT(in) :: step        ! Current time step
  INTEGER, INTENT(in) :: dim_p       ! PE-local dimension of model state
  INTEGER, INTENT(out) :: dim_obs_f  ! PE-local dimension of observation vector
  INTEGER, INTENT(in) :: dim_ens     ! Size of ensemble
  REAL, INTENT(inout) :: state_p(dim_p)        ! PE-local model state
  REAL, INTENT(inout) :: Ainv(dim_ens, dim_ens)      ! Inverse of matrix U
  REAL, INTENT(inout) :: ens_p(dim_p, dim_ens) ! PE-local ensemble matrix
  INTEGER, INTENT(in) :: screen      ! Verbosity flag
  INTEGER, INTENT(inout) :: flag     ! Status flag

! ! External subroutines 
! ! (PDAF-internal names, real names are defined in the call to PDAF)
  EXTERNAL :: U_obs_op_f, &    ! Observation operator
       U_init_dim_obs_f, &     ! Initialize dimension of observation vector
       U_get_obs_f, &          ! Provide observation vector to user
       U_init_obserr_f, &      ! Initialize vector of observation error standard deviations
       U_prepoststep           ! User supplied pre/poststep routine

! !CALLING SEQUENCE:
! Called by:               PDAF_put_state_letkf
! Calls: U_prepoststep
! Calls: U_init_dim_obs
! Calls: U_obs_op
! Calls: U_init_obs
! Calls: PDAF_generate_rndmat
! Calls: PDAF_timeit
! Calls: PDAF_memcount
! Calls: MPI_reduce
!EOP

! *** local variables ***
  INTEGER :: i, member, row        ! Counters
  INTEGER, SAVE :: first = 0       ! Flag whether routine is called first time
  REAL    :: invdimens             ! Inverse global ensemble size
  INTEGER :: minusStep             ! Time step counter
  REAL, ALLOCATABLE :: Hx_f(:)     ! HX for PE-local model state
  REAL, ALLOCATABLE :: rndvec(:)   ! random vector to perturb observations
  REAL, ALLOCATABLE :: rms_obs(:)  ! vector of observation error standard deviations
  INTEGER, SAVE :: iseed(4)        ! seed array for random number routine
  INTEGER :: seedset = 1           ! Choice of seed set for random numbers


! ******************************
! *** Set random number seed ***
! ******************************

  IF (first == 0) THEN
     IF (seedset == 2) THEN
        iseed(1)=1
        iseed(2)=5
        iseed(3)=7
        iseed(4)=9
     ELSE IF (seedset == 3) THEN
        iseed(1)=2
        iseed(2)=5
        iseed(3)=7
        iseed(4)=9
     ELSE IF (seedset == 4) THEN
        iseed(1)=1
        iseed(2)=6
        iseed(3)=7
        iseed(4)=9
     ELSE IF (seedset == 5) THEN
        iseed(1)=1
        iseed(2)=5
        iseed(3)=8
        iseed(4)=9
     ELSE IF (seedset == 6) THEN
        iseed(1)=2
        iseed(2)=5
        iseed(3)=8
        iseed(4)=9
     ELSE IF (seedset == 7) THEN
        iseed(1)=2
        iseed(2)=6
        iseed(3)=8
        iseed(4)=9
     ELSE IF (seedset == 8) THEN
        iseed(1)=2
        iseed(2)=6
        iseed(3)=8
        iseed(4)=11
     ELSE IF (seedset == 9) THEN
        iseed(1)=3
        iseed(2)=6
        iseed(3)=8
        iseed(4)=11
     ELSE IF (seedset == 10) THEN
        iseed(1)=3
        iseed(2)=7
        iseed(3)=8
        iseed(4)=11
     ELSE IF (seedset == 11) THEN
        iseed(1)=13
        iseed(2)=7
        iseed(3)=8
        iseed(4)=11
     ELSE IF (seedset == 12) THEN
        iseed(1)=13
        iseed(2)=11
        iseed(3)=8
        iseed(4)=11
     ELSE IF (seedset == 13) THEN
        iseed(1)=13
        iseed(2)=13
        iseed(3)=8
        iseed(4)=11
     ELSE IF (seedset == 14) THEN
        iseed(1)=13
        iseed(2)=13
        iseed(3)=17
        iseed(4)=11
     ELSE IF (seedset == 15) THEN
        iseed(1)=13
        iseed(2)=13
        iseed(3)=19
        iseed(4)=11
     ELSE IF (seedset == 16) THEN
        iseed(1)=15
        iseed(2)=13
        iseed(3)=19
        iseed(4)=11
     ELSE IF (seedset == 17) THEN
        iseed(1)=15
        iseed(2)=135
        iseed(3)=19
        iseed(4)=11
     ELSE IF (seedset == 18) THEN
        iseed(1)=19
        iseed(2)=135
        iseed(3)=19
        iseed(4)=11
     ELSE IF (seedset == 19) THEN
        iseed(1)=19
        iseed(2)=135
        iseed(3)=19
        iseed(4)=17
     ELSE IF (seedset == 20) THEN
        iseed(1)=15
        iseed(2)=15
        iseed(3)=47
        iseed(4)=17
     ELSE
        ! Standard seed
        iseed(1) = 1000
        iseed(2) = 2034
        iseed(3) = 0
        iseed(4) = 3
     END IF
     first = 2
  END IF

! *************************************
! *** Prestep for forecast ensemble ***
! *************************************


  CALL PDAF_timeit(5, 'new')
  minusStep = - step  ! Indicate forecast by negative time step number
  IF (mype == 0 .AND. screen > 0) THEN
     WRITE (*, '(a, 5x, a, i7)') 'PDAF', 'Call pre-post routine after forecast; step ', step
  ENDIF
  CALL U_prepoststep(minusStep, dim_p, dim_ens, dim_ens_l, dim_obs_f, &
       state_p, Ainv, ens_p, flag)
  CALL PDAF_timeit(5, 'old')

  IF (mype == 0 .AND. screen > 0) THEN
     IF (screen > 1) THEN
        WRITE (*, '(a, 5x, a, F10.3, 1x, a)') &
             'PDAF ', '--- duration of prestep:', PDAF_time_temp(5), 's'
     END IF
     WRITE (*, '(a, 55a)') 'PDAF generate observations ', ('-', i = 1, 42)
  END IF


! ********************************************************************
! *** Generate synthetic observations                              ***
! *** Steps:                                                       ***
! ***   1. Get observed state by applying the observation operator ***
! ***   2. Generate Gaussian random number                         ***
! ***   3. Get vector of observation error standard deviations     ***
! ***   4. Add random noise scaled by obs. error to observed state ***
! ***   5. Provide observation vector to user                      ***
! ********************************************************************

  CALL PDAF_timeit(3, 'new')

! *** 1. Get full observed state vector ***

  ! Get full observation dimension
  CALL U_init_dim_obs_f(step, dim_obs_f)

  IF (screen > 0) THEN
     IF (screen<=2 .AND. mype == 0) THEN
        WRITE (*, '(a, 5x, a, i6, a, i10)') &
             'PDAF', '--- PE-Domain:', mype, &
             ' dimension of PE-local full obs. vector', dim_obs_f
     ELSE IF (screen>2) THEN
        WRITE (*, '(a, 5x, a, i6, a, i10)') &
             'PDAF', '--- PE-Domain:', mype, &
             ' dimension of PE-local full obs. vector', dim_obs_f
     END IF
  END IF

  CALL PDAF_timeit(11, 'new')

  ! Compute mean forecast state
  state_p = 0.0
  invdimens = 1.0 / REAL(dim_ens)
  DO member = 1, dim_ens
     DO row = 1, dim_p
        state_p(row) = state_p(row) + invdimens * ens_p(row, member)
     END DO
  END DO

  CALL PDAF_timeit(11, 'old')

  CALL PDAF_timeit(33, 'new')

  ! Apply observation operator to get Hx for full DIM_OBS_F region on PE-local domain
  ALLOCATE(Hx_f(dim_obs_f))
  IF (first == 0) CALL PDAF_memcount(3, 'r', dim_obs_f)

  obs_member = 1  ! Set this for completeness

  CALL U_obs_op_f(step, dim_p, dim_obs_f, state_p, Hx_f)


! *** 2. Initialize random vector ***

  ALLOCATE(rndvec(dim_obs_f))
  IF (first == 0) CALL PDAF_memcount(3, 'r', dim_ens)

  CALL larnvTYPE(3, iseed, dim_obs_f, rndvec)

! *** 3. Get vector of observation errors ***

  ALLOCATE(rms_obs(dim_obs_f))
  IF (first == 0) CALL PDAF_memcount(3, 'r', dim_ens)

  CALL U_init_obserr_f(step, dim_obs_f, Hx_f, rms_obs)

! *** 4. Generate perturbed observations ***
  DO i = 1, dim_obs_f
     Hx_f(i) = Hx_f(i) + rms_obs(i)*rndvec(i)
  END DO

! *** Provide observation vector to user ***
  CALL U_get_obs_f(step, dim_obs_f, Hx_f)

  CALL PDAF_timeit(33, 'old')

  CALL PDAF_timeit(3, 'old')

  IF (mype == 0 .AND. screen > 0) THEN
     WRITE (*, '(a, 55a)') 'PDAF Forecast ', ('-', i = 1, 55)
  END IF


! ********************
! *** Finishing up ***
! ********************

  DEALLOCATE(Hx_f, rndvec, rms_obs)

  IF (first == 0) first = 1

END SUBROUTINE PDAF_gen_obs