ad_misfit.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE ad_misfit_mod
 
#if  defined ADJOINT && defined FOUR_DVAR  && defined OBSERVATIONS && \
    !defined WEAK_CONSTRAINT
!
!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 computes the model minus observations adjoint misfit   !
!  forcing for each state variable.                                    !
!                                                                      !
!  The observation screening and quality control variable "ObsScale"   !
!  is not modified in this routine. Their values are the ones set in   !
!  obs_write.F during the running of the nonlinear model.              !
!                                                                      !
!=======================================================================
!
      implicit none
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE ad_misfit (ng, tile, model)
!***********************************************************************
!
      USE mod_param
      USE mod_grid
      USE mod_ocean
      USE mod_stepping
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
!
!  Local variable declarations.
!
# include "tile.h"
!
      CALL ad_misfit_tile (ng, tile, model,                             &
     &                     lbi, ubi, lbj, ubj,                          &
     &                     imins, imaxs, jmins, jmaxs,                  &
     &                     kstp(ng),                                    &
# ifdef SOLVE3D
     &                     nstp(ng),                                    &
# endif
# ifdef MASKING
     &                     grid(ng) % rmask,                            &
     &                     grid(ng) % umask,                            &
     &                     grid(ng) % vmask,                            &
# endif
# ifdef SOLVE3D
     &                     grid(ng) % z_r,                              &
     &                     grid(ng) % z_v,                              &
     &                     ocean(ng) % ad_u,                            &
     &                     ocean(ng) % ad_v,                            &
     &                     ocean(ng) % ad_t,                            &
# endif
     &                     ocean(ng) % ad_ubar,                         &
     &                     ocean(ng) % ad_vbar,                         &
     &                     ocean(ng) % ad_zeta)
      RETURN
      END SUBROUTINE ad_misfit
!
!***********************************************************************
      SUBROUTINE ad_misfit_tile (ng, tile, model,                       &
     &                           LBi, UBi, LBj, UBj,                    &
     &                           IminS, ImaxS, JminS, JmaxS,            &
     &                           kstp,                                  &
# ifdef SOLVE3D
     &                           nstp,                                  &
# endif
# ifdef MASKING
     &                           rmask, umask, vmask,                   &
# endif
# ifdef SOLVE3D
     &                           z_r, z_v,                              &
     &                           ad_u, ad_v, ad_t,                      &
# endif
     &                           ad_ubar, ad_vbar, ad_zeta)
!***********************************************************************
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
# ifdef DISTRIBUTE
      USE distribute_mod, ONLY : mp_collect
# endif
      USE ad_extract_obs_mod, ONLY : ad_extract_obs2d
# ifdef SOLVE3D
      USE ad_extract_obs_mod, ONLY : ad_extract_obs3d
# endif
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
      integer, intent(in) :: LBi, UBi, LBj, UBj
      integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
      integer, intent(in) :: kstp
# ifdef SOLVE3D
      integer, intent(in) :: nstp
# endif
!
# ifdef ASSUMED_SHAPE
#  ifdef MASKING
      real(r8), intent(in) :: rmask(LBi:,LBj:)
      real(r8), intent(in) :: umask(LBi:,LBj:)
      real(r8), intent(in) :: vmask(LBi:,LBj:)
#  endif
#  ifdef SOLVE3D
      real(r8), intent(in) :: z_r(LBi:,LBj:,:)
      real(r8), intent(inout) :: z_v(LBi:,LBj:,:)
      real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
      real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
#  endif
      real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
      real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
      real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
# else
#  ifdef MASKING
      real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
#  endif
#  ifdef SOLVE3D
      real(r8), intent(in) :: z_r(LBi:UBi,LBj:UBj,N(ng))
      real(r8), intent(inout) :: z_v(LBi:UBi,LBj:UBj,N(ng))
      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_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
#  endif
      real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
      real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
      real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
# endif
!
!  Local variable declarations.
!
      integer :: ObsSum, ObsVoid
      integer :: i, ic, ie, iobs, is, kfrc
# ifdef SOLVE3D
      integer :: itrc, j, k, nfrc
# endif
      real(r8) :: angle
 
      real(r8), parameter :: IniVal = 0.0_r8
 
      real(r8) :: ad_uradial(Mobs), ad_vradial(Mobs)
 
