ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
obs_cost.F File Reference
#include "cppdefs.h"
Include dependency graph for obs_cost.F:

Go to the source code of this file.

Functions/Subroutines

subroutine obs_cost (ng, model)
 

Function/Subroutine Documentation

◆ obs_cost()

subroutine obs_cost ( integer, intent(in) ng,
integer, intent(in) model )

Definition at line 3 of file obs_cost.F.

4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11
12# ifdef WEAK_CONSTRAINT
13# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
14 defined tl_r4dvar
15! !
16! This routine computes the data penalty function directly in during !
17! runs of the representer model: !
18! !
19# else
20! !
21! This routine computes the data penalty function directly in during !
22! runs of the nonlinear model: !
23! !
24# endif
25! Jdata = transpose(H X - Xo) * O^(-1) * (H X - Xo) !
26! !
27! H : observation operator (linearized if incremental) !
28! Xo : observations vector !
29! H X : representer model at observation points !
30! O : observations error covariance !
31# else
32! !
33! This routine computes the observation cost function (Jo) as the !
34! misfit (squared difference) between the model and observations. !
35! !
36! If conventional strong contraint 4D-Var: !
37! !
38! Jo = 1/2 transpose(H X - Xo) * O^(-1) * (H X - Xo) !
39! !
40! or if incremental strong contraint 4D-Var: !
41! !
42! Jo = 1/2 transpose(H deltaX - d) * O^(-1) * (H deltaX - d) !
43! !
44! where !
45! !
46! d = Xo - H Xb !
47! !
48! d : innovation vector !
49! H : observation operator (linearized if incremental) !
50! H Xb : background at observation points previous forecast) !
51! Xo : observations vector !
52! H X : nonlinear model at observation points !
53! H deltaX : increment at observation point !
54! O : observations error covariance !
55# endif
56! !
57!=======================================================================
58!
59 USE mod_param
60 USE mod_parallel
61 USE mod_fourdvar
62 USE mod_scalars
63!
64 implicit none
65!
66! Imported variable declarations.
67!
68 integer, intent(in) :: ng, model
69!
70! Local variable declarations.
71!
72 integer :: NSUB, iobs, ivar
73
74 real(r8) :: cff, cff1
75
76 real(r8), dimension(0:NobsVar(ng)) :: my_ObsCost
77# if defined DATALESS_LOOPS && \
78 (defined r4dvar || defined r4dvar_ana_sensitivity || \
79 defined tl_r4dvar)
80 real(r8), dimension(0:NobsVar(ng)) :: my_ObsCost1
81# endif
82!
83!-----------------------------------------------------------------------
84! Compute observation misfit cost function (ObsCost).
85!-----------------------------------------------------------------------
86
87# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
88 defined tl_r4dvar
89!
90! Compute data penalty function.
91!
92 IF (model.eq.irpm) THEN
93 DO ivar=0,nobsvar(ng)
94 my_obscost(ivar)=0.0_r8
95# ifdef DATALESS_LOOPS
96 my_obscost1(ivar)=0.0_r8
97# endif
98 END DO
99 DO iobs=nstrobs(ng),nendobs(ng)
100 ivar=obstype2state(obstype(iobs))
101 IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8).and. &
102 & (obserr(iobs).ne.0.0_r8)) THEN
103 cff=obsscale(iobs)*(tlmodval(iobs)-obsval(iobs))**2/ &
104 & obserr(iobs)
105# ifdef DATALESS_LOOPS
106 cff1=obsscale(iobs)*(nlmodval(iobs)-obsval(iobs))**2/ &
107 & obserr(iobs)
108# endif
109 my_obscost(0)=my_obscost(0)+cff
110 my_obscost(ivar)=my_obscost(ivar)+cff
111# ifdef DATALESS_LOOPS
112 my_obscost1(0)=my_obscost1(0)+cff1
113 my_obscost1(ivar)=my_obscost1(ivar)+cff1
114# endif
115 END IF
116 END DO
117 END IF
118
119# elif defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY || \
120 defined sp4dvar || defined tl_rbl4dvar
121!
122! Compute nonlinear model data penalty function.
123!
124 IF (model.eq.inlm) THEN
125 DO ivar=0,nobsvar(ng)
126 my_obscost(ivar)=0.0_r8
127 END DO
128 DO iobs=nstrobs(ng),nendobs(ng)
129 ivar=obstype2state(obstype(iobs))
130 IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8).and. &
131 & (obserr(iobs).ne.0.0_r8)) THEN
132 cff=obsscale(iobs)*(nlmodval(iobs)-obsval(iobs))**2/ &
133 & obserr(iobs)
134 my_obscost(0)=my_obscost(0)+cff
135 my_obscost(ivar)=my_obscost(ivar)+cff
136 END IF
137 END DO
138 END IF
139# else
140!
141! Compute tangent linear model cost function.
142!
143 IF (model.eq.itlm) THEN
144 DO ivar=0,nobsvar(ng)
145 my_obscost(ivar)=0.0_r8
146 END DO
147 DO iobs=1,nobs(ng)
148 ivar=obstype2state(obstype(iobs))
149 IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8)) THEN
150 cff=0.5_r8*obsscale(iobs)*obserr(iobs)* &
151 & (nlmodval(iobs)+tlmodval(iobs)-obsval(iobs))**2
152 my_obscost(0)=my_obscost(0)+cff
153 my_obscost(ivar)=my_obscost(ivar)+cff
154 END IF
155 END DO
156!
157! Compute nonlinear model cost function.
158!
159 ELSE IF (model.eq.inlm) THEN
160 DO ivar=0,nobsvar(ng)
161 my_obscost(ivar)=0.0_r8
162 END DO
163 DO iobs=1,nobs(ng)
164 ivar=obstype2state(obstype(iobs))
165 IF ((ivar.gt.0).and.(obsscale(iobs).gt.0.0_r8)) THEN
166 cff=0.5_r8*obsscale(iobs)*obserr(iobs)* &
167 & (nlmodval(iobs)-obsval(iobs))**2
168 my_obscost(0)=my_obscost(0)+cff
169 my_obscost(ivar)=my_obscost(ivar)+cff
170 END IF
171 END DO
172 END IF
173# endif
174!
175!-----------------------------------------------------------------------
176! Load global values. Notice that there is not need for a global
177! reduction here since all the threads have the same copy of all
178! the vectors used.
179!-----------------------------------------------------------------------
180!
181# ifdef DISTRIBUTE
182 nsub=1 ! distributed-memory
183# else
184 nsub=ntilex(ng)*ntilee(ng)
185# endif
186!$OMP CRITICAL (COST_FUN)
187 tile_count=tile_count+1
188 IF (tile_count.eq.nsub) THEN
189 tile_count=0
190# if defined R4DVAR || defined R4DVAR_ANA_SENSITIVITY || \
191 defined tl_r4dvar
192 IF (model.eq.irpm) THEN
193 DO ivar=0,nobsvar(ng)
194 fourdvar(ng)%DataPenalty(ivar)=my_obscost(ivar)+ &
195 & fourdvar(ng)%DataPenalty(ivar)
196# ifdef DATALESS_LOOPS
197 fourdvar(ng)%NLPenalty(ivar)=my_obscost1(ivar)+ &
198 & fourdvar(ng)%NLPenalty(ivar)
199# endif
200 END DO
201 END IF
202# elif defined RBL4DVAR || defined RBL4DVAR_ANA_SENSITIVITY || \
203 defined sp4dvar || defined tl_rbl4dvar
204 IF (model.eq.inlm) THEN
205 DO ivar=0,nobsvar(ng)
206 fourdvar(ng)%NLPenalty(ivar)=my_obscost(ivar)+ &
207 & fourdvar(ng)%NLPenalty(ivar)
208 END DO
209 END IF
210# else
211 IF (model.eq.itlm) THEN
212 DO ivar=0,nobsvar(ng)
213 fourdvar(ng)%ObsCost(ivar)=fourdvar(ng)%ObsCost(ivar)+ &
214 & my_obscost(ivar)
215 END DO
216 ELSE IF (model.eq.inlm) THEN
217 DO ivar=0,nobsvar(ng)
218 fourdvar(ng)%NLobsCost(ivar)=fourdvar(ng)%NLobsCost(ivar)+ &
219 & my_obscost(ivar)
220 END DO
221 END IF
222# endif
223 END IF
224!$OMP END CRITICAL (COST_FUN)
225
226# ifndef WEAK_CONSTRAINT
227!
228! If start of minimization, set cost function scales used to report
229! normalized values.
230!
231 IF ((nrun.eq.1).and.(model.eq.itlm)) THEN
232 DO ivar=0,nobsvar(ng)
233 fourdvar(ng)%CostNorm(ivar)=fourdvar(ng)%ObsCost(ivar)
234 END DO
235 END IF
236!
237! Save initial inner loop cost function.
238!
239 IF ((inner.eq.0).and.(model.eq.itlm)) THEN
240 fourdvar(ng)%Cost0(outer)=fourdvar(ng)%ObsCost(0)
241 END IF
242# endif
243
244 RETURN
mod_fourdvar
Definition mod_fourdvar.F:2
mod_fourdvar::fourdvar
type(t_fourdvar), dimension(:), allocatable fourdvar
Definition mod_fourdvar.F:211
mod_fourdvar::obstype2state
integer, dimension(:), allocatable obstype2state
Definition mod_fourdvar.F:219
mod_fourdvar::nobs
integer, dimension(:), allocatable nobs
Definition mod_fourdvar.F:421
mod_fourdvar::nobsvar
integer, dimension(:), allocatable nobsvar
Definition mod_fourdvar.F:372
mod_fourdvar::obsval
real(r8), dimension(:), allocatable obsval
Definition mod_fourdvar.F:232
mod_fourdvar::obsscale
real(r8), dimension(:), allocatable obsscale
Definition mod_fourdvar.F:227
mod_fourdvar::obserr
real(r8), dimension(:), allocatable obserr
Definition mod_fourdvar.F:225
mod_fourdvar::obstype
integer, dimension(:), allocatable obstype
Definition mod_fourdvar.F:217
mod_fourdvar::nlmodval
real(r8), dimension(:), allocatable nlmodval
Definition mod_fourdvar.F:240
mod_fourdvar::nstrobs
integer, dimension(:), allocatable nstrobs
Definition mod_fourdvar.F:425
mod_fourdvar::nendobs
integer, dimension(:), allocatable nendobs
Definition mod_fourdvar.F:426
mod_parallel
Definition mod_parallel.F:2
mod_parallel::tile_count
integer tile_count
Definition mod_parallel.F:89
mod_param
Definition mod_param.F:2
mod_param::inlm
integer, parameter inlm
Definition mod_param.F:662
mod_param::irpm
integer, parameter irpm
Definition mod_param.F:664
mod_param::ntilex
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
mod_param::ntilee
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
mod_param::itlm
integer, parameter itlm
Definition mod_param.F:663
mod_scalars
Definition mod_scalars.F:2
mod_scalars::nrun
integer nrun
Definition mod_scalars.F:355
mod_scalars::inner
integer inner
Definition mod_scalars.F:360
mod_scalars::outer
integer outer
Definition mod_scalars.F:361

References mod_fourdvar::fourdvar, mod_param::inlm, mod_scalars::inner, mod_param::irpm, mod_param::itlm, mod_fourdvar::nendobs, mod_fourdvar::nlmodval, mod_fourdvar::nobs, mod_fourdvar::nobsvar, mod_scalars::nrun, mod_fourdvar::nstrobs, mod_param::ntilee, mod_param::ntilex, mod_fourdvar::obserr, mod_fourdvar::obsscale, mod_fourdvar::obstype, mod_fourdvar::obstype2state, mod_fourdvar::obsval, mod_scalars::outer, and mod_parallel::tile_count.

Referenced by output(), rp_output(), and tl_output().

Here is the caller graph for this function: