ad_extract_obs.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE ad_extract_obs_mod
#if defined ADJOINT && defined FOUR_DVAR && defined OBSERVATIONS
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This subroutine computes the adjoint misfit forcing terms at the    !
!  observations locations (Xobs,Yobs,Zobs).  The extraction is done    !
!  via linear interpolation.  The (Xobs,Yobs) locations  must be in    !
!  grid coordinates. Zobs can be in grid coordinates (Zobs >= 0) or    !
!  actual depths (Zobs < 0), if applicable.                            !
!                                                                      !
!  On Input:                                                           !
!                                                                      !
!     ng         Nested grid number.                                   !
!     Imin       Global I-coordinate lower bound threshold for         !
!                  requested state field type.                         !
!     Imax       Global I-coordinate upper bound threshold for         !
!                  requested state field type.                         !
!     Jmin       Global J-coordinate lower bound threshold for         !
!                  requested state field type.                         !
!     Jmax       Global J-coordinate upper bound threshold for         !
!                  requested state field type.                         !
!     LBi        I-dimension Lower bound.                              !
!     UBi        I-dimension Upper bound.                              !
!     LBj        J-dimension Lower bound.                              !
!     UBj        J-dimension Upper bound.                              !
!     LBk        K-dimension Lower bound.                              !
!     UBk        K-dimension Upper bound.                              !
!     ifield     State field identification to process.                !
!     Mobs       Observation dimension in the calling program.         !
!     NobsSTR    Starting observation to process.                      !
!     NobsEND    Last observations to process.                         !
!     Xmin       Global minimum fractional I-coordinate to consider.   !
!     Xmax       Global maximum fractional I-coordinate to consider.   !
!     Ymin       Global minimum fractional J-coordinate to consider.   !
!     Ymax       Global maximum fractional J-coordinate to consider.   !
!     time       Current model time (secods).                          !
!     dt         Model baroclinic time-step (seconds).                 !
!     ObsType    Observations type.                                    !
!     Tobs       Observations time (days).                             !
!     Xobs       Observations X-locations (grid coordinates).          !
!     Yobs       Observations Y-locations (grid coordinates).          !
!     Zobs       Observations Z-locations (grid coordinates or meters).!
!     ad_Aobs    Adjoint forcing term at observations locations.       !
!     Adepth     Depths (meter; negative).                             !
!     Amask      Land-sea masking.                                     !
!                                                                      !
!  On Output:                                                          !
!                                                                      !
!     ObsVetting Observation screenning flag, 0: reject or 1: accept.  !
!     ad_A       Adjoint equations forcing array (2D or 3D).           !
!     ad_Aobs    Adjoint forcing term at observations locations.       !
!     Zobs       Observations Z-locations (grid coordinates).          !
!                                                                      !
!  The interpolation weights matrix, Hmat(1:8), is as follows:         !
!                                                                      !
!                               8____________7                         !
!                               /.          /| (i2,j2,k2)              !
!                              / .         / |                         !
!                            5/___________/6 |                         !
!                             |  .        |  |                         !
!                             |  .        |  |         Grid Cell       !
!                             | 4.........|..|3                        !
!                             | .         |  /                         !
!                             |.          | /                          !
!                  (i1,j1,k1) |___________|/                           !
!                             1           2                            !
!                                                                      !
!  Notice that the indices i2 and j2 are reset when observations are   !
!  located exactly at the eastern and/or northern boundaries. This is  !
!  needed to avoid out-of-range array computations.                    !
!                                                                      !
!  All the observations are assumed to in fractional coordinates with  !
!  respect to RHO-points:                                              !
!                                                                      !
!                                                                      !
!  M      r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r      !
!         :                                                     :      !
!  Mm+.5  v  p++v++p++v++p++v++p++v++p++v++p++v++p++v++p++v++p  v      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!  Mm     r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!  Mm-.5  v  p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p  v      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!         r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!         v  p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p  v      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!         r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!  2.5    v  p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p  v      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!  2.0    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!  1.5    v  p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p  v      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!  1.0    r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r  u  r      !
!         :  +     |     |     |     |     |     |     |     +  :      !
!  0.5    v  p++v++p++v++p++v++p++v++p++v++p++v++p++v++p++v++p  v      !
!         :                                                     :      !
!  0.0    r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r      !
!                                                                      !
!           0.5   1.5   2.5                          Lm-.5 Lm+.5       !
!                                                                      !
!        0.0   1.0   2.0                                  Lm    L      !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
 
      implicit none
 
      PUBLIC :: ad_extract_obs2d