# include "set_bounds.h"
!
!=======================================================================
!  Compute model minus observations adjoint misfit forcing.
# ifdef DISTRIBUTE
!  Notice that ad_mp_exchange*d calls are not required here
!  (see weak constraint routine ad_htobs.F) because we are
!  using the adjoint state arrays directly which will be
!  exchanged during the initialization step.
# endif
!
!  The processing flag used to reject (ObsVetting=0) or accept
!  (ObsVetting=1) observations is computed here but it is never
!  used. The observation screening and quality control variable
!  (ObsScale) is only computed in routine obs_write.
!=======================================================================
!
      IF (processobs(ng)) THEN
!
!  Set adjoint time index to update.
!
        kfrc=kstp
# ifdef SOLVE3D
        nfrc=nstp
# endif
!
!  Initialize observation reject/accept processing flag.
!
        DO iobs=1,mobs
          obsvetting(iobs)=inival
        END DO
!
!  Compute adjoint forcing terms at observation locations.
!
        DO iobs=1,nobs(ng)
          admodval(iobs)=obserr(iobs)*                                  &
     &                   (nlmodval(iobs)+tlmodval(iobs)-obsval(iobs))
        END DO
!
!-----------------------------------------------------------------------
!  Load observation reject and accept flag into screening variable.
!-----------------------------------------------------------------------
!
!  In the primal formulation, the current time window (or survey)
!  observations are loaded to the working arrays using local indexes
!  (array elements 1:Nobs) as opposed to global indexes in the dual
!  formulation. Recall that the screening variable ObsScale is computed
!  only once to facilitate Background Quality Control on the first pass
!  of the NLM. Therefore, we need to load its values from the saved
!  global arrays values, ObsScaleGlobal.
!
        ic=0
        DO iobs=nstrobs(ng),nendobs(ng)
          ic=ic+1
          obsscale(ic)=obsscaleglobal(iobs)
        END DO
 
# ifdef BGQC
!
!-----------------------------------------------------------------------
!  Reject observations that fail background quality control check.
!-----------------------------------------------------------------------
!
        DO iobs=1,nobs(ng)
          admodval(iobs)=obsscale(iobs)*admodval(iobs)
        END DO
# endif
!
!-----------------------------------------------------------------------
!  Free-surface.
!-----------------------------------------------------------------------
!
        IF (fourdvar(ng)%ObsCount(isfsur).gt.0) THEN
          CALL ad_extract_obs2d (ng, 0, lm(ng)+1, 0, mm(ng)+1,          &
     &                           lbi, ubi, lbj, ubj,                    &
     &                           obsstate2type(isfsur),                 &
     &                           mobs, 1, nobs(ng),                     &
     &                           rxmin(ng), rxmax(ng),                  &
     &                           rymin(ng), rymax(ng),                  &
     &                           time(ng), dt(ng),                      &
     &                           obstype, obsvetting,                   &
     &                           tobs, xobs, yobs,                      &
     &                           ad_zeta(:,:,kfrc),                     &
# ifdef MASKING
     &                           rmask,                                 &
# endif
     &                           admodval)
        END IF
!
!-----------------------------------------------------------------------
!  2D u-momentum component.
!-----------------------------------------------------------------------
!
        IF (fourdvar(ng)%ObsCount(isubar).gt.0) THEN
          CALL ad_extract_obs2d (ng, 1, lm(ng)+1, 0, mm(ng)+1,          &
     &                           lbi, ubi, lbj, ubj,                    &
     &                           obsstate2type(isubar),                 &
     &                           mobs, 1, nobs(ng),                     &
     &                           uxmin(ng), uxmax(ng),                  &
     &                           uymin(ng), uymax(ng),                  &
     &                           time(ng), dt(ng),                      &
     &                           obstype, obsvetting,                   &
     &                           tobs, xobs, yobs,                      &
     &                           ad_ubar(:,:,kfrc),                     &
# ifdef MASKING
     &                           umask,                                 &
# endif
     &                           admodval)
        END IF
!
!-----------------------------------------------------------------------
!  2D v-momentum component.
!-----------------------------------------------------------------------
!
        IF (fourdvar(ng)%ObsCount(isvbar).gt.0) THEN
          CALL ad_extract_obs2d (ng, 0, lm(ng)+1, 1, mm(ng)+1,          &
     &                           lbi, ubi, lbj, ubj,                    &
     &                           obsstate2type(isvbar),                 &
     &                           mobs, 1, nobs(ng),                     &
     &                           vxmin(ng), vxmax(ng),                  &
     &                           vymin(ng), vymax(ng),                  &
     &                           time(ng), dt(ng),                      &
     &                           obstype, obsvetting,                   &
     &                           tobs, xobs, yobs,                      &
     &                           ad_vbar(:,:,kfrc),                     &
# ifdef MASKING
     &                           vmask,                                 &
# endif
     &                           admodval)
        END IF
 
