init_pdaf.F90 Source File


Source Code

!-------------------------------------------------------------------------------------------
!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_pdaf.F90: TSMP-PDAF implementation of routine
!               'init_pdaf' (PDAF online coupling)
!-------------------------------------------------------------------------------------------

!$Id: init_pdaf.F90 1444 2013-10-04 10:54:08Z lnerger $
!BOP
!
! !ROUTINE: init_pdaf - Interface routine to call initialization of PDAF
!
! !INTERFACE:
SUBROUTINE init_pdaf()

! !DESCRIPTION:
! This routine collects the initialization of variables for PDAF.
! In addition, the initialization routine PDAF_init is called
! such that the internal initialization of PDAF is performed.
! This variant is for the online mode of PDAF.
!
! This routine is generic. However, it assumes a constant observation
! error (rms_obs). Further, with parallelization the local state
! dimension dim_state_p is used.
!
! !TSMP-PDAF-DESCRIPTION:
! This routine initializes a pointer to the state vector that is set
! by component-model-specific routines in `initialize_tsmp`. 
!
! This routine sets the local and global state vector dimension.
!
! Debug output for this routine is turned on by preprocessor flag
! `PDAF_DEBUG`.
!
! !REVISION HISTORY:
! 2008-10 - Lars Nerger - Initial code
! Later revisions - see svn log
!
! !USES:
!   USE mod_model, &             ! Model variables
!        ONLY: nx, ny, nx_p
  
  USE mod_parallel_pdaf, &     ! Parallelization variables for
    ! assimilation
        ONLY: n_modeltasks, task_id, COMM_filter, COMM_couple, filterpe, &
        abort_parallel, &
        mype_world, COMM_model, npes_model, &
        mpi_success, mpi_comm_world, mpi_integer, mype_model
  USE mod_assimilation, &      ! Variables for assimilation
        ONLY: dim_state_p, dim_state, screen, filtertype, subtype, toffset,&
        dim_ens, rms_obs, model_error, model_err_amp, incremental, &
        type_forget, forget, dim_bias, rank_analysis_enkf, &
        locweight, cradius, sradius, filename, &
        type_trans, type_sqrt, delt_obs, toffset, dim_state_p_count, &
        dim_lag, &
        type_winf, limit_winf, &
        type_hyb, hyb_gamma, hyb_kappa, &
        pf_res_type, pf_noise_type, pf_noise_amp
  USE mod_tsmp, &
        ONLY: pf_statevecsize, nprocpf, tag_model_parflow, tag_model_clm, nprocclm, pf_statevec, pf_statevec_fortran, &
        idx_map_subvec2state, idx_map_subvec2state_fortran, model
#if defined CLMSA
  ! kuw: get access to clm variables
#ifndef CLMFIVE    
  USE shr_kind_mod , only : r8 => shr_kind_r8
  USE clm_atmlnd   , only : clm_l2a, atm_l2a, clm_mapl2a
  USE clmtype      , only : clm3, nameg
  USE subgridAveMod, only : p2g, c2g
  USE domainMod    , only : latlon_type
  USE clm_varpar   , only : nlevsoi
  USE decompMod    , only : get_proc_global, get_proc_bounds, adecomp
  USE spmdGathScatMod , only : gather_data_to_master
  USE spmdMod      , only : masterproc
#endif
  use enkf_clm_mod, only: clm_statevecsize
#endif
  ! kuw end

  use, intrinsic :: iso_c_binding

  IMPLICIT NONE

! !CALLING SEQUENCE:
! Called by: main
! Calls: init_pdaf_parse
! Calls: init_pdaf_info
! Calls: PDAF_init
! Calls: PDAF_get_state
! Calls: PDAF_set_debug_flag
!EOP

! Local variables
  INTEGER :: filter_param_i(7) ! Integer parameter array for filter
  REAL    :: filter_param_r(2) ! Real parameter array for filter
  INTEGER :: status_pdaf       ! PDAF status flag
  INTEGER :: doexit, steps     ! Not used in this implementation
  REAL    :: timenow           ! Not used in this implementation
  integer :: ierror

  ! External subroutines
  EXTERNAL :: init_ens         ! Ensemble initialization
  EXTERNAL :: 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
    prepoststep_ens_pdaf            ! User supplied pre/poststep routine
  integer :: i, j
  !kuw: dimensions for clm variables
  integer :: numg           ! total number of gridcells across all processors
  integer :: numl           ! total number of landunits across all processors
  integer :: numc           ! total number of columns across all processors
  integer :: nump           ! total number of pfts across all processors
  integer :: begg,endg      ! local beg/end gridcells gdc
  integer :: begl,endl      ! local beg/end landunits
  integer :: begc,endc      ! local beg/end columns
  integer :: begp,endp      ! local beg/end pfts
  !kuw end
  CHARACTER (LEN = 30) :: fn    !TSMP-PDAF: function name for output of index array

! ***************************
! ***   Initialize PDAF   ***
! ***************************

  IF (mype_world == 0) THEN
     WRITE (*,'(/1x,a)') 'INITIALIZE PDAF - ONLINE MODE'
  END IF