# ifdef SOLVE3D
      PUBLIC :: ad_extract_obs3d
# endif
 
      CONTAINS
!
!***********************************************************************
      SUBROUTINE ad_extract_obs2d (ng, Imin, Imax, Jmin, Jmax,          &
     &                             LBi, UBi, LBj, UBj,                  &
     &                             ifield, Mobs, NobsSTR, NobsEND,      &
     &                             Xmin, Xmax, Ymin, Ymax,              &
     &                             time, dt,                            &
     &                             ObsType, ObsVetting,                 &
     &                             Tobs, Xobs, Yobs,                    &
     &                             ad_A,                                &
# ifdef MASKING
     &                             Amask,                               &
# endif
     &                             ad_Aobs)
!***********************************************************************
!
      USE mod_ncparam,  ONLY : isubar, isvbar
      USE mod_fourdvar, ONLY : obsstate2type
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, imin, imax, jmin, jmax
      integer, intent(in) :: lbi, ubi, lbj, ubj
      integer, intent(in) :: ifield, mobs, nobsstr, nobsend
!
      real(r8), intent(in) :: xmin, xmax, ymin, ymax
      real(dp), intent(in) :: time, dt
!
# ifdef ASSUMED_SHAPE
      integer, intent(in) :: obstype(:)
 
      real(dp), intent(in) :: tobs(:)
      real(r8), intent(in) :: xobs(:)
      real(r8), intent(in) :: yobs(:)
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
#  endif
      real(r8), intent(inout) :: obsvetting(:)
      real(r8), intent(inout) :: ad_aobs(:)
      real(r8), intent(inout) :: ad_a(lbi:,lbj:)
 
# else
      integer, intent(in) :: obstype(mobs)
 
      real(dp), intent(in) :: tobs(mobs)
      real(r8), intent(in) :: xobs(mobs)
      real(r8), intent(in) :: yobs(mobs)
#  ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
#  endif
      real(r8), intent(inout) :: obsvetting(mobs)
      real(r8), intent(inout) :: ad_aobs(mobs)
      real(r8), intent(inout) :: ad_a(lbi:ubi,lbj:ubj)
# endif
!
!  Local variable declarations.
!
      integer :: ic, iobs, i1, i2, j1, j2
 
      real(dp) :: timelb, timeub
 
      real(r8) ::  p1, p2, q1, q2, wsum
 
      real(r8), dimension(8) :: hmat
!
!-----------------------------------------------------------------------
!  Interpolate from requested 2D state field when appropriate.
!-----------------------------------------------------------------------
!
      timelb=(time-0.5_dp*dt)/86400.0_dp
      timeub=(time+0.5_dp*dt)/86400.0_dp