# ifdef SOLVE3D
!
!-----------------------------------------------------------------------
!  3D u-momentum component.
!-----------------------------------------------------------------------
!
        IF (fourdvar(ng)%ObsCount(isuvel).gt.0) THEN
          DO k=1,n(ng)
            DO j=jstr-1,jend+1
              DO i=istru-1,iend+1
                z_v(i,j,k)=0.5_r8*(z_r(i-1,j,k)+                        &
     &                             z_r(i  ,j,k))
              END DO
            END DO
          END DO
          CALL ad_extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,          &
     &                           lbi, ubi, lbj, ubj, 1, n(ng),          &
     &                           obsstate2type(isuvel),                 &
     &                           mobs, 1, nobs(ng),                     &
     &                           uxmin(ng), uxmax(ng),                  &
     &                           uymin(ng), uymax(ng),                  &
     &                           time(ng), dt(ng),                      &
     &                           obstype, obsvetting,                   &
     &                           tobs, xobs, yobs, zobs,                &
     &                           ad_u(:,:,:,nfrc),                      &
     &                           z_v,                                   &
#  ifdef MASKING
     &                           umask,                                 &
#  endif
     &                           admodval)
        END IF
!
!-----------------------------------------------------------------------
!  3D v-momentum component.
!-----------------------------------------------------------------------
!
        IF (fourdvar(ng)%ObsCount(isvvel).gt.0) THEN
          DO k=1,n(ng)
            DO j=jstrv-1,jend+1
              DO i=istr-1,iend+1
                z_v(i,j,k)=0.5_r8*(z_r(i,j-1,k)+                        &
     &                             z_r(i,j  ,k))
              END DO
            END DO
          END DO
          CALL ad_extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,          &
     &                           lbi, ubi, lbj, ubj, 1, n(ng),          &
     &                           obsstate2type(isvvel),                 &
     &                           mobs, 1, nobs(ng),                     &
     &                           vxmin(ng), vxmax(ng),                  &
     &                           vymin(ng), vymax(ng),                  &
     &                           time(ng), dt(ng),                      &
     &                           obstype, obsvetting,                   &
     &                           tobs, xobs, yobs, zobs,                &
     &                           ad_v(:,:,:,nfrc),                      &
     &                           z_v,                                   &
#  ifdef MASKING
     &                           vmask,                                 &
#  endif
     &                           admodval)
        END IF
!
!-----------------------------------------------------------------------
!  Radial Velocity.  The observations are in terms of radial speed and
!  angle (stored in obs_meta).  The observation angle converts the
!  velocity components to geographical EAST and North components.
#  ifdef RADIAL_ANGLE_CCW_EAST
!  The radial velocity observations are processed as magnitude and
!  heading angle (obs_meta; radians) in the math convention: an
!  azimuth that is counterclockwise from TRUE East.
!
!  In curvilinear coordinates, the radial forward problem is:
!
!     radial = u * COS(obs_meta - angler) + v * SIN(obs_meta - angler)
!
!  In the adjoint, we get:
!
!        ad_vradial = ad_vradial + ADmodVal * SIN(obs_meta - angler)
!        ad_uradial = ad_uradial + ADmodVal * COS(obs_meta - angler)
#  else
!  By default, the radial velocity observations are processed as
!  magnitude and heading angle (obs_meta; radians) in the navigation
!  convention: an azimuth that is clockwise from TRUE North.
!
!  In curvilinear coordinates, the radial forward problem is:
!
!     radial = u * SIN(obs_meta + angler) + v * COS(obs_meta + angler)
!
!  In the adjoint, we get:
!
!        ad_vradial = ad_vradial + ADmodVal * COS(obs_meta + angler)
!        ad_uradial = ad_uradial + ADmodVal * SIN(obs_meta + angler)
#  endif
!-----------------------------------------------------------------------
!
        IF (fourdvar(ng)%ObsCount(isradial).gt.0) THEN
          DO iobs=1,nobs(ng)
            ad_uradial(iobs)=inival
            ad_vradial(iobs)=inival
          END DO
          DO iobs=1,nobs(ng)
            IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