! *** Pointer initialization for ParFlow-type state vector ***
  if (model == tag_model_parflow) then
    ! Parflow: Initialize Fortran-pointer on pf_statevec
    call C_F_POINTER(pf_statevec, pf_statevec_fortran, [pf_statevecsize])

    ! Parflow: Initialize Fortran-pointer on idx_mapping_subvec2state
    call C_F_POINTER(idx_map_subvec2state, idx_map_subvec2state_fortran, [pf_statevecsize])
  end if

#ifdef PDAF_DEBUG
  ! Index array only needed for a single modeltask
  IF(filterpe) THEN

    ! idx_map_subvec2state is only set for ParFlow processes
    IF(model == tag_model_parflow) THEN
      ! Debug output: index manipulation array subvec->state
      WRITE(fn, "(a,i5.5,a)") "idx_map_subvec2state_", mype_world, ".txt"
      OPEN(unit=71, file=fn, action="write")
      DO i = 1, pf_statevecsize
        WRITE (71,"(i10)") idx_map_subvec2state_fortran(i)
      END DO
      CLOSE(71)
    END IF

  END IF
#endif

! *** Define state dimension ***

! *** Setting local state vector dimension ***
  if (model == tag_model_parflow) then
    ! Parflow component, setting local state dimension `dim_state_p`
    ! and later dim_state from `pf_statevecsize` from `initialize_tsmp
    ! -> parflow_oasis_init`.
    dim_state_p = pf_statevecsize  ! Local state dimension
  else
    ! CLM/COSMO component, setting dummy dim_state_p and dim_state
    dim_state_p = 1  ! Local state dimension
  end if

#if defined CLMSA
  if (model == tag_model_clm) then

    ! CLM component: setting local state dimension from
    ! `clm_statevecsize` from `initialize_tsmp -> clm(5)_init ->
    ! define_clm_statevec`
    dim_state_p = clm_statevecsize

  end if
#endif

! *** Setting global state vector dimension ***
  IF (allocated(dim_state_p_count)) deallocate(dim_state_p_count)
  allocate(dim_state_p_count(npes_model))
  call MPI_Gather(dim_state_p, 1, MPI_INTEGER, dim_state_p_count, 1, MPI_INTEGER, 0, COMM_model, ierror)

#ifdef PDAF_DEBUG
  ! Debug output: local state dimension array
  if (mype_model == 0) WRITE(*, '(a,x,a,i5,x,a,x)', advance="no") "TSMP-PDAF-debug", "mype(w)=", mype_world, "init_pdaf: dim_state_p_count in modified:"
  if (mype_model == 0) WRITE(*, *) dim_state_p_count
#endif

  if (mype_model == 0) then
    dim_state = sum(dim_state_p_count)
  end if
  call MPI_BCAST(dim_state, 1, MPI_INTEGER, 0, COMM_model, IERROR)

#ifdef PDAF_DEBUG
  ! Debug output: global state dimension
  WRITE(*, '(a,x,a,i5,x,a,x,i9)') "TSMP-PDAF-debug", "mype(w)=", mype_world, "init_pdaf: my local state vector dimension dim_state_p:", dim_state_p
  WRITE(*, '(a,x,a,i5,x,a,2x,i9)') "TSMP-PDAF-debug", "mype(w)=", mype_world, "init_pdaf: my global state vector dimension dim_state:", dim_state
#endif

  call MPI_Barrier(MPI_COMM_WORLD, ierror)

! **********************************************************
! ***   CONTROL OF PDAF - used in call to PDAF_init      ***
! **********************************************************

! *** IO options ***
  screen      = 2    ! Write screen output (1) for output, (2) add timings

! *** Ensemble size ***
  dim_ens = n_modeltasks        ! Size of ensemble for all ensemble filters

! *** Options for filter method

  ! ++++++++++++++++++++++++++++++++++++++++++++++++++
  ! +++ For available options see MOD_ASSIMILATION +++
  ! ++++++++++++++++++++++++++++++++++++++++++++++++++

  filtertype = 2     ! Type of filter
  subtype = 0        ! Subtype of filter

  forget  = 1.0      ! Forgetting factor value for inflation
  type_forget = 0    ! Type of forgetting factor

  type_trans = 0     ! Type of ensemble transformation (deterministic or random)
  type_sqrt = 0      ! SEIK/LSEIK/ESTKF/LESTKF: Type of transform matrix square-root
  incremental = 0    ! SEIK/LSEIK: (1) to perform incremental updating

  !EnKF
  rank_analysis_enkf = 0  ! EnKF: rank to be considered for inversion of HPH in analysis step

  ! NETF/LNETF/PF
  type_winf = 0      ! NETF/LNETF/PF: Type of weights inflation
  limit_winf = 0.0   ! NETF/LNETF/PF: Limit for weights inflation

  ! LKNETF
  type_hyb = 0       ! LKNETF: Type of hybrid weight
  hyb_gamma =  1.0   ! LKNETF: Hybrid filter weight for state (1.0: LETKF, 0.0: LNETF)
  hyb_kappa = 30.0   ! LKNETF: Hybrid norm for using skewness and kurtosis (type_hyb 3 or 4)

  ! PF
  pf_res_type = 1    ! PF: Resampling type for particle filter
  pf_noise_type = 0  ! PF: Type of pertubing noise
  pf_noise_amp = 0.0 ! PF: Noise amplitude for particle filter


