adsen_force.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE adsen_force_mod
 
#if defined ADJOINT          && \
   (defined ad_sensitivity   || defined i4dvar_ana_sensitivity || \
    defined opt_observations || defined sensitivity_4dvar)
!
!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                                               !
!=======================================================================
!                                                                      !
!  This routine forces the adjoint state with the functional whose     !
!  sensitivity is required.                                            !
!                                                                      !
!=======================================================================
!
      implicit none
 
      PRIVATE
      PUBLIC  :: adsen_force
 
      CONTAINS
!
!***********************************************************************
      SUBROUTINE adsen_force (ng, tile)
!***********************************************************************
!
      USE mod_param
      USE mod_clima
# ifdef SOLVE3D
      USE mod_coupling
# endif
      USE mod_grid
      USE mod_ocean
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
!
!  Local variable declarations.
!
# include "tile.h"
!
      CALL adsen_force_tile (ng, tile,                                  &
     &                       lbi, ubi, lbj, ubj,                        &
     &                       imins, imaxs, jmins, jmaxs,                &
     &                       knew(ng),                                  &
# ifdef SOLVE3D
     &                       nnew(ng), nstp(ng),                        &
# endif
     &                       grid(ng) % Rscope,                         &
     &                       grid(ng) % Uscope,                         &
     &                       grid(ng) % Vscope,                         &
# ifdef SOLVE3D
     &                       clima(ng) % u_ads,                         &
     &                       clima(ng) % v_ads,                         &
     &                       clima(ng) % wvel_ads,                      &
     &                       clima(ng) % t_ads,                         &
# endif
     &                       clima(ng) % ubar_ads,                      &
     &                       clima(ng) % vbar_ads,                      &
     &                       clima(ng) % zeta_ads,                      &
# ifdef SOLVE3D
     &                       ocean(ng) % ad_u,                          &
     &                       ocean(ng) % ad_v,                          &
     &                       ocean(ng) % ad_wvel,                       &
     &                       ocean(ng) % ad_t,                          &
     &                       coupling(ng) % ad_Zt_avg1,                 &
# else
     &                       ocean(ng) % ad_zeta,                       &
# endif
     &                       ocean(ng) % ad_ubar,                       &
     &                       ocean(ng) % ad_vbar)
 
      RETURN
      END SUBROUTINE adsen_force
!
!***********************************************************************
      SUBROUTINE adsen_force_tile (ng, tile,                            &
     &                             LBi, UBi, LBj, UBj,                  &
     &                             IminS, ImaxS, JminS, JmaxS,          &
     &                             knew,                                &
# ifdef SOLVE3D
     &                             nnew, nstp,                          &
# endif
     &                             Rscope, Uscope, Vscope,              &
# ifdef SOLVE3D
     &                             u_ads, v_ads, wvel_ads,              &
     &                             t_ads,                               &
# endif
     &                             ubar_ads, vbar_ads, zeta_ads,        &
# ifdef SOLVE3D
     &                             ad_u, ad_v, ad_wvel,                 &
     &                             ad_t, ad_Zt_avg1,                    &
# else
     &                             ad_zeta,                             &
# endif
     &                             ad_ubar, ad_vbar)
!***********************************************************************
!
      USE mod_param
      USE mod_parallel
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: knew
# ifdef SOLVE3D
      integer, intent(in) :: nnew, nstp
# endif
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: Rscope(LBi:,LBj:)
      real(r8), intent(in) :: Uscope(LBi:,LBj:)
      real(r8), intent(in) :: Vscope(LBi:,LBj:)
#  ifdef SOLVE3D
      real(r8), intent(in) :: u_ads(LBi:,LBj:,:)
      real(r8), intent(in) :: v_ads(LBi:,LBj:,:)
      real(r8), intent(in) :: wvel_ads(LBi:,LBj:,:)
      real(r8), intent(in) :: t_ads(LBi:,LBj:,:,:)
#  endif
      real(r8), intent(in) :: ubar_ads(LBi:,LBj:)
      real(r8), intent(in) :: vbar_ads(LBi:,LBj:)
      real(r8), intent(in) :: zeta_ads(LBi:,LBj:)
#  ifdef SOLVE3D
      real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: ad_wvel(LBi:,LBj:,:)
      real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
      real(r8), intent(inout) :: ad_Zt_avg1(LBi:,LBj:)
#  else
      real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
#  endif
      real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
      real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
