obs_cost.F

Go to the documentation of this file.

#include "cppdefs.h"
#if defined FOUR_DVAR && defined OBSERVATIONS
      SUBROUTINE obs_cost (ng, model)
!
!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                                               !
!=======================================================================
 
# ifdef WEAK_CONSTRAINT
#  if defined R4DVAR    || defined R4DVAR_ANA_SENSITIVITY || \
      defined tl_r4dvar
!                                                                      !
!  This routine computes the data penalty function directly in during  !
!  runs of the representer model:                                      !
!                                                                      !
#  else
!                                                                      !
!  This routine computes the data penalty function directly in during  !
!  runs of the nonlinear model:                                        !
!                                                                      !
#  endif
!         Jdata = transpose(H X - Xo) * O^(-1) * (H X - Xo)            !
!                                                                      !
!         H  : observation operator (linearized if incremental)        !
!         Xo : observations vector                                     !
!       H X  : representer model at observation points                 !
!         O  : observations error covariance                           !
# else
!                                                                      !
!  This routine computes the observation cost function (Jo) as the     !
!  misfit (squared difference) between the model and observations.     !
!                                                                      !
!  If conventional strong contraint 4D-Var:                            !
!                                                                      !
!         Jo = 1/2 transpose(H X - Xo) * O^(-1) * (H X - Xo)           !
!                                                                      !
!  or if incremental strong contraint 4D-Var:                          !
!                                                                      !
!         Jo = 1/2 transpose(H deltaX - d) * O^(-1) * (H deltaX - d)   !
!                                                                      !
!  where                                                               !
!                                                                      !
!          d = Xo - H Xb                                               !
!                                                                      !
!         d  : innovation vector                                       !
!         H  : observation operator (linearized if incremental)        !
!       H Xb : background at observation points previous forecast)     !
!         Xo : observations vector                                     !
!       H X  : nonlinear model at observation points                   !
!  H deltaX  : increment at observation point                          !
!         O  : observations error covariance                           !
# endif
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_scalars
!
      implicit none
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, model
!
!  Local variable declarations.
!
      integer :: NSUB, iobs, ivar
 
      real(r8) ::  cff, cff1
 
      real(r8), dimension(0:NobsVar(ng)) :: my_ObsCost
# if defined DATALESS_LOOPS && \
    (defined r4dvar         || defined r4dvar_ana_sensitivity || \
     defined tl_r4dvar)
      real(r8), dimension(0:NobsVar(ng)) :: my_ObsCost1
# endif
!
!-----------------------------------------------------------------------
!  Compute observation misfit cost function (ObsCost).
!-----------------------------------------------------------------------
 
# if defined R4DVAR    || defined R4DVAR_ANA_SENSITIVITY || \
     defined tl_r4dvar
!
!  Compute data penalty function.
!
      IF (model.eq.irpm) THEN
        DO ivar=0,nobsvar(ng)
          my_obscost(ivar)=0.0_r8
#  ifdef DATALESS_LOOPS
          my_obscost1(ivar)=0.0_r8
#  endif
        END DO
        DO iobs=nstrobs(ng),nendobs(ng)
          ivar=obstype2state(obstype(iobs))
          IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8).and.           &
     &        (obserr(iobs).ne.0.0_r8)) THEN
            cff=obsscale(iobs)*(tlmodval(iobs)-obsval(iobs))**2/        &
     &          obserr(iobs)
#  ifdef DATALESS_LOOPS
            cff1=obsscale(iobs)*(nlmodval(iobs)-obsval(iobs))**2/       &
     &           obserr(iobs)
#  endif
            my_obscost(0)=my_obscost(0)+cff
            my_obscost(ivar)=my_obscost(ivar)+cff
#  ifdef DATALESS_LOOPS
            my_obscost1(0)=my_obscost1(0)+cff1
            my_obscost1(ivar)=my_obscost1(ivar)+cff1
#  endif
          END IF
        END DO
      END IF
 
# elif defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY || \
       defined sp4dvar  || defined tl_rbl4dvar
!
!  Compute nonlinear model data penalty function.
!
      IF (model.eq.inlm) THEN
        DO ivar=0,nobsvar(ng)
          my_obscost(ivar)=0.0_r8
        END DO
        DO iobs=nstrobs(ng),nendobs(ng)
          ivar=obstype2state(obstype(iobs))
          IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8).and.           &
     &        (obserr(iobs).ne.0.0_r8)) THEN
            cff=obsscale(iobs)*(nlmodval(iobs)-obsval(iobs))**2/        &
     &          obserr(iobs)
            my_obscost(0)=my_obscost(0)+cff
            my_obscost(ivar)=my_obscost(ivar)+cff
          END IF
        END DO
      END IF
# else
!
!  Compute tangent linear model cost function.
!
      IF (model.eq.itlm) THEN
        DO ivar=0,nobsvar(ng)
          my_obscost(ivar)=0.0_r8
        END DO
        DO iobs=1,nobs(ng)
          ivar=obstype2state(obstype(iobs))
          IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8)) THEN
            cff=0.5_r8*obsscale(iobs)*obserr(iobs)*                     &
     &          (nlmodval(iobs)+tlmodval(iobs)-obsval(iobs))**2
            my_obscost(0)=my_obscost(0)+cff
            my_obscost(ivar)=my_obscost(ivar)+cff
          END IF
        END DO
!
!  Compute nonlinear model cost function.
!
      ELSE IF (model.eq.inlm) THEN
        DO ivar=0,nobsvar(ng)
          my_obscost(ivar)=0.0_r8
        END DO
        DO iobs=1,nobs(ng)
          ivar=obstype2state(obstype(iobs))
          IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8)) THEN
            cff=0.5_r8*obsscale(iobs)*obserr(iobs)*                     &
     &          (nlmodval(iobs)-obsval(iobs))**2
            my_obscost(0)=my_obscost(0)+cff
            my_obscost(ivar)=my_obscost(ivar)+cff
          END IF
        END DO
      END IF
# endif
!
!-----------------------------------------------------------------------
!  Load global values.  Notice that there is not need for a global
!  reduction here since all the threads have the same copy of all
!  the vectors used.
!-----------------------------------------------------------------------
!
# ifdef DISTRIBUTE
      nsub=1                             ! distributed-memory
# else
      nsub=ntilex(ng)*ntilee(ng)
# endif
!$OMP CRITICAL (COST_FUN)
      tile_count=tile_count+1
      IF (tile_count.eq.nsub) THEN
        tile_count=0
# if defined R4DVAR    || defined R4DVAR_ANA_SENSITIVITY || \
     defined tl_r4dvar
        IF (model.eq.irpm) THEN
          DO ivar=0,nobsvar(ng)
            fourdvar(ng)%DataPenalty(ivar)=my_obscost(ivar)+            &
     &                                    fourdvar(ng)%DataPenalty(ivar)
#  ifdef DATALESS_LOOPS
            fourdvar(ng)%NLPenalty(ivar)=my_obscost1(ivar)+             &
     &                                   fourdvar(ng)%NLPenalty(ivar)
#  endif
          END DO
        END IF
# elif defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY || \
       defined sp4dvar  || defined tl_rbl4dvar
        IF (model.eq.inlm) THEN
          DO ivar=0,nobsvar(ng)
            fourdvar(ng)%NLPenalty(ivar)=my_obscost(ivar)+              &
     &                                   fourdvar(ng)%NLPenalty(ivar)
          END DO
        END IF
# else
        IF (model.eq.itlm) THEN
          DO ivar=0,nobsvar(ng)
            fourdvar(ng)%ObsCost(ivar)=fourdvar(ng)%ObsCost(ivar)+      &
     &                                 my_obscost(ivar)
          END DO
        ELSE IF (model.eq.inlm) THEN
          DO ivar=0,nobsvar(ng)
            fourdvar(ng)%NLobsCost(ivar)=fourdvar(ng)%NLobsCost(ivar)+  &
     &                                   my_obscost(ivar)
          END DO
        END IF
# endif
      END IF
!$OMP END CRITICAL (COST_FUN)
 
# ifndef WEAK_CONSTRAINT
!
!  If start of minimization, set cost function scales used to report
!  normalized values.
!
      IF ((nrun.eq.1).and.(model.eq.itlm)) THEN
        DO ivar=0,nobsvar(ng)
          fourdvar(ng)%CostNorm(ivar)=fourdvar(ng)%ObsCost(ivar)
        END DO
      END IF
!
!  Save initial inner loop cost function.
!
      IF ((inner.eq.0).and.(model.eq.itlm)) THEN
        fourdvar(ng)%Cost0(outer)=fourdvar(ng)%ObsCost(0)
      END IF
# endif
 
      RETURN
      END SUBROUTINE obs_cost
#else
      SUBROUTINE obs_cost
      RETURN
      END SUBROUTINE obs_cost
#endif