!
      DO iobs=nobsstr,nobsend
        IF ((obstype(iobs).eq.ifield).and.                              &
     &      ((timelb.le.tobs(iobs)).and.(tobs(iobs).lt.timeub)).and.    &
     &      ((xmin.le.xobs(iobs)).and.(xobs(iobs).lt.xmax)).and.        &
     &      ((ymin.le.yobs(iobs)).and.(yobs(iobs).lt.ymax))) THEN
          IF (obstype(iobs).eq.obsstate2type(isubar)) THEN
            i1=int(xobs(iobs)+0.5_r8)        ! 2D U-grid type variable
            j1=int(yobs(iobs))
          ELSE IF (obstype(iobs).eq.obsstate2type(isvbar)) THEN
            i1=int(xobs(iobs))               ! 2D V-grid type variable
            j1=int(yobs(iobs)+0.5_r8)
          ELSE
            i1=int(xobs(iobs))               ! 2D RHO-grid type variable
            j1=int(yobs(iobs))
          END IF
          i2=i1+1
          j2=j1+1
          IF (i2.gt.imax) THEN
            i2=i1                 ! Observation at the eastern boundary
          END IF
          IF (j2.gt.jmax) THEN
            j2=j1                 ! Observation at the northern boundary
          END IF
          p2=real(i2-i1,r8)*(xobs(iobs)-real(i1,r8))
          q2=real(j2-j1,r8)*(yobs(iobs)-real(j1,r8))
          p1=1.0_r8-p2
          q1=1.0_r8-q2
          hmat(1)=p1*q1
          hmat(2)=p2*q1
          hmat(3)=p2*q2
          hmat(4)=p1*q2
# ifdef MASKING
          hmat(1)=hmat(1)*amask(i1,j1)
          hmat(2)=hmat(2)*amask(i2,j1)
          hmat(3)=hmat(3)*amask(i2,j2)
          hmat(4)=hmat(4)*amask(i1,j2)
          wsum=0.0_r8
          DO ic=1,4
            wsum=wsum+hmat(ic)
          END DO
          IF (wsum.gt.0.0_r8) THEN
            wsum=1.0_r8/wsum
            DO ic=1,4
              hmat(ic)=hmat(ic)*wsum
            END DO
          END IF
# endif
!^        Aobs(iobs)=Hmat(1)*A(i1,j1)+                                  &
!^   &               Hmat(2)*A(i2,j1)+                                  &
!^   &               Hmat(3)*A(i2,j2)+                                  &
!^   &               Hmat(4)*A(i1,j2)
!^
          ad_a(i1,j1)=ad_a(i1,j1)+hmat(1)*ad_aobs(iobs)
          ad_a(i2,j1)=ad_a(i2,j1)+hmat(2)*ad_aobs(iobs)
          ad_a(i2,j2)=ad_a(i2,j2)+hmat(3)*ad_aobs(iobs)
          ad_a(i1,j2)=ad_a(i1,j2)+hmat(4)*ad_aobs(iobs)
          ad_aobs(iobs)=0.0_r8
# ifdef MASKING
          IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
# else
          obsvetting(iobs)=1.0_r8
# endif
        END IF
      END DO
 
      RETURN
      END SUBROUTINE ad_extract_obs2d
 
# ifdef SOLVE3D
!
!***********************************************************************
      SUBROUTINE ad_extract_obs3d (ng, Imin, Imax, Jmin, Jmax,          &
     &                             LBi, UBi, LBj, UBj, LBk, UBk,        &
     &                             ifield, Mobs, NobsSTR, NobsEND,      &
     &                             Xmin, Xmax, Ymin, Ymax,              &
     &                             time, dt,                            &
     &                             ObsType, ObsVetting,                 &
     &                             Tobs, Xobs, Yobs, Zobs,              &
     &                             ad_A, Adepth,                        &
#  ifdef MASKING
     &                             Amask,                               &
#  endif
     &                             ad_Aobs)
!***********************************************************************
!
      USE mod_param
!
      USE mod_ncparam,  ONLY : isuvel, isvvel
      USE mod_fourdvar, ONLY : obsstate2type
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, imin, imax, jmin, jmax
      integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk
      integer, intent(in) :: ifield, mobs, nobsstr, nobsend
!
      real(r8), intent(in) :: xmin, xmax, ymin, ymax
      real(dp), intent(in) :: time, dt
!
#  ifdef ASSUMED_SHAPE
      integer, intent(in) :: obstype(:)
 
      real(dp), intent(in) :: tobs(:)
      real(r8), intent(in) :: xobs(:)
      real(r8), intent(in) :: yobs(:)
      real(r8), intent(in) :: adepth(lbi:,lbj:,lbk:)
#   ifdef MASKING
      real(r8), intent(in) :: amask(lbi:,lbj:)