#  ifdef RADIAL_ANGLE_CCW_EAST
#   ifdef CURVGRID
              angle=obsmeta(iobs)-obsangler(iobs)
              ad_uradial(iobs)=ad_uradial(iobs)+                        &
     &                         admodval(iobs)*cos(angle)
              ad_vradial(iobs)=ad_vradial(iobs)+                        &
     &                         admodval(iobs)*sin(angle)
#   else
              ad_uradial(iobs)=ad_uradial(iobs)+                        &
     &                         admodval(iobs)*cos(obsmeta(iobs))
              ad_vradial(iobs)=ad_vradial(iobs)+                        &
     &                         admodval(iobs)*sin(obsmeta(iobs))
#   endif
#  else
#   ifdef CURVGRID
              angle=obsmeta(iobs)+obsangler(iobs)
              ad_uradial(iobs)=ad_uradial(iobs)+                        &
     &                         admodval(iobs)*sin(angle)
              ad_vradial(iobs)=ad_vradial(iobs)+                        &
     &                         admodval(iobs)*cos(angle)
#   else
              ad_uradial(iobs)=ad_uradial(iobs)+                        &
     &                         admodval(iobs)*sin(obsmeta(iobs))
              ad_vradial(iobs)=ad_vradial(iobs)+                        &
     &                         admodval(iobs)*cos(obsmeta(iobs))
#   endif
#  endif
            END IF
          END DO
          DO k=1,n(ng)
            DO j=jstrv-1,jend+1
              DO i=istr-1,iend+1
                z_v(i,j,k)=0.5_r8*(z_r(i,j-1,k)+                        &
     &                             z_r(i,j  ,k))
              END DO
            END DO
          END DO
          CALL ad_extract_obs3d (ng, 0, lm(ng)+1, 1, mm(ng)+1,          &
     &                           lbi, ubi, lbj, ubj, 1, n(ng),          &
     &                           obsstate2type(isradial),               &
     &                           mobs, 1, nobs(ng),                     &
     &                           vxmin(ng), vxmax(ng),                  &
     &                           vymin(ng), vymax(ng),                  &
     &                           time(ng), dt(ng),                      &
     &                           obstype, obsvetting,                   &
     &                           tobs, xobs, yobs, zobs,                &
     &                           ad_v(:,:,:,nfrc), z_v,                 &
#  ifdef MASKING
     &                           vmask,                                 &
#  endif
     &                           ad_vradial)
          DO k=1,n(ng)
            DO j=jstr-1,jend+1
              DO i=istru-1,iend+1
                z_v(i,j,k)=0.5_r8*(z_r(i-1,j,k)+                        &
     &                             z_r(i  ,j,k))
              END DO
            END DO
          END DO
          CALL ad_extract_obs3d (ng, 1, lm(ng)+1, 0, mm(ng)+1,          &
     &                           lbi, ubi, lbj, ubj, 1, n(ng),          &
     &                           obsstate2type(isradial),               &
     &                           mobs, 1, nobs(ng),                     &
     &                           uxmin(ng), uxmax(ng),                  &
     &                           uymin(ng), uymax(ng),                  &
     &                           time(ng), dt(ng),                      &
     &                           obstype, obsvetting,                   &
     &                           tobs, xobs, yobs, zobs,                &
     &                           ad_u(:,:,:,nfrc), z_v,                 &
#  ifdef MASKING
     &                           umask,                                 &
#  endif
     &                           ad_uradial)
        END IF
!
!-----------------------------------------------------------------------
!  Tracer type variables.
!-----------------------------------------------------------------------
!
        DO itrc=1,nt(ng)
          IF (fourdvar(ng)%ObsCount(istvar(itrc)).gt.0) THEN
            CALL ad_extract_obs3d (ng, 0, lm(ng)+1, 0, mm(ng)+1,        &
     &                             lbi, ubi, lbj, ubj, 1, n(ng),        &
     &                             obsstate2type(istvar(itrc)),         &
     &                             mobs, 1, nobs(ng),                   &
     &                             rxmin(ng), rxmax(ng),                &
     &                             rymin(ng), rymax(ng),                &
     &                             time(ng), dt(ng),                    &
     &                             obstype, obsvetting,                 &
     &                             tobs, xobs, yobs, zobs,              &
     &                             ad_t(:,:,:,nfrc,itrc),               &
     &                             z_r,                                 &
#  ifdef MASKING
     &                             rmask,                               &
#  endif
     &                             admodval)
          END IF
        END DO
