#include "cppdefs.h"
#if defined FOUR_DVAR && defined OBSERVATIONS
      SUBROUTINE obs_cost (ng, model)
!
!svn $Id: obs_cost.F 889 2018-02-10 03:32:52Z arango $
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2019 The ROMS/TOMS Group       Andrew M. Moore   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.txt                                              !
!=======================================================================

# ifdef WEAK_CONSTRAINT
#  if defined TL_W4DVAR || defined W4DVAR || defined W4DVAR_SENSITIVITY
!                                                                      !
!  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 4DVAR:                             !
!                                                                      !
!         Jo = 1/2 transpose(H X - Xo) * O^(-1) * (H X - Xo)           !
!                                                                      !
!  or if incremental strong contraint 4DVAR:                           !
!                                                                      !
!         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 TL_W4DVAR || defined W4DVAR || defined W4DVAR_SENSITIVITY)
      real(r8), dimension(0:NobsVar(ng)) :: my_ObsCost1
# endif
!
!-----------------------------------------------------------------------
!  Compute observation misfit cost function (ObsCost).
!-----------------------------------------------------------------------

# if defined TL_W4DVAR || defined W4DVAR || defined W4DVAR_SENSITIVITY
!
!  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 TL_W4DPSAS          || defined W4DPSAS || \
       defined W4DPSAS_SENSITIVITY
!
!  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 TL_W4DVAR || defined W4DVAR || defined W4DVAR_SENSITIVITY
        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 TL_W4DPSAS          || defined W4DPSAS || \
       defined W4DPSAS_SENSITIVITY
        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