PDAF_3dvar_analysis_cvt.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_3dvar_analysis_cvt --- 3DVAR with CVT
!
! !INTERFACE:
SUBROUTINE PDAF_3dvar_analysis_cvt(step, dim_p, dim_obs_p, dim_cvec, &
     state_p, state_inc_p, &
     U_init_dim_obs, U_obs_op, U_init_obs, U_prodRinvA, &
     U_cvt, U_cvt_adj, U_obs_op_lin, U_obs_op_adj, &
     screen, incremental, type_opt, flag)

! !DESCRIPTION:
! Analysis step of incremental 3DVAR with control
! variable transformation.
!
! Variant for domain decomposed states.
!
! !  This is a core routine of PDAF and
!    should not be changed by the user   !
!
! !REVISION HISTORY:
! 2021-03 - 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
  USE PDAF_memcounting, &
       ONLY: PDAF_memcount
  USE PDAF_mod_filtermpi, &
       ONLY: mype
  USE PDAF_mod_filter, &
       ONLY: obs_member, debug
  USE PDAFomi, &
       ONLY: omi_omit_obs => omit_obs

  IMPLICIT NONE

! !ARGUMENTS:
  INTEGER, INTENT(in) :: step         ! Current time step
  INTEGER, INTENT(in) :: dim_p        ! PE-local dimension of model state
  INTEGER, INTENT(out) :: dim_obs_p   ! PE-local dimension of observation vector
  INTEGER, INTENT(in) :: dim_cvec   ! Size of control vector
  REAL, INTENT(out)   :: state_p(dim_p)          ! on exit: PE-local forecast state
  REAL, INTENT(inout) :: state_inc_p(dim_p)      ! PE-local state analysis increment
  INTEGER, INTENT(in) :: screen       ! Verbosity flag
  INTEGER, INTENT(in) :: incremental  ! Control incremental updating
  INTEGER, INTENT(in) :: type_opt     ! Type of minimizer for 3DVar
  INTEGER, INTENT(inout) :: flag      ! Status flag

! ! External subroutines 
! ! (PDAF-internal names, real names are defined in the call to PDAF)
  EXTERNAL :: U_init_dim_obs, & ! Initialize dimension of observation vector
       U_obs_op, &              ! Observation operator
       U_init_obs, &            ! Initialize observation vector
       U_prodRinvA, &           ! Provide product R^-1 A
       U_cvt, &                 ! Apply control vector transform matrix to control vector
       U_cvt_adj, &             ! Apply adjoint control vector transform matrix
       U_obs_op_lin, &          ! Linearized observation operator
       U_obs_op_adj             ! Adjoint observation operator

! !CALLING SEQUENCE:
! Called by: PDAF_3dvar_update
! Calls: U_init_dim_obs
! Calls: U_obs_op
! Calls: U_init_obs
! Calls: PDAF_timeit
! Calls: PDAF_memcount
!EOP

! *** local variables ***
  INTEGER, SAVE :: allocflag = 0       ! Flag whether first time allocation is done
  REAL, ALLOCATABLE :: obs_p(:)        ! PE-local observation vector
  REAL, ALLOCATABLE :: dy_p(:)         ! PE-local observation background residual
  REAL, ALLOCATABLE :: v_p(:)          ! PE-local analysis increment vector
  INTEGER :: opt_parallel              ! Whether to run solver with decomposed control vector


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

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_3dvar_analysis -- START'

  IF (mype == 0 .AND. screen > 0) THEN
     IF (type_opt==1) THEN
        WRITE (*, '(a, 5x, a)') 'PDAF', '--- solver: LBFGS' 
     ELSEIF (type_opt==2) THEN
        WRITE (*, '(a, 5x, a)') 'PDAF', '--- solver: CG+' 
     ELSEIF (type_opt==3) THEN
        WRITE (*, '(a, 5x, a)') 'PDAF', '--- solver: plain CG' 
     ELSEIF (type_opt==12) THEN
        WRITE (*, '(a, 5x, a)') 'PDAF', '--- solver: CG+ parallelized' 
     ELSEIF (type_opt==13) THEN
        WRITE (*, '(a, 5x, a)') 'PDAF', '--- solver: plain CG parallelized' 
     END IF
  END IF


! *********************************
! *** Get observation dimension ***
! *********************************

  IF (debug>0) THEN
     WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, '  dim_p', dim_p
     IF (incremental<2) &
          WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_3dvar_analysis -- call init_dim_obs'
  END IF

  CALL PDAF_timeit(43, 'new')
  CALL U_init_dim_obs(step, dim_obs_p)
  CALL PDAF_timeit(43, 'old')
  
  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, '  dim_obs_p', dim_obs_p

  IF (screen > 2) THEN
     WRITE (*, '(a, 5x, a13, 1x, i6, 1x, a, i10)') &
          'PDAF', '--- PE-domain', mype, 'dimension of observation vector', dim_obs_p
  END IF


  haveobsB: IF (dim_obs_p > 0) THEN