# endif
# ifdef DISTRIBUTE
!
!-----------------------------------------------------------------------
!  For debugging purposes, collect all observations reject/accept
!  processing flag.
!-----------------------------------------------------------------------
!
        CALL mp_collect (ng, model, mobs, inival, obsvetting)
# endif
!
!-----------------------------------------------------------------------
!  Set counters for the number of rejected observations for each state
!  variable. Although unnecessary, the counters are recomputed here to
!  check if "ObsScale" changed from its initial values.
!-----------------------------------------------------------------------
!
        DO iobs=1,nobs(ng)
          IF (obsscale(iobs).lt.1.0) THEN
            IF (obstype(iobs).eq.obsstate2type(isfsur)) THEN
              fourdvar(ng)%ObsReject(isfsur)=                           &
     &                              fourdvar(ng)%ObsReject(isfsur)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isubar)) THEN
              fourdvar(ng)%ObsReject(isubar)=                           &
     &                              fourdvar(ng)%ObsReject(isubar)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isvbar)) THEN
              fourdvar(ng)%ObsReject(isvbar)=                           &
     &                              fourdvar(ng)%ObsReject(isvbar)+1
# ifdef SOLVE3D
            ELSE IF (obstype(iobs).eq.obsstate2type(isuvel)) THEN
              fourdvar(ng)%ObsReject(isuvel)=                           &
     &                              fourdvar(ng)%ObsReject(isuvel)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isvvel)) THEN
              fourdvar(ng)%ObsReject(isvvel)=                           &
     &                              fourdvar(ng)%ObsReject(isvvel)+1
            ELSE IF (obstype(iobs).eq.obsstate2type(isradial)) THEN
              fourdvar(ng)%ObsReject(isradial)=                         &
     &                              fourdvar(ng)%ObsReject(isradial)+1
            ELSE
              DO itrc=1,nt(ng)
                IF (obstype(iobs).eq.obsstate2type(istvar(itrc))) THEN
                  i=istvar(itrc)
                  fourdvar(ng)%ObsReject(i)=fourdvar(ng)%ObsReject(i)+1
                END IF
              END DO
# endif
            END IF
          END IF
        END DO
!
!  Load total available and rejected observations into structure
!  array.
!
        DO i=1,nobsvar(ng)
          fourdvar(ng)%ObsCount(0)=fourdvar(ng)%ObsCount(0)+            &
     &                             fourdvar(ng)%ObsCount(i)
          fourdvar(ng)%ObsReject(0)=fourdvar(ng)%ObsReject(0)+          &
     &                              fourdvar(ng)%ObsReject(i)
        END DO
!
!  Report.
!
        IF (domain(ng)%SouthWest_Test(tile)) THEN
          IF (master) THEN
            obssum=0
            obsvoid=0
            is=nstrobs(ng)
            DO i=1,nobsvar(ng)
              IF (fourdvar(ng)%ObsCount(i).gt.0) THEN
                ie=is+fourdvar(ng)%ObsCount(i)-1
                WRITE (stdout,10) trim(obsname(i)), is, ie,             &
     &                            ie-is+1, fourdvar(ng)%ObsReject(i)
                is=ie+1
                obssum=obssum+fourdvar(ng)%ObsCount(i)
                obsvoid=obsvoid+fourdvar(ng)%ObsReject(i)
              END IF
            END DO
            WRITE (stdout,20) obssum, obsvoid,                          &
     &                        fourdvar(ng)%ObsCount(0),                 &
     &                        fourdvar(ng)%ObsReject(0)
            WRITE (stdout,30) time_code(ng), nstrobs(ng), nendobs(ng),  &
     &                        iic(ng)
  10        FORMAT (10x,a,t25,4(1x,i10))
  20        FORMAT (/,10x,'Total',t47,2(1x,i10),                        &
     &              /,10x,'Obs Tally',t47,2(1x,i10),/)
  30        FORMAT (2x,' AD_MISFIT       - Added observations misfit ', &
     &              'forcing,',t75,a,/,22x,'(Observation ',             &
     &              'records = ',i0,' - ',i0,', iic = ',i0,')')
          END IF
        END IF
      END IF
      RETURN
      END SUBROUTINE ad_misfit_tile
#endif
      END MODULE ad_misfit_mod