propagator_afte.h

Go to the documentation of this file.

      MODULE propagator_mod
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group            Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  Adjoint Finite Time Eigenvalues Propagator:                         !
!                                                                      !
!  This routine is used during the computation of the eigenvectors of  !
!  the adjoint propagator, transpose[R(t,0)]. They are computed in an  !
!  analogous way to those of  R(t,0).  A  single  integration  of  an  !
!  arbitrary perturbation state vector  "u" backward in time over the  !
!  interval [t,0] by the adjoint model: transpose[R(t,0)]*u.           !
!                                                                      !
!   Reference:                                                         !
!                                                                      !
!     Moore, A.M. et al., 2004: A comprehensive ocean prediction and   !
!       analysis system based on the tangent linear and adjoint of a   !
!       regional ocean model, Ocean Modelling, 7, 227-258.             !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
!
      implicit none
!
      PRIVATE
      PUBLIC  :: propagator_afte
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE propagator_afte (RunInterval, state, ad_state)
!***********************************************************************
!
      USE mod_param
      USE mod_parallel
#ifdef SOLVE3D
      USE mod_coupling
#endif
      USE mod_iounits
      USE mod_ocean
      USE mod_scalars
      USE mod_stepping
!
      USE close_io_mod,   ONLY : close_inp
      USE dotproduct_mod, ONLY : ad_statenorm
      USE packing_mod,    ONLY : ad_unpack, ad_pack
#ifdef SOLVE3D
      USE set_depth_mod,  ONLY : set_depth
#endif
      USE strings_mod,    ONLY : founderror
!
!  Imported variable declarations.
!
      real(dp), intent(in) :: runinterval
!
      TYPE (t_gst), intent(in) :: state(ngrids)
      TYPE (t_gst), intent(inout) :: ad_state(ngrids)
!
!  Local variable declarations.
!
#ifdef SOLVE3D
      logical :: firstpass = .true.
#endif
!
      integer :: ng, tile
!
      real(r8) :: statenorm(ngrids)
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
!=======================================================================
!  Forward integration of the tangent linear model.
!=======================================================================
!
      nrun=nrun+1
      IF (master) THEN
        DO ng=1,ngrids
          WRITE (stdout,10) ' PROPAGATOR - Grid: ', ng,                 &
     &                      ',  Iteration: ', nrun,                     &
     &                      ',  number converged RITZ values: ',        &
     &                      nconv(ng)
        END DO
      END IF
!
!  Initialize time stepping indices and counters.
!
      DO ng=1,ngrids
        iif(ng)=1
        indx1(ng)=1
        kstp(ng)=1
        krhs(ng)=3
        knew(ng)=2
        predictor_2d_step(ng)=.false.
!
        iic(ng)=0
        nstp(ng)=1
        nrhs(ng)=1
        nnew(ng)=2
!
        synchro_flag(ng)=.true.
        tdays(ng)=dstart+dt(ng)*real(ntimes(ng),r8)*sec2day
        time(ng)=tdays(ng)*day2sec
        ntstart(ng)=ntimes(ng)+1
        ntend(ng)=1
        ntfirst(ng)=ntend(ng)
      END DO
!
!-----------------------------------------------------------------------
!  Clear adjoint state variables.  There is not need to clean the basic
!  state arrays since they were zeroth out at initialization and bottom
!  of previous iteration.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, iadm)
        END DO
      END DO
 
#ifdef SOLVE3D
!
!-----------------------------------------------------------------------
!  Compute basic state initial level thicknesses used for state norm
!  scaling. It uses zero time averaged free-surface (rest state).
!  Therefore, the norm scaling is time invariant.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=last_tile(ng),first_tile(ng),-1
          CALL set_depth (ng, tile, iadm)
        END DO
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Unpack adjoint initial conditions from state vector.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ad_unpack (ng, tile, nstr(ng), nend(ng),                 &
     &                    state(ng)%vector)
        END DO
      END DO
!
!-----------------------------------------------------------------------
!  Compute initial adjoint state dot product norm.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=last_tile(ng),first_tile(ng),-1
          CALL ad_statenorm (ng, tile, knew(ng), nstp(ng),              &
     &                       statenorm(ng))
        END DO
        IF (master) THEN
          WRITE (stdout,20) ' PROPAGATOR - Grid: ', ng,                 &
     &                      ',  Adjoint Initial Norm: ', statenorm(ng)
        END IF
      END DO
!
!-----------------------------------------------------------------------
!  Read in initial forcing, climatology and assimilation data from
!  input NetCDF files.  It loads the first relevant data record for
!  the time-interpolation between snapshots.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        CALL close_inp (ng, iadm)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL ad_get_idata (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
 
        CALL ad_get_data (ng)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
      END DO
!
!-----------------------------------------------------------------------
!  Time-step the adjoint model backwards.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        IF (master) THEN
          WRITE (stdout,30) 'AD', ng, ntstart(ng), ntend(ng)
        END IF
        time(ng)=time(ng)+dt(ng)
        iic(ng)=ntstart(ng)+1
      END DO
 
#ifdef SOLVE3D
      CALL ad_main3d (runinterval)
#else
      CALL ad_main2d (runinterval)
#endif
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!-----------------------------------------------------------------------
!  Clear nonlinear state (basic state) variables and insure that the
!  time averaged free-surface is zero for scaling below and next
!  iteration.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL initialize_ocean (ng, tile, inlm)
#ifdef SOLVE3D
          CALL initialize_coupling (ng, tile, 0)
#endif
        END DO
      END DO
 
#ifdef SOLVE3D
!
!-----------------------------------------------------------------------
!  Compute basic state final level thicknesses used for state norm
!  scaling. It uses zero time averaged free-surface (rest state).
!  Therefore, the norm scaling is time invariant.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=last_tile(ng),first_tile(ng),-1
          CALL set_depth (ng, tile, iadm)
        END DO
      END DO
#endif
!
!-----------------------------------------------------------------------
!  Compute final adjoint state dot product norm.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=first_tile(ng),last_tile(ng),+1
          CALL ad_statenorm (ng, tile, knew(ng), nstp(ng),              &
     &                       statenorm(ng))
        END DO
        IF (master) THEN
          WRITE (stdout,20) ' PROPAGATOR - Grid: ', ng,                 &
     &                      ',  Adjoint   Final Norm: ', statenorm(ng)
        END IF
      END DO
!
!-----------------------------------------------------------------------
!  Pack final adjoint solution into adjoint state vector.
!-----------------------------------------------------------------------
!
      DO ng=1,ngrids
        DO tile=last_tile(ng),first_tile(ng),-1
          CALL ad_pack (ng, tile, nstr(ng), nend(ng),                   &
     &                  ad_state(ng)%vector)
        END DO
      END DO
!
 10   FORMAT (/,a,i2.2,a,i3.3,a,i3.3/)
 20   FORMAT (/,a,i2.2,a,1p,e15.6,/)
 30   FORMAT (/,1x,a,1x,'ROMS: started time-stepping:',                 &
     &        ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')')
!
      RETURN
      END SUBROUTINE propagator_afte
 
      END MODULE propagator_mod