! *********************************************************************
! ***   Settings for analysis steps  - used in call-back routines   ***
! *********************************************************************

! *** Forecast length (time interval between analysis steps) ***
  delt_obs = 2     ! Number of time steps between analysis/assimilation steps

! *** specifications for observations ***
  rms_obs = 0.5    ! Observation error standard deviation

! *** Localization settings
  locweight = 0     ! Type of localizating weighting
                    !   (0) constant weight of 1
                    !   (1) exponentially decreasing with SRADIUS
                    !   (2) use 5th-order polynomial
                    !   (3) regulated localization of R with mean error variance
                    !   (4) regulated localization of R with single-point error variance
  cradius = 0.0     ! Cut-off radius for observation domain in local filters
  sradius = cradius ! Support radius for 5th-order polynomial
                    ! or radius for 1/e for exponential weighting

! *** File names
  filename = 'output.dat'

! *** TSMP-PDAF-specific inputs
  !kuw: add smoother support
  dim_lag = 0
  !kuw end

  ! hcp
  toffset = 0      ! offset of time steps shifting all analysis/assimilation steps
  ! hcp end
  
! ***********************************
! *** Some optional functionality ***
! ***********************************

! *** Parse command line options   ***
! *** This is optional, but useful ***

  call init_pdaf_parse()


! *** Initial Screen output ***
! *** This is optional      ***

  IF (mype_world == 0) call init_pdaf_info()

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

! *****************************************************
! *** Call PDAF initialization routine on all PEs.  ***
! ***                                               ***
! *** Here, the full selection of filters is        ***
! *** implemented. In a real implementation, one    ***
! *** reduce this to selected filters.              ***
! ***                                               ***
! *** For all filters, first the arrays of integer  ***
! *** and real number parameters are initialized.   ***
! *** Subsequently, PDAF_init is called.            ***
! *****************************************************

  whichinit: IF (filtertype == 2 .or. filtertype==8) THEN
    ! *** EnKF with Monte Carlo init ***
    filter_param_i(1) = dim_state_p ! State dimension
    filter_param_i(2) = dim_ens     ! Size of ensemble
    filter_param_i(3) = rank_analysis_enkf ! Rank of speudo-inverse in analysis
    filter_param_i(4) = incremental ! Whether to perform incremental analysis
    filter_param_i(5) = 0           ! Smoother lag (not implemented here)
    !kuw: add smoother support
    filter_param_i(5) = dim_lag     ! Smoother lag (not implemented here)
    !kuw end
    filter_param_r(1) = forget      ! Forgetting factor

    !hcp 0-> toffset
    CALL PDAF_init(filtertype, subtype, toffset, &
      filter_param_i, 6,&
      filter_param_r, 2, &
      COMM_model, COMM_filter, COMM_couple, &
      task_id, n_modeltasks, filterpe, init_ens, &
      screen, status_pdaf)
  ELSE
    ! *** All other filters                       ***
    ! *** SEIK, LSEIK, ETKF, LETKF, ESTKF, LESTKF ***
    filter_param_i(1) = dim_state_p ! State dimension
    filter_param_i(2) = dim_ens     ! Size of ensemble
    filter_param_i(3) = 0           ! Smoother lag (not implemented here)
    filter_param_i(4) = incremental ! Whether to perform incremental analysis
    filter_param_i(5) = type_forget ! Type of forgetting factor
    filter_param_i(6) = type_trans  ! Type of ensemble transformation
    filter_param_i(7) = type_sqrt   ! Type of transform square-root (SEIK-sub4/ESTKF)
    filter_param_r(1) = forget      ! Forgetting factor

    !hcp 0-> toffset
    CALL PDAF_init(filtertype, subtype, toffset, &
      filter_param_i, 7,&
      filter_param_r, 2, &
      COMM_model, COMM_filter, COMM_couple, &
      task_id, n_modeltasks, filterpe, init_ens, &
      screen, status_pdaf)
  END IF whichinit


! *** Check whether initialization of PDAF was successful ***
  IF (status_pdaf /= 0) THEN
     WRITE (*,'(/1x,a6,i3,a43,i4,a1/)') &
          'ERROR ', status_pdaf, &
          ' in initialization of PDAF - stopping! (PE ', mype_world,')'
     CALL abort_parallel()
  END IF

! ******************************'***
! *** Prepare ensemble forecasts ***
! ******************************'***

  CALL PDAF_get_state(steps, timenow, doexit, next_observation_pdaf, &
       distribute_state_pdaf, prepoststep_ens_pdaf, status_pdaf)

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

  if (mype_world == 0 .and. screen > 2) then
    print *, "TSMP-PDAF INITIALIZE PDAF FINISHED"
  end if

END SUBROUTINE init_pdaf