# else
      real(r8), intent(in) :: Rscope(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Uscope(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: Vscope(LBi:UBi,LBj:UBj)
#  ifdef SOLVE3D
      real(r8), intent(in) :: u_ads(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: v_ads(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: wvel_ads(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(in) :: t_ads(LBi:UBi,LBj:UBj,N(ng),NT(ng))
#  endif
      real(r8), intent(in) :: ubar_ads(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: vbar_ads(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: zeta_ads(LBi:UBi,LBj:UBj)
#  ifdef SOLVE3D
      real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
      real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
      real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,0:N(ng))
      real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
      real(r8), intent(inout) :: ad_Zt_avg1(LBi:UBi,LBj:UBj)
#  else
      real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
#  endif
      real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
      real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
# endif
!
!  Local variable declarations.
!
      integer :: Kfrc, Nfrc, i, itrc, j, k
 
      real(r8) :: adFac
 
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Initialize adjoint staye with the functional whose sensitivity is
!  required.  Use functional loaded into first record of climatological
!  arrays.
!-----------------------------------------------------------------------
!
      IF (iic(ng).eq.ntend(ng)) THEN
        kfrc=knew
        nfrc=nstp
      ELSE
        kfrc=1
        nfrc=nnew
      END IF
 
# ifdef AD_IMPULSE
!
!  Impose the forcing in the adjoint model as impulses in a manner that
!  is consistent with the definition of the sensitivity functional.
!  Notice that nTLM is used to control the forcing since nADJ is
!  used to control strong versus weak constraint so cannot be
!  arbitrarily defined.
!
      adfac=0.0_r8
#  ifdef I4DVAR_ANA_SENSITIVITY
      IF ((mod(iic(ng)-1,nhis(ng)).eq.0).and.                           &
     &    (dends(ng).ge.tdays(ng)).and.(tdays(ng).ge.dstrs(ng))) THEN
#  else
      IF ((mod(iic(ng)-1,ntlm(ng)).eq.0).and.                           &
     &    (dends(ng).ge.tdays(ng)).and.(tdays(ng).ge.dstrs(ng))) THEN
#  endif
          adfac=1.0_r8
          IF (master) THEN
            WRITE (stdout,10) iic(ng)-1
 10         FORMAT (2x,'ADSEN_FORCE      - forcing Adjoint model at',   &
     &                 ' TimeStep: ', i0)
          END IF
      END IF
# else
      adfac=1.0_r8
# endif
!
!  Free-surface.
!
      IF (scalars(ng)%Lstate(isfsur)) THEN
# ifdef SOLVE3D
        DO j=jstrr,jendr
          DO i=istrr,iendr
            ad_zt_avg1(i,j)=ad_zt_avg1(i,j)+                            &
     &                      adfac*zeta_ads(i,j)*rscope(i,j)
          END DO
        END DO
# else
        DO j=jstrr,jendr
          DO i=istrr,iendr
            ad_zeta(i,j,kfrc)=ad_zeta(i,j,kfrc)+                        &
     &                        adfac*zeta_ads(i,j)*rscope(i,j)
          END DO
        END DO
# endif
      END IF
!
!  2D Momentum.
!
      IF (scalars(ng)%Lstate(isubar)) THEN
        DO j=jstrr,jendr
          DO i=istr,iendr
            ad_ubar(i,j,kfrc)=ad_ubar(i,j,kfrc)+                        &
     &                        adfac*ubar_ads(i,j)*uscope(i,j)
          END DO
        END DO
      END IF
!
      IF (scalars(ng)%Lstate(isvbar)) THEN
        DO j=jstr,jendr
          DO i=istrr,iendr
            ad_vbar(i,j,kfrc)=ad_vbar(i,j,kfrc)+                        &
     &                        adfac*vbar_ads(i,j)*vscope(i,j)
          END DO
        END DO
      END IF
# ifdef SOLVE3D
!
!  3D Momentum.
!
      IF (scalars(ng)%Lstate(isuvel)) THEN
        DO k=kstrs(ng),kends(ng)
          DO j=jstrr,jendr
            DO i=istr,iendr
              ad_u(i,j,k,nfrc)=ad_u(i,j,k,nfrc)+                        &
     &                         adfac*u_ads(i,j,k)*uscope(i,j)
            END DO
          END DO
        END DO
      END IF
!
      IF (scalars(ng)%Lstate(isvvel)) THEN
        DO k=kstrs(ng),kends(ng)
          DO j=jstr,jendr
            DO i=istrr,iendr
              ad_v(i,j,k,nfrc)=ad_v(i,j,k,nfrc)+                        &
     &                         adfac*v_ads(i,j,k)*vscope(i,j)
            END DO
          END DO
        END DO
      END IF
!
!  Vertical velocity.
!
!  Notice that vertical velocity will not be forced at the top in the
!  current code. To do this, uncomment the next line and ensure that
!  KendSb=N in the ocean input file.
!
      IF (scalars(ng)%Lstate(isvvel)) THEN
!!      DO k=KstrS(ng),KendS(ng)+1             ! if forced at the top
        DO k=kstrs(ng),kends(ng)
          DO j=jstrr,jendr
            DO i=istrr,iendr
              ad_wvel(i,j,k)=ad_wvel(i,j,k)+                            &
     &                       adfac*wvel_ads(i,j,k)*rscope(i,j)
            END DO
          END DO
        END DO
      END IF
!
!  Tracers.
!
      DO itrc=1,nt(ng)
        IF (scalars(ng)%Lstate(istvar(itrc))) THEN
          DO k=kstrs(ng),kends(ng)
            DO j=jstrr,jendr
              DO i=istrr,iendr
                ad_t(i,j,k,nfrc,itrc)=ad_t(i,j,k,nfrc,itrc)+            &
     &                                adfac*t_ads(i,j,k,itrc)*          &
     &                                rscope(i,j)
              END DO
            END DO
          END DO
        END IF
      END DO
# endif
 
      RETURN
      END SUBROUTINE adsen_force_tile
#endif
      END MODULE adsen_force_mod