! 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_get_state --- Interface to control ensemble integration ! ! !INTERFACE: SUBROUTINE PDAF_get_state(steps, time, doexit, U_next_observation, U_distribute_state, & U_prepoststep, outflag) ! !DESCRIPTION: ! Interface routine called from the model before the ! forecast of each ensemble state to transfer data ! from PDAF to the model. For the parallelization ! this involves transfer from filter PEs to model PEs.\\ ! At the beginning of a forecast phase sub-ensembles ! are distributed to the model tasks. During the ! forecast phase each state vector of a sub-ensemble ! is transferred to the model fields by U\_dist\_state. ! ! This version is for PDAF with domain-decomposition. ! ! ! This is a core routine of PDAF and ! should not be changed by the user ! ! ! !REVISION HISTORY: ! 2003-07 - 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 mpi USE PDAF_communicate_ens, & ONLY: PDAF_scatter_ens USE PDAF_timer, & ONLY: PDAF_timeit, PDAF_time_temp USE PDAF_mod_filter, & ONLY: dim_p, dim_eof, dim_ens, local_dim_ens, dim_obs, nsteps, & step_obs, step, member_get, member_put=>member, member_save, subtype_filter, & ensemblefilter, initevol, epsilon, state, eofV, eofU, & firsttime, end_forecast, screen, flag, dim_lag, sens, & cnt_maxlag, cnt_steps, debug, offline_mode USE PDAF_mod_filtermpi, & ONLY: mype_world, mype_model, task_id, statetask, filterpe, & modelpe, dim_eof_l, dim_ens_l IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(inout) :: steps ! Flag and number of time steps REAL, INTENT(out) :: time ! current model time INTEGER, INTENT(inout) :: doexit ! Whether to exit from forecasts INTEGER, INTENT(inout) :: outflag ! Status flag ! ! External subroutines ! ! (PDAF-internal names, real names are defined in the call to PDAF) EXTERNAL :: U_next_observation, & ! Routine to provide time step, time and dimension ! of next observation U_distribute_state, & ! Routine to distribute a state vector U_prepoststep ! User supplied pre/poststep routine ! !CALLING SEQUENCE: ! Called by: model code ! Calls: PDAF_timeit ! Calls: PDAF_time_temp ! Calls: U_prepoststep ! Calls: U_next_observation ! Calls: U_distribute_state ! Calls: MPI_bcast (MPI) !EOP ! local variables INTEGER :: i, j ! Counters LOGICAL :: central_state ! Perform evolution of only central state in SEEK ! ********************** ! *** INITIALIZATION *** ! ********************** IF (debug>0) & WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_get_state -- START' firstt: IF (firsttime == 1 .AND. .NOT.offline_mode) THEN ! Screen output IF (mype_world == 0 .AND. screen > 0) THEN WRITE (*, '(//a5, 64a)') 'PDAF ',('-', i = 1, 64) WRITE (*, '(a, 20x, a)') 'PDAF', '+++++ ASSIMILATION +++++' WRITE (*, '(a5, 64a)') 'PDAF ', ('-', i = 1, 64) END IF ! *** Initial prestep IF (filterpe) THEN CALL PDAF_timeit(5, 'new') IF (mype_world == 0 .AND. screen > 0) THEN WRITE (*, '(a, 5x, a)') 'PDAF','Call pre-post routine at initial time' ENDIF typef: IF (ensemblefilter) THEN ! *** Ensemble-based filter *** CALL U_prepoststep(step_obs, dim_p, dim_ens, dim_ens_l, dim_obs, & state, eofU, eofV, flag) ELSE ! *** Mode-based filter (SEEK) *** CALL U_prepoststep(step_obs, dim_p, dim_eof, dim_eof_l, dim_obs, & state, eofU, eofV, flag) END IF typef IF (debug>0) & WRITE (*,*) '++ PDAF-debug get_state: ', debug, 'prepoststep completed' CALL PDAF_timeit(5, 'old') END IF IF (mype_world == 0 .AND. screen > 0) THEN IF (screen >= 2) THEN WRITE (*, '(a, 5x, a, F10.3, 1x, a)') & 'PDAF','--- duration of prestep:', PDAF_time_temp(5), 's' END IF WRITE (*, '(a, 55a)') 'PDAF Forecast ',('-', i = 1, 55) END IF ! *** Unset first time flag firsttime = 0 END IF firstt ! *** Initialize local ensemble sizes *** IF (ensemblefilter) THEN local_dim_ens = dim_ens_l IF (debug>0) & WRITE (*,*) '++ PDAF-debug get_state: ', debug, 'task-local ensemble size', & local_dim_ens ELSE ! For SEEK take care of central state IF (task_id /= statetask) THEN local_dim_ens = dim_eof_l IF (debug>0) & WRITE (*,*) '++ PDAF-debug get_state: ', debug, 'task-local ensemble size', & local_dim_ens ELSE local_dim_ens = dim_eof_l + 1 IF (debug>0) & WRITE (*,*) '++ PDAF-debug get_state: ', debug, 'task-local ensemble size', & local_dim_ens, 'including central state' END IF END IF IF (offline_mode) THEN ! For offline mode of PDAF skip initialization ! of forecast and loop through ensemble states initevol = 0 nsteps = 0 step_obs = 0 step = 0 END IF ! ********************************* ! *** Initialize forecast phase *** ! ********************************* evolinit: IF (initevol == 1) THEN ! *** call timer CALL PDAF_timeit(2, 'new') ! ******************************************** ! *** Get information on next observations *** ! ******************************************** CALL U_next_observation(step_obs, nsteps, end_forecast, time) IF (debug>0) & WRITE (*,*) '++ PDAF-debug get_state: ', debug, & 'user input from next_observation: nsteps, doexit, time', & nsteps, end_forecast, time ! *** Set time step of next observation *** step_obs = step + nsteps - 1 IF (debug>0) & WRITE (*,*) '++ PDAF-debug get_state: ', debug, 'time step of next observation', & step_obs ! *** Initialize ensemble of error states for SEEK *** SEEK1: IF ((.NOT.ensemblefilter) .AND. filterpe & .AND. (subtype_filter == 0 .OR. subtype_filter ==1)) THEN DO j = 1, dim_eof DO i = 1, dim_p eofV(i, j) = state(i) + epsilon * eofV(i, j) END DO END DO END IF SEEK1 ! ********************************************************** ! *** Shift states in ensemble array in case of smoother *** ! ********************************************************** IF (nsteps > 0 .AND. dim_lag > 0) THEN CALL PDAF_smoother_shift(dim_p, dim_ens, dim_lag, eofV, sens, & cnt_maxlag, screen) END IF ! ***************************************************** ! *** Distribute ensembles and forecast information *** ! ***************************************************** doevol: IF (nsteps > 0 .AND. end_forecast /= 1) THEN ! *** More forecasts to be done *** ! *** Send sub-ensembles to each task which *** ! *** not includes the filter PEs *** IF ((mype_world == 0) .AND. (screen > 0)) THEN IF (subtype_filter == 2 .OR. subtype_filter == 3) THEN ! Output for fixed-basis (only SEEK/SEIK/LSEIK/ESTKF/LESTKF) WRITE (*, '(a, 5x, a)') 'PDAF', 'Fixed basis - evolve only state estimate' ELSE IF (ensemblefilter) THEN ! Ensemble-based filters WRITE (*, '(a, 5x, a)') 'PDAF', 'Evolve state ensemble' ELSE ! Mode-based: SEEK with evolution of basis WRITE (*, '(a, 5x, a)') 'PDAF', 'Evolve state and error space basis' END IF END IF END IF ! *** Distribute sub-ensembles to all model PEs 0 *** ! *** *** ! *** This is used if multiple model tasks are used. *** ! call timer CALL PDAF_timeit(19, 'new') CALL PDAF_timeit(49, 'new') CALL PDAF_scatter_ens(dim_p, dim_ens_l, eofV, state, screen) ! call timer CALL PDAF_timeit(19, 'old') CALL PDAF_timeit(49, 'old') END IF doevol ENSF1: IF (ensemblefilter .AND. filterpe & .AND. (subtype_filter == 2 .OR. subtype_filter == 3)) THEN ! *** Compute ensemble mean state *** state = 0.0 DO j = 1, dim_ens DO i = 1, dim_p state(i) = state(i) + eofV(i, j) END DO END DO state = state / REAL(dim_ens) ! *** Remove mean from ensemble members *** DO j = 1, dim_ens DO i = 1, dim_p eofV(i, j) = eofV(i, j) - state(i) END DO END DO END IF ENSF1 ! *** Set INIT flag *** initevol = 2 ELSE IF (initevol == 2) THEN ! Routine is called just after the first ensemble member is evolved initevol=0 END IF evolinit ! ******************************************************** ! *** Initialize state variables for ensemble forecast *** ! ******************************************************** doevol1: IF (nsteps > 0) THEN IF (ensemblefilter) THEN IF (subtype_filter/=2 .AND. subtype_filter/=3) THEN IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Evolve member ', member_get, & 'in task ', task_id ELSE IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Evolve ensemble mean state ', & 'in task ', task_id END IF ELSE IF ((task_id == statetask) .AND. (member_get == dim_eof_l + 1)) THEN IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Evolve central state ', & 'in task ', task_id ELSE IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Evolve member ', member_get, & 'in task ', task_id END IF END IF ! *** Distribute state fields within model communicators *** IF (ensemblefilter) THEN IF (subtype_filter/=2 .AND. subtype_filter/=3) THEN IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Distribute state fields', & ' in task ', task_id, ', member ', member_get ELSE IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Distribute state fields', & ' in task ', task_id, ', ensemble mean state ' END IF ELSE IF ((task_id == statetask) .AND. (member_get == dim_eof_l + 1)) THEN IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Distribute state fields', & ' in task ', task_id, ', central state ' ELSE IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' Distribute state fields', & ' in task ', task_id, ', member ', member_get END IF END IF ! *** call timer CALL PDAF_timeit(40, 'new') filtertype: IF (ensemblefilter) THEN modelpes: IF (modelpe) THEN ! Store member index for PDAF_get_memberid member_save = member_get IF (subtype_filter/=2 .AND. subtype_filter/=3) THEN ! Dynamic ensemble filter with ensemble forecast ! distribute ensemble state CALL U_distribute_state(dim_p, eofV(1:dim_p, member_get)) IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' task: ', task_id, & ' evolve sub-member ', member_get ! Increment member counter member_get = member_get + 1 ELSE ! Ensemble filter with fixed error-space basis ! (Option ONLY for SEIK/LSEIK) ! set member to maximum member_get=dim_ens_l member_put=dim_ens_l ! distribute and evolve ensemble mean state CALL U_distribute_state(dim_p, state) IF (debug > 0 .AND. modelpe .AND. mype_model==0) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' task: ', task_id, & ' evolve ensemble mean state ' END IF END IF modelpes ! Reset member counting IF (member_get == local_dim_ens+1) member_get = 1 ELSE ! Mode-based filter (SEEK) !!!! Set CENTRAL_STATE = .TRUE. for evolution of only the central !!!! state for computation of free evolution. !!!! In this case also set FREE_EVOLUTION=.TRUE. in PDAF-D_SEEK_UPDATE central_state = .FALSE. IF (central_state) THEN WRITE (*,*) 'PDAF-NOTICE: EVOLVE ONLY CENTRAL STATE FOR FREE EVOLUTION !!!' member_get = dim_eof_l + 1 member_put = dim_eof_l + 1 END IF ! For fixed basis SFEK set member to maximum IF (subtype_filter == 2 .OR. subtype_filter == 3) THEN member_get = dim_eof_l + 1 member_put = dim_eof_l + 1 END IF modelpesB: IF (modelpe) THEN ! Store member index for PDAF_get_memberid member_save = member_get IF ((task_id == statetask) .AND. (member_get == dim_eof_l + 1)) THEN ! distribute central state CALL U_distribute_state(dim_p, state) IF (debug > 0 .AND. filterpe) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' task: ',task_id, & ' evolve central state' ! Reset member counting IF (member_get == dim_eof_l+1) member_get = 1 ELSE ! distribute ensemble state CALL U_distribute_state(dim_p, eofV(1:dim_p, member_get)) IF (debug > 0 .AND. filterpe) & WRITE (*,*) '++ PDAF-debug get_state:', debug, ' task: ',task_id, & ' evolve sub-member ',member_get ! Increment member counter member_get = member_get + 1 END IF END IF modelpesB ! Reset member counting IF (task_id /= statetask .AND. member_get == dim_eof_l+1) member_get = 1 END IF filtertype ! *** call timer CALL PDAF_timeit(40, 'old') END IF doevol1 ! ******************** ! *** finishing up *** ! ******************** ! Set time step counter for next forecast phase ! (only used with PDAF_assimilate_X in fully parallel variant) cnt_steps = 0 ! Set variables for exiting the routine steps = nsteps doexit = end_forecast outflag = flag IF (debug>0) THEN WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_get_state -- END' END IF END SUBROUTINE PDAF_get_state