ROMS
Loading...
Searching...
No Matches
extract_obs.F
Go to the documentation of this file.
1#include "cppdefs.h"
3#if (defined FOUR_DVAR || defined VERIFICATION) && defined OBSERVATIONS
4!
5!git $Id$
6!================================================== Hernan G. Arango ===
7! Copyright (c) 2002-2025 The ROMS Group !
8! Licensed under a MIT/X style license !
9! See License_ROMS.md !
10!=======================================================================
11! !
12! This routine extracts model data at the requested observations !
13! positions (Xobs,Yobs,Zobs). The extraction is done via linear !
14! interpolation. The (Xobs,Yobs) positions must be in fractional !
15! grid coordinates. Zobs can be in fractional grid coordinates !
16! (Zobs >= 0) or actual depths (Zobs < 0), if applicable. !
17! !
18! On Input: !
19! !
20! ng Nested grid number. !
21! Imin Global I-coordinate lower bound threshold for !
22! requested state field type. !
23! Imax Global I-coordinate upper bound threshold for !
24! requested state field type. !
25! Jmin Global J-coordinate lower bound threshold for !
26! requested state field type. !
27! Jmax Global J-coordinate upper bound threshold for !
28! requested state field type. !
29! LBi I-dimension Lower bound. !
30! UBi I-dimension Upper bound. !
31! LBj J-dimension Lower bound. !
32! UBj J-dimension Upper bound. !
33! LBk K-dimension Lower bound. !
34! UBk K-dimension Upper bound. !
35! ifield State field identification to process. !
36! Mobs Observation dimension in the calling program. !
37! NobsSTR Starting observation to process. !
38! NobsEND Last observations to process. !
39! Xmin Global minimum fractional I-coordinate to consider. !
40! Xmax Global maximum fractional I-coordinate to consider. !
41! Ymin Global minimum fractional J-coordinate to consider. !
42! Ymax Global maximum fractional J-coordinate to consider. !
43! time Current model time (seconds). !
44! dt Model baroclinic time-step (seconds). !
45! ObsType Observations type. !
46! Tobs Observations time (days). !
47! Xobs Observations X-locations (grid coordinates). !
48! Yobs Observations Y-locations (grid coordinates). !
49! Zobs Observations Z-locations (grid coordinates or meters).!
50! A Model array (2D or 3D) to process. !
51! Adepth Depths (meter; negative). !
52! Amask Land-sea masking. !
53! !
54! On Output: !
55! !
56! ObsVetting Observation screenning flag, 0: reject or 1: accept. !
57! Aobs Extracted model values at observation positions. !
58! Zobs Observations Z-locations (grid coordinates). !
59! !
60! The interpolation weights matrix, Hmat(1:8), is as follows: !
61! !
62! 8____________7 !
63! /. /| (i2,j2,k2) !
64! / . / | !
65! 5/___________/6 | !
66! | . | | !
67! | . | | Grid Cell !
68! | 4.........|..|3 !
69! | . | / !
70! |. | / !
71! (i1,j1,k1) |___________|/ !
72! 1 2 !
73! !
74! Notice that the indices i2 and j2 are reset when observations are !
75! located exactly at the eastern and/or northern boundaries. This is !
76! needed to avoid out-of-range array computations. !
77! !
78! All the observations are assumed to in fractional coordinates with !
79! respect to RHO-points: !
80! !
81! !
82! M r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r !
83! : : !
84! Mm+.5 v p++v++p++v++p++v++p++v++p++v++p++v++p++v++p++v++p v !
85! : + | | | | | | | + : !
86! Mm r u r u r u r u r u r u r u r u r u r !
87! : + | | | | | | | + : !
88! Mm-.5 v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v !
89! : + | | | | | | | + : !
90! r u r u r u r u r u r u r u r u r u r !
91! : + | | | | | | | + : !
92! v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v !
93! : + | | | | | | | + : !
94! r u r u r u r u r u r u r u r u r u r !
95! : + | | | | | | | + : !
96! 2.5 v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v !
97! : + | | | | | | | + : !
98! 2.0 r u r u r u r u r u r u r u r u r u r !
99! : + | | | | | | | + : !
100! 1.5 v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v !
101! : + | | | | | | | + : !
102! 1.0 r u r u r u r u r u r u r u r u r u r !
103! : + | | | | | | | + : !
104! 0.5 v p++v++p++v++p++v++p++v++p++v++p++v++p++v++p++v++p v !
105! : : !
106! 0.0 r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r !
107! !
108! 0.5 1.5 2.5 Lm-.5 Lm+.5 !
109! !
110! 0.0 1.0 2.0 Lm L !
111! !
112!=======================================================================
113!
114 USE mod_kinds
115
116 implicit none
117
118 PUBLIC extract_obs2d
119# ifdef SOLVE3D
120 PUBLIC extract_obs3d
121# endif
122
123 CONTAINS
124!
125!***********************************************************************
126 SUBROUTINE extract_obs2d (ng, Imin, Imax, Jmin, Jmax, &
127 & LBi, UBi, LBj, UBj, &
128 & ifield, Mobs, NobsSTR, NobsEND, &
129 & Xmin, Xmax, Ymin, Ymax, &
130 & time, dt, &
131 & ObsType, ObsVetting, &
132 & Tobs, Xobs, Yobs, &
133 & A, &
134# ifdef MASKING
135 & Amask, &
136# endif
137 & Aobs)
138!***********************************************************************
139!
140 USE mod_ncparam, ONLY : isubar, isvbar
141 USE mod_fourdvar, ONLY : obsstate2type
142!
143! Imported variable declarations.
144!
145 integer, intent(in) :: ng, imin, imax, jmin, jmax
146 integer, intent(in) :: lbi, ubi, lbj, ubj
147 integer, intent(in) :: ifield, mobs, nobsstr, nobsend
148!
149 real(r8), intent(in) :: xmin, xmax, ymin, ymax
150 real(dp), intent(in) :: time, dt
151!
152# ifdef ASSUMED_SHAPE
153 integer, intent(in) :: obstype(:)
154
155 real(dp), intent(in) :: tobs(:)
156 real(r8), intent(in) :: xobs(:)
157 real(r8), intent(in) :: yobs(:)
158 real(r8), intent(in) :: a(lbi:,lbj:)
159# ifdef MASKING
160 real(r8), intent(in) :: amask(lbi:,lbj:)
161# endif
162
163 real(r8), intent(inout) :: obsvetting(:)
164 real(r8), intent(inout) :: aobs(:)
165# else
166 integer, intent(in) :: obstype(mobs)
167
168 real(dp), intent(in) :: tobs(mobs)
169 real(r8), intent(in) :: xobs(mobs)
170 real(r8), intent(in) :: yobs(mobs)
171 real(r8), intent(in) :: a(lbi:ubi,lbj:ubj)
172# ifdef MASKING
173 real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
174# endif
175 real(r8), intent(inout) :: obsvetting(mobs)
176 real(r8), intent(inout) :: aobs(mobs)
177# endif
178!
179! Local variable declarations.
180!
181 integer :: ic, iobs, i1, i2, j1, j2
182
183 real(dp) :: timelb, timeub
184
185 real(r8) :: p1, p2, q1, q2, wsum
186
187 real(r8), dimension(8) :: hmat
188!
189!-----------------------------------------------------------------------
190! Interpolate from requested 2D state field when appropriate.
191!-----------------------------------------------------------------------
192!
193 timelb=(time-0.5_dp*dt)/86400.0_dp
194 timeub=(time+0.5_dp*dt)/86400.0_dp
195!
196 DO iobs=nobsstr,nobsend
197 IF ((obstype(iobs).eq.ifield).and. &
198 & ((timelb.le.tobs(iobs)).and.(tobs(iobs).lt.timeub)).and. &
199 & ((xmin.le.xobs(iobs)).and.(xobs(iobs).lt.xmax)).and. &
200 & ((ymin.le.yobs(iobs)).and.(yobs(iobs).lt.ymax))) THEN
201 IF (obstype(iobs).eq.obsstate2type(isubar)) THEN
202 i1=int(xobs(iobs)+0.5_r8) ! 2D U-grid type variable
203 j1=int(yobs(iobs))
204 ELSE IF (obstype(iobs).eq.obsstate2type(isvbar)) THEN
205 i1=int(xobs(iobs)) ! 2D V-grid type variable
206 j1=int(yobs(iobs)+0.5_r8)
207 ELSE
208 i1=int(xobs(iobs)) ! 2D RHO-grid type variable
209 j1=int(yobs(iobs))
210 END IF
211 i2=i1+1
212 j2=j1+1
213 IF (i2.gt.imax) THEN
214 i2=i1 ! Observation at the eastern boundary
215 END IF
216 IF (j2.gt.jmax) THEN
217 j2=j1 ! Observation at the northern boundary
218 END IF
219 p2=real(i2-i1,r8)*(xobs(iobs)-real(i1,r8))
220 q2=real(j2-j1,r8)*(yobs(iobs)-real(j1,r8))
221 p1=1.0_r8-p2
222 q1=1.0_r8-q2
223 hmat(1)=p1*q1
224 hmat(2)=p2*q1
225 hmat(3)=p2*q2
226 hmat(4)=p1*q2
227# ifdef MASKING
228 hmat(1)=hmat(1)*amask(i1,j1)
229 hmat(2)=hmat(2)*amask(i2,j1)
230 hmat(3)=hmat(3)*amask(i2,j2)
231 hmat(4)=hmat(4)*amask(i1,j2)
232 wsum=0.0_r8
233 DO ic=1,4
234 wsum=wsum+hmat(ic)
235 END DO
236 IF (wsum.gt.0.0_r8) THEN
237 wsum=1.0_r8/wsum
238 DO ic=1,4
239 hmat(ic)=hmat(ic)*wsum
240 END DO
241 END IF
242# endif
243 aobs(iobs)=hmat(1)*a(i1,j1)+ &
244 & hmat(2)*a(i2,j1)+ &
245 & hmat(3)*a(i2,j2)+ &
246 & hmat(4)*a(i1,j2)
247# ifdef MASKING
248 IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
249# else
250 obsvetting(iobs)=1.0_r8
251# endif
252 END IF
253 END DO
254
255 RETURN
256 END SUBROUTINE extract_obs2d
257
258# ifdef SOLVE3D
259!
260!***********************************************************************
261 SUBROUTINE extract_obs3d (ng, Imin, Imax, Jmin, Jmax, &
262 & LBi, UBi, LBj, UBj, LBk, UBk, &
263 & ifield, Mobs, NobsSTR, NobsEND, &
264 & Xmin, Xmax, Ymin, Ymax, &
265 & time, dt, &
266 & ObsType, ObsVetting, &
267 & Tobs, Xobs, Yobs, Zobs, &
268 & A, Adepth, &
269# ifdef MASKING
270 & Amask, &
271# endif
272 & Aobs)
273!***********************************************************************
274!
275 USE mod_param
276!
277 USE mod_ncparam, ONLY : isuvel, isvvel
278 USE mod_fourdvar, ONLY : obsstate2type
279!
280! Imported variable declarations.
281!
282 integer, intent(in) :: ng, imin, imax, jmin, jmax
283 integer, intent(in) :: lbi, ubi, lbj, ubj, lbk, ubk
284 integer, intent(in) :: ifield, mobs, nobsstr, nobsend
285!
286 real(r8), intent(in) :: xmin, xmax, ymin, ymax
287 real(dp), intent(in) :: time, dt
288!
289# ifdef ASSUMED_SHAPE
290 integer, intent(in) :: obstype(:)
291
292 real(dp), intent(in) :: tobs(:)
293 real(r8), intent(in) :: xobs(:)
294 real(r8), intent(in) :: yobs(:)
295 real(r8), intent(in) :: a(lbi:,lbj:,lbk:)
296 real(r8), intent(in) :: adepth(lbi:,lbj:,lbk:)
297# ifdef MASKING
298 real(r8), intent(in) :: amask(lbi:,lbj:)
299# endif
300 real(r8), intent(inout) :: obsvetting(:)
301 real(r8), intent(inout) :: zobs(:)
302 real(r8), intent(inout) :: aobs(:)
303# else
304 integer, intent(in) :: obstype(mobs)
305
306 real(dp), intent(in) :: tobs(mobs)
307 real(r8), intent(in) :: xobs(mobs)
308 real(r8), intent(in) :: yobs(mobs)
309 real(r8), intent(in) :: a(lbi:ubi,lbj:ubj,lbk:ubk)
310 real(r8), intent(in) :: adepth(lbi:ubi,lbj:ubj,lbk:ubk)
311# ifdef MASKING
312 real(r8), intent(in) :: amask(lbi:ubi,lbj:ubj)
313# endif
314 real(r8), intent(inout) :: obsvetting(mobs)
315 real(r8), intent(inout) :: zobs(mobs)
316 real(r8), intent(inout) :: aobs(mobs)
317# endif
318!
319! Local variable declarations.
320!
321 integer :: i, ic, iobs, i1, i2, j1, j2, k, k1, k2
322
323 real(dp) :: timelb, timeub
324
325 real(r8) :: zbot, ztop, dz, p1, p2, q1, q2, r1, r2
326 real(r8) :: w11, w12, w21, w22, wsum
327
328 real(r8), dimension(8) :: hmat
329!
330!-----------------------------------------------------------------------
331! Interpolate from requested 3D state field.
332!-----------------------------------------------------------------------
333!
334 timelb=(time-0.5_dp*dt)/86400.0_dp
335 timeub=(time+0.5_dp*dt)/86400.0_dp
336!
337 DO iobs=nobsstr,nobsend
338 IF ((obstype(iobs).eq.ifield).and. &
339 & ((timelb.le.tobs(iobs)).and.(tobs(iobs).lt.timeub)).and. &
340 & ((xmin.le.xobs(iobs)).and.(xobs(iobs).lt.xmax)).and. &
341 & ((ymin.le.yobs(iobs)).and.(yobs(iobs).lt.ymax))) THEN
342 IF (obstype(iobs).eq.obsstate2type(isuvel)) THEN
343 i1=int(xobs(iobs)+0.5_r8) ! 3D U-grid type variable
344 j1=int(yobs(iobs))
345 ELSE IF (obstype(iobs).eq.obsstate2type(isvvel)) THEN
346 i1=int(xobs(iobs)) ! 3D V-grid type variable
347 j1=int(yobs(iobs)+0.5_r8)
348 ELSE
349 i1=int(xobs(iobs)) ! 3D RHO-grid type variable
350 j1=int(yobs(iobs))
351 END IF
352 i2=i1+1
353 j2=j1+1
354 IF (i2.gt.imax) THEN
355 i2=i1 ! Observation at the eastern boundary
356 END IF
357 IF (j2.gt.jmax) THEN
358 j2=j1 ! Observation at the northern boundary
359 END IF
360 p2=real(i2-i1,r8)*(xobs(iobs)-real(i1,r8))
361 q2=real(j2-j1,r8)*(yobs(iobs)-real(j1,r8))
362 p1=1.0_r8-p2
363 q1=1.0_r8-q2
364 w11=p1*q1
365 w21=p2*q1
366 w22=p2*q2
367 w12=p1*q2
368 IF (zobs(iobs).gt.0.0_r8) THEN
369 k1=max(1,int(zobs(iobs))) ! Positions in fractional
370 k2=min(int(zobs(iobs))+1,n(ng)) ! levels
371 r2=real(k2-k1,r8)*(zobs(iobs)-real(k1,r8))
372 r1=1.0_r8-r2
373 ELSE
374 ztop=adepth(i1,j1,n(ng))
375 zbot=adepth(i1,j1,1 )
376 IF (zobs(iobs).ge.ztop) THEN
377 k1=n(ng) ! If shallower, assign to
378 k2=n(ng) ! top grid cell. The
379 r1=1.0_r8 ! observation is located
380 r2=0.0_r8 ! on the upper cell half
381 zobs(iobs)=real(n(ng),r8) ! above its middle depth.
382 ELSE IF (zbot.ge.zobs(iobs)) THEN
383 r1=0.0_r8 ! If deeper, ignore.
384 r2=0.0_r8
385 obsvetting(iobs)=0.0_r8
386 ELSE
387 DO k=n(ng),2,-1 ! Otherwise, interpolate
388 ztop=adepth(i1,j1,k ) ! to fractional level
389 zbot=adepth(i1,j1,k-1)
390 IF ((ztop.gt.zobs(iobs)).and.(zobs(iobs).ge.zbot)) THEN
391 k1=k-1
392 k2=k
393 END IF
394 END DO
395 dz=adepth(i1,j1,k2)-adepth(i1,j1,k1)
396 r2=(zobs(iobs)-adepth(i1,j1,k1))/dz
397 r1=1.0_r8-r2
398 zobs(iobs)=real(k1,r8)+r2 ! overwrite
399 END IF
400 END IF
401 IF ((r1+r2).gt.0.0_r8) THEN
402 hmat(1)=w11*r1
403 hmat(2)=w21*r1
404 hmat(3)=w22*r1
405 hmat(4)=w12*r1
406 hmat(5)=w11*r2
407 hmat(6)=w21*r2
408 hmat(7)=w22*r2
409 hmat(8)=w12*r2
410# ifdef MASKING
411 hmat(1)=hmat(1)*amask(i1,j1)
412 hmat(2)=hmat(2)*amask(i2,j1)
413 hmat(3)=hmat(3)*amask(i2,j2)
414 hmat(4)=hmat(4)*amask(i1,j2)
415 hmat(5)=hmat(5)*amask(i1,j1)
416 hmat(6)=hmat(6)*amask(i2,j1)
417 hmat(7)=hmat(7)*amask(i2,j2)
418 hmat(8)=hmat(8)*amask(i1,j2)
419 wsum=0.0_r8
420 DO ic=1,8
421 wsum=wsum+hmat(ic)
422 END DO
423 IF (wsum.gt.0.0_r8) THEN
424 wsum=1.0_r8/wsum
425 DO ic=1,8
426 hmat(ic)=hmat(ic)*wsum
427 END DO
428 END IF
429# endif
430 aobs(iobs)=hmat(1)*a(i1,j1,k1)+ &
431 & hmat(2)*a(i2,j1,k1)+ &
432 & hmat(3)*a(i2,j2,k1)+ &
433 & hmat(4)*a(i1,j2,k1)+ &
434 & hmat(5)*a(i1,j1,k2)+ &
435 & hmat(6)*a(i2,j1,k2)+ &
436 & hmat(7)*a(i2,j2,k2)+ &
437 & hmat(8)*a(i1,j2,k2)
438# ifdef MASKING
439 IF (wsum.gt.0.0_r8) obsvetting(iobs)=1.0_r8
440# else
441 obsvetting(iobs)=1.0_r8
442# endif
443# ifndef ALLOW_BOTTOM_OBS
444!
445! Reject observations that lie in the lower bottom grid cell (k=1) to
446! avoid clustering due shallowing of bathymetry during smoothing and
447! coarse level half-thickness (-h < Zobs < Adepth(:,:,1)) in deep
448! water.
449!
450 IF ((zobs(iobs).gt.0.0_r8).and.(zobs(iobs).le.1.0_r8)) THEN
451 obsvetting(iobs)=0.0_r8
452 END IF
453# endif
454 END IF
455 END IF
456 END DO
457!
458 RETURN
459 END SUBROUTINE extract_obs3d
460# endif
461#endif
462 END MODULE extract_obs_mod
subroutine, public 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, a, amask, aobs)
subroutine, public 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, a, adepth, amask, aobs)
real(r8), dimension(:), allocatable obsvetting
integer, dimension(:), allocatable obstype
real(r8), dimension(:), allocatable zobs
real(dp), dimension(:), allocatable tobs
real(r8), dimension(:), allocatable xobs
integer, dimension(:), allocatable obsstate2type
real(r8), dimension(:), allocatable yobs
integer, parameter r8
Definition mod_kinds.F:28
integer, parameter dp
Definition mod_kinds.F:25
integer isvvel
integer isvbar
integer isuvel
integer isubar
integer, dimension(:), allocatable n
Definition mod_param.F:479