! *******************************
! *** Background innovation   ***
! ***      d = y - H xb       ***
! *******************************

     CALL PDAF_timeit(12, 'new')
  
     ! *** Observation background innovation ***

     ! Get observed state estimate
     ALLOCATE(dy_p(dim_obs_p))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_p)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_3dvar_analysis -- call obs_op'

     obs_member = 0 ! Store member index (0 for central state)
     CALL PDAF_timeit(44, 'new')
     CALL U_obs_op(step, dim_p, dim_obs_p, state_p, dy_p)
     CALL PDAF_timeit(44, 'old')

     ! get observation vector
     ALLOCATE(obs_p(dim_obs_p))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_p)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_3dvar_analysis -- call init_obs'

     CALL PDAF_timeit(50, 'new')
     CALL U_init_obs(step, dim_obs_p, obs_p)
     CALL PDAF_timeit(50, 'old')

     ! Get residual as difference of observation and observed state
     CALL PDAF_timeit(51, 'new')
     dy_p = obs_p - dy_p
     CALL PDAF_timeit(51, 'old')

     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, &
             'innovation d(1:min(dim_obs_p,6))', dy_p(1:min(dim_p,6))
        WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, &
             'MIN/MAX of innovation', MINVAL(dy_p), MAXVAL(dy_p)
     END IF

     ! Omit observations with too high innovation
     IF (omi_omit_obs)  THEN
        CALL PDAF_timeit(51, 'new')
        CALL PDAFomi_omit_by_inno_cb(dim_obs_p, dy_p, obs_p)
        CALL PDAF_timeit(51, 'old')
     END IF

     CALL PDAF_timeit(12, 'old')


! ****************************
! ***   Iterative solving  ***
! ****************************

     CALL PDAF_timeit(52, 'new')

     opt_parallel = 0

     ! Prepare control vector for optimization
     ALLOCATE(v_p(dim_cvec))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_cvec)
     v_p = 0.0

     ! Choose solver
     opt: IF (type_opt==1) THEN

        ! LBFGS solver
        CALL PDAF_3dvar_optim_lbfgs(step, dim_p, dim_cvec, dim_obs_p, &
             obs_p, dy_p, v_p, U_prodRinvA, &
             U_cvt, U_cvt_adj, U_obs_op_lin, U_obs_op_adj, &
             opt_parallel, screen)

     ELSEIF (type_opt==2 .OR. type_opt==12) THEN

        IF (type_opt==12) opt_parallel = 1

        ! CG+ solver
        CALL PDAF_3dvar_optim_cgplus(step, dim_p, dim_cvec, dim_obs_p, &
             obs_p, dy_p, v_p, U_prodRinvA, &
             U_cvt, U_cvt_adj, U_obs_op_lin, U_obs_op_adj, &
             opt_parallel, screen)

     ELSEIF (type_opt==3 .OR. type_opt==13) THEN

        IF (type_opt==13) opt_parallel = 1

        ! CG solver
        CALL PDAF_3dvar_optim_cg(step, dim_p, dim_cvec, dim_obs_p, &
             obs_p, dy_p, v_p, U_prodRinvA, &
             U_cvt, U_cvt_adj, U_obs_op_lin, U_obs_op_adj, &
             opt_parallel, screen)

     ELSE
        ! Further solvers - not implemented
        WRITE (*,'(a)') 'PDAF-ERROR: Invalid solver chosen!'
        flag=10
     END IF opt

     CALL PDAF_timeit(52, 'old')


! ****************************************************
! ***   Solving completed: Compute state increment ***
! ****************************************************

     CALL PDAF_timeit(14, 'new')

     ! State increment: Apply V to control vector v_p
     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, &
             'control vector (1:min(dim_cvec,6))', v_p(1:min(dim_cvec,6))
        WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, &
             'MIN/MAX of control vector', MINVAL(v_p), MAXVAL(v_p)
        WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_3dvar_analysis -- call cvt for final state increment'
     END IF

     CALL PDAF_timeit(49, 'new')
     CALL U_cvt(0, dim_p, dim_cvec, v_p, state_inc_p)
     CALL PDAF_timeit(49, 'old')

     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, &
             'state vector increment (1:min(dim_cvec,p,6))', state_inc_p(1:min(dim_p,6))
        WRITE (*,*) '++ PDAF-debug PDAF_3dvar_analysis:', debug, &
             'MIN/MAX of state vector increment', MINVAL(state_inc_p), MAXVAL(state_inc_p)
     END IF

     CALL PDAF_timeit(51, 'new')
     IF (incremental<1) THEN
        ! Add analysis increment to state vector
        state_p = state_p + state_inc_p
     END IF
     CALL PDAF_timeit(51, 'old')

     CALL PDAF_timeit(14, 'old')

  END IF haveobsB


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

  IF (dim_obs_p > 0) THEN
     DEALLOCATE(obs_p, dy_p)
     DEALLOCATE(v_p)
  END IF

  IF (allocflag == 0) allocflag = 1

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_3dvar_analysis -- END'

END SUBROUTINE PDAF_3dvar_analysis_cvt