#   endif
      real(r8), intent(inout) :: obsvetting(:)
      real(r8), intent(inout) :: ad_aobs(:)
      real(r8), intent(inout) :: zobs(:)
      real(r8), intent(inout) :: ad_a(lbi:,lbj:,lbk:)
#  else
      integer, intent(in) :: obstype(mobs)
 
      real(dp), intent(in) :: tobs(mobs)
      real(r8), intent(in) :: xobs(mobs)
      real(r8), intent(in) :: yobs(mobs)
      real(r8), intent(in) :: adepth(lbi:ubi,lbj:ubj,lbk:ubk)
#   ifdef MASKING
      real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
#   endif
      real(r8), intent(inout) :: obsvetting(mobs)
      real(r8), intent(inout) :: ad_aobs(mobs)
      real(r8), intent(inout) :: zobs(mobs)
      real(r8), intent(inout) :: ad_a(lbi:ubi,lbj:ubj,lbk:ubk)
#  endif
!
!  Local variable declarations.
!
      integer :: i, ic, iobs, i1, i2, j1, j2, k, k1, k2
 
      real(dp) :: timelb, timeub
 
      real(r8) :: zbot, ztop, dz, p1, p2, q1, q2, r1, r2
      real(r8) :: w11, w12, w21, w22, wsum
 
      real(r8), dimension(8) :: hmat
!
!-----------------------------------------------------------------------
!  Interpolate from requested 3D state field.
!-----------------------------------------------------------------------
!
      timelb=(time-0.5_dp*dt)/86400.0_dp
      timeub=(time+0.5_dp*dt)/86400.0_dp
!
      DO iobs=nobsstr,nobsend
        IF ((obstype(iobs).eq.ifield).and.                              &
     &      ((timelb.le.tobs(iobs)).and.(tobs(iobs).lt.timeub)).and.    &
     &      ((xmin.le.xobs(iobs)).and.(xobs(iobs).lt.xmax)).and.        &
     &      ((ymin.le.yobs(iobs)).and.(yobs(iobs).lt.ymax))) THEN
          IF (obstype(iobs).eq.obsstate2type(isuvel)) THEN
            i1=int(xobs(iobs)+0.5_r8)        ! 3D U-grid type variable
            j1=int(yobs(iobs))
          ELSE IF (obstype(iobs).eq.obsstate2type(isvvel)) THEN
            i1=int(xobs(iobs))               ! 3D V-grid type variable
            j1=int(yobs(iobs)+0.5_r8)
          ELSE
            i1=int(xobs(iobs))               ! 3D RHO-grid type variable
            j1=int(yobs(iobs))
          END IF
          i2=i1+1
          j2=j1+1
          IF (i2.gt.imax) THEN
            i2=i1                 ! Observation at the eastern boundary
          END IF
          IF (j2.gt.jmax) THEN
            j2=j1                 ! Observation at the northern boundary
          END IF
          p2=real(i2-i1,r8)*(xobs(iobs)-real(i1,r8))
          q2=real(j2-j1,r8)*(yobs(iobs)-real(j1,r8))
          p1=1.0_r8-p2
          q1=1.0_r8-q2
          w11=p1*q1
          w21=p2*q1
          w22=p2*q2
          w12=p1*q2
          IF (zobs(iobs).gt.0.0_r8) THEN
            k1=max(1,int(zobs(iobs)))          ! Positions in fractional
            k2=min(k1+1,n(ng))                 ! levels
            r2=real(k2-k1,r8)*(zobs(iobs)-real(k1,r8))
            r1=1.0_r8-r2
          ELSE
            ztop=adepth(i1,j1,n(ng))
            zbot=adepth(i1,j1,1    )
            IF (zobs(iobs).ge.ztop) THEN
              k1=n(ng)                         ! If shallower, assign to
              k2=n(ng)                         ! top grid cell. The
              r1=1.0_r8                        ! observation is located
              r2=0.0_r8                        ! on the upper cell half
              zobs(iobs)=real(n(ng),r8)        ! above its middle depth.
            ELSE IF (zbot.ge.zobs(iobs)) THEN
              r1=0.0_r8                        ! If deeper, ignore.
              r2=0.0_r8
              obsvetting(iobs)=0.0_r8
            ELSE
              DO k=n(ng),2,-1                  ! Otherwise, interpolate
                ztop=adepth(i1,j1,k  )         ! to fractional level
                zbot=adepth(i1,j1,k-1)
                IF ((ztop.gt.zobs(iobs)).and.(zobs(iobs).ge.zbot)) THEN
                  k1=k-1
                  k2=k
                END IF
              END DO
              dz=adepth(i1,j1,k2)-adepth(i1,j1,k1)
              r2=(zobs(iobs)-adepth(i1,j1,k1))/dz
              r1=1.0_r8-r2
              zobs(iobs)=real(k1,r8)+r2        ! overwrite
            END IF
          END IF
          IF ((r1+r2).gt.0.0_r8) THEN
            hmat(1)=w11*r1
            hmat(2)=w21*r1
            hmat(3)=w22*r1
            hmat(4)=w12*r1
            hmat(5)=w11*r2
            hmat(6)=w21*r2
            hmat(7)=w22*r2
            hmat(8)=w12*r2
#  ifdef MASKING
            hmat(1)=hmat(1)*amask(i1,j1)
            hmat(2)=hmat(2)*amask(i2,j1)
            hmat(3)=hmat(3)*amask(i2,j2)
            hmat(4)=hmat(4)*amask(i1,j2)
            hmat(5)=hmat(5)*amask(i1,j1)
            hmat(6)=hmat(6)*amask(i2,j1)
            hmat(7)=hmat(7)*amask(i2,j2)
            hmat(8)=hmat(8)*amask(i1,j2)
            wsum=0.0_r8
            DO ic=1,8
              wsum=wsum+hmat(ic)
            END DO
            IF (wsum.gt.0.0_r8) THEN
              wsum=1.0_r8/wsum
              DO ic=1,8
                hmat(ic)=hmat(ic)*wsum
              END DO
            END IF
#  endif
!^          Aobs(iobs)=Hmat(1)*A(i1,j1,k1)+                             &
!^   &                 Hmat(2)*A(i2,j1,k1)+                             &
!^   &                 Hmat(3)*A(i2,j2,k1)+                             &
!^   &                 Hmat(4)*A(i1,j2,k1)+                             &
!^   &                 Hmat(5)*A(i1,j1,k2)+                             &
!^   &                 Hmat(6)*A(i2,j1,k2)+                             &
!^   &                 Hmat(7)*A(i2,j2,k2)+                             &
!^   &                 Hmat(8)*A(i1,j2,k2)
!^
            ad_a(i1,j1,k1)=ad_a(i1,j1,k1)+hmat(1)*ad_aobs(iobs)
            ad_a(i2,j1,k1)=ad_a(i2,j1,k1)+hmat(2)*ad_aobs(iobs)
            ad_a(i2,j2,k1)=ad_a(i2,j2,k1)+hmat(3)*ad_aobs(iobs)
            ad_a(i1,j2,k1)=ad_a(i1,j2,k1)+hmat(4)*ad_aobs(iobs)
            ad_a(i1,j1,k2)=ad_a(i1,j1,k2)+hmat(5)*ad_aobs(iobs)
            ad_a(i2,j1,k2)=ad_a(i2,j1,k2)+hmat(6)*ad_aobs(iobs)
            ad_a(i2,j2,k2)=ad_a(i2,j2,k2)+hmat(7)*ad_aobs(iobs)
            ad_a(i1,j2,k2)=ad_a(i1,j2,k2)+hmat(8)*ad_aobs(iobs)
            ad_aobs(iobs)=0.0_r8
#  ifdef MASKING
            IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
#  else
            obsvetting(iobs)=1.0_r8
#  endif
          END IF
        END IF
      END DO
 
      RETURN
      END SUBROUTINE ad_extract_obs3d
# endif
#endif
      END MODULE ad_extract_obs_mod