1MODULE ensembleObservations_mod
2 ! MODULE ensembleObservations_mod (prefix='eob' category='6. High-level data objects')
3 !
4 !:Purpose: Store and manipulate ensemble of quanitites in observation space.
5 ! This module uses the kdtree2 module for efficiently finding the
6 ! nearest observations within the local volume.
7 !
8 use kdTree2_mod
9 use columnData_mod
10 use tovsNL_mod
11 use rttov_types, only: rttov_transmission, rttov_profile
12 use parkind1, only: jpim, jprb
13 use midasMpi_mod
14 use oceanMask_mod
15 use obsSpaceData_mod
16 use randomNumber_mod
17 use mathPhysConstants_mod
18 use utilities_mod
19 use earthConstants_mod
20 use bufr_mod
21 use codePrecision_mod
22 use codtyp_mod
23 use obsfamilylist_mod
24 use varnamelist_mod
25 implicit none
26 save
27 private
28
29 ! public types
30 public :: struct_eob
31
32 ! public procedures
33 public :: eob_init, eob_allocate, eob_deallocate, eob_allGather, eob_getLocalBodyIndices
34 public :: eob_setYb, eob_setYa, eob_setDeterYb, eob_setLatLonObs, eob_setObsErrInv
35 public :: eob_setHPHT, eob_calcAndRemoveMeanYb, eob_setVertLocation, eob_setAssFlag, eob_copy, eob_zero
36 public :: eob_calcRandPert, eob_setSigiSigo, eob_setTypeVertCoord, eob_setSimObsVal
37 public :: eob_backgroundCheck, eob_huberNorm, eob_rejectRadNearSfc, eob_setMeanOMP
38 public :: eob_removeObsNearLand, eob_readFromFiles, eob_writeToFiles
39
40 ! public variables
41 public :: eob_simObsAssim
42
43 integer, parameter :: maxNumLocalObsSearch = 500000
44 integer, external :: get_max_rss
45 logical :: eob_simObsAssim, psvObsAssim
46 integer :: numSimObsFam
47 integer :: numPsvObsFam
48 integer :: numSimCodTyp(ofl_numFamily), numPsvCodTyp(ofl_numFamily)
49 integer :: numSimVarNum(vnl_numvarmax), numPsvVarNum(vnl_numvarmax)
50 integer, allocatable :: simCodTyp(:,:), psvCodTyp(:,:)
51 integer, allocatable :: simVarNum(:,:), psvVarNum(:,:)
52
53 type struct_eob
54 logical :: allocated = .false.
55 logical :: meanRemoved = .false.
56 integer :: numMembers ! number of ensemble members
57 integer :: numObs ! number of observations
58 integer :: fileMemberIndex1 = 1 ! first member number in ensemble set
59 character(len=20) :: typeVertCoord = 'undefined' ! 'logPressure' or 'depth'
60 type(struct_obs), pointer :: obsSpaceData ! pointer to obsSpaceData object
61 real(8), allocatable :: lat(:), lon(:) ! lat/lon of observation
62 real(8), allocatable :: vertLocation(:) ! in ln(pres) or meters, used for localization
63 real(8), allocatable :: obsErrInv(:) ! inverse of obs error variances
64 real(8), allocatable :: obsErrInv_sim(:) ! like obsErrInv, used when simulating observations
65 real(4), allocatable :: Yb_r4(:,:) ! background ensemble perturbation in obs space
66 real(4), allocatable :: Ya_r4(:,:) ! analysis ensemble perturbation in obs space
67 real(4), allocatable :: randPert_r4(:,:) ! unbiased random perturbations with covariance equal to R
68 real(8), allocatable :: meanYb(:) ! ensemble mean background state in obs space
69 real(8), allocatable :: deterYb(:) ! deterministic background state in obs space
70 real(8), allocatable :: obsValue(:) ! the observed value
71 integer, allocatable :: assFlag(:) ! assimilation flag
72 end type struct_eob
73
74 type(kdtree2), pointer :: tree => null()
75
76 ! namelist variables
77 character(len=2) :: simObsFamily(ofl_numFamily) ! observation families for simulation
78 character(len=2) :: psvObsFamily(ofl_numFamily) ! observation families for passive assimilation
79 character(len=codtyp_name_length) :: simCodTypName(ofl_numFamily,codtyp_maxNumber) ! codtyp names for sim. obs families
80 character(len=codtyp_name_length) :: psvCodTypName(ofl_numFamily,codtyp_maxNumber) ! codtyp names for psv. obs families
81 character(len=4) :: simVarName(ofl_numFamily,vnl_numvarmax) ! varName(s) for sim. obs families
82 character(len=4) :: psvVarName(ofl_numFamily,vnl_numvarmax) ! varName(s) for psv. obs families
83 namelist /NAMENSOBS/simObsFamily, psvObsFamily, simCodTypName, psvCodTypName, simVarName, psvVarName
84
85
86CONTAINS
87
88 !--------------------------------------------------------------------------
89 ! eob_init
90 !--------------------------------------------------------------------------
91 subroutine eob_init()
92 !
93 !: Purpose: This subroutine reads the namelist section NAMENSOBS for this module.
94 !
95 implicit none
96
97 ! Locals:
98 integer :: nulnam, ierr, obsfamIndex, codtypIndex, varnumIndex
99 integer, external :: fnom, fclos
100 logical, save :: eob_initialized = .false.
101
102 if (eob_initialized) return
103 write(*,*) 'eob_init: starting'
104 eob_initialized = .true.
105
106 ! default values for namelist variables
107 simObsFamily(:) = ''
108 simCodTypName(:,:) = ''
109 simVarName(:,:) = ''
110 psvObsFamily(:) = ''
111 psvCodTypName(:,:) = ''
112 psvVarName(:,:) = ''
113
114 ! for tracking the number of non-empty chars in namelist variable arrays;
115 ! these are used in loops in various subroutines
116 numSimObsFam = 0
117 numPsvObsFam = 0
118 numSimCodTyp(:) = 0
119 numPsvCodTyp(:) = 0
120 numSimVarNum(:) = 0
121 numPsvVarNum(:) = 0
122
123 ! read namelist
124 if (utl_isNamelistPresent('namensobs','./flnml')) then
125 nulnam=0
126 ierr=fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
127 read(nulnam,nml=namensobs,iostat=ierr)
128 if (ierr /= 0) call utl_abort('eob_init: Error reading namelist namensobs')
129 ierr=fclos(nulnam)
130 do obsfamIndex = 1, ofl_numFamily
131 if (trim(simObsFamily(obsfamIndex)) /= '') then
132 numSimObsFam = numSimObsFam + 1
133 end if
134 if (trim(psvObsFamily(obsfamIndex)) /= '') then
135 numPsvObsFam = numPsvObsFam + 1
136 end if
137 end do
138
139 do obsfamIndex = 1, ofl_numFamily
140 ! Simulation functionality section
141 if (trim(simObsFamily(obsfamIndex)) /= '') then
142 ! simulated observation family specified for current
143 ! obsfamIndex; check to see if any codtyp names specified
144 do codtypIndex = 1, codtyp_maxNumber
145 if (trim(simCodTypName(obsfamIndex,codtypIndex)) /= '') then
146 numSimCodTyp(obsfamIndex) = numSimCodTyp(obsfamIndex) + 1
147 if (.not. allocated(simCodTyp)) then
148 allocate(simCodTyp(numSimObsFam,codtyp_maxNumber))
149 simCodTyp(:,:) = -999
150 end if
151 ! store CodTyp for simulated obs family
152 simCodTyp(obsfamIndex,codtypIndex) = codtyp_get_codtyp(simCodTypName(obsfamIndex,codtypIndex))
153 end if
154 end do ! codtypIndex
155 ! also check to see if any varnames specified
156 do varnumIndex = 1, vnl_numvarmax
157 if (trim(simVarName(obsfamIndex,varnumIndex)) /= '') then
158 numSimVarNum(obsfamIndex) = numSimVarNum(obsfamIndex) + 1
159 if (.not. allocated(simVarNum)) then
160 allocate(simVarNum(numSimObsFam,vnl_numvarmax))
161 simVarNum(:,:) = -999
162 end if
163 ! store VarNum for simulated obs family
164 simVarNum(obsfamIndex,varnumIndex) = vnl_varnumFromVarName(simVarName(obsfamIndex,varnumIndex))
165 end if
166 end do ! varnumIndex
167 end if ! simObsFamily
168
169 ! Passive functionality section
170 if (trim(psvObsFamily(obsfamIndex)) /= '') then
171 ! passive observation family specified for current
172 ! obsfamIndex; check to see if any codtyp names specified
173 do codtypIndex = 1, codtyp_maxNumber
174 if (trim(psvCodTypName(obsfamIndex,codtypIndex)) /= '') then
175 numPsvCodTyp(obsfamIndex) = numPsvCodTyp(obsfamIndex) + 1
176 if (.not. allocated(psvCodTyp)) then
177 allocate(psvCodTyp(numPsvObsFam,codtyp_maxNumber))
178 psvCodTyp(:,:) = -999
179 end if
180 ! store CodTyp for passive obs family
181 psvCodTyp(obsfamIndex,codtypIndex) = codtyp_get_codtyp(psvCodTypName(obsfamIndex,codtypIndex))
182 end if
183 end do ! codtypIndex
184 ! also check to see if any varnames specified
185 do varnumIndex = 1, vnl_numvarmax
186 if (trim(psvVarName(obsfamIndex,varnumIndex)) /= '') then
187 numPsvVarNum(obsfamIndex) = numPsvVarNum(obsfamIndex) + 1
188 if (.not. allocated(PsvVarNum)) then
189 allocate(psvVarNum(numPsvObsFam,vnl_numvarmax))
190 psvVarNum(:,:) = -999
191 end if
192 ! store VarNum for passive obs family
193 psvVarNum(obsfamIndex,varnumIndex) = vnl_varnumFromVarName(psvVarName(obsfamIndex,varnumIndex))
194 end if
195 end do ! varnumIndex
196 end if ! psvObsFamily
197
198 end do ! obsFamIndex
199 else
200 write(*,*)
201 write(*,*) 'eob_init: namensobs is missing in the namelist. The default value will be taken.'
202 end if
203
204 eob_simObsAssim = numSimObsFam > 0
205 psvObsAssim = numPsvObsFam > 0
206
207 write(*,*) 'eob_init: eob_simObsAssim = ', eob_simObsAssim
208 write(*,*) 'eob_init: psvObsAssim = ', psvObsAssim
209
210 end subroutine eob_init
211
212 !--------------------------------------------------------------------------
213 ! eob_allocate
214 !--------------------------------------------------------------------------
215 subroutine eob_allocate(ensObs, numMembers, numObs, obsSpaceData, &
216 fileMemberIndex1_opt)
217 !
218 !:Purpose: Allocate an ensObs object
219 !
220 implicit none
221
222 ! Arguments:
223 type(struct_eob) , intent(inout) :: ensObs
224 integer , intent(in) :: numMembers
225 integer , intent(in) :: numObs
226 type(struct_obs), target, intent(in) :: obsSpaceData
227 integer, optional , intent(in) :: fileMemberIndex1_opt
228
229 if (ensObs%allocated) then
230 write(*,*) 'eob_allocate: this object is already allocated, deallocating first.'
231 call eob_deallocate(ensObs)
232 end if
233
234 if (present(fileMemberIndex1_opt)) ensObs%fileMemberIndex1 = fileMemberIndex1_opt
235
236 call eob_init()
237
238 ensObs%obsSpaceData => obsSpaceData
239 ensObs%numMembers = numMembers
240 ensObs%numObs = numObs
241
242 allocate(ensObs%lat(ensObs%numObs))
243 allocate(ensObs%lon(ensObs%numObs))
244 allocate(ensObs%vertLocation(ensObs%numObs))
245 allocate(ensObs%obsValue(ensObs%numObs))
246 allocate(ensObs%obsErrInv(ensObs%numObs))
247 if (eob_simObsAssim) allocate(ensObs%obsErrInv_sim(ensObs%numObs))
248 allocate(ensObs%Yb_r4(ensObs%numMembers,ensObs%numObs))
249 allocate(ensObs%meanYb(ensObs%numObs))
250 allocate(ensObs%deterYb(ensObs%numObs))
251 allocate(ensObs%assFlag(ensObs%numObs))
252
253 ensObs%allocated = .true.
254
255 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
256
257 end subroutine eob_allocate
258
259 !--------------------------------------------------------------------------
260 ! eob_deallocate
261 !--------------------------------------------------------------------------
262 subroutine eob_deallocate(ensObs)
263 implicit none
264
265 ! Arguments:
266 type(struct_eob), intent(inout) :: ensObs
267
268 if (.not. ensObs%allocated) return
269
270 deallocate(ensObs%lat)
271 deallocate(ensObs%lon)
272 deallocate(ensObs%vertLocation)
273 deallocate(ensObs%obsValue)
274 deallocate(ensObs%obsErrInv)
275 if (allocated(ensObs%obsErrInv_sim)) deallocate(ensObs%obsErrInv_sim)
276 deallocate(ensObs%Yb_r4)
277 if (allocated(ensObs%Ya_r4)) deallocate(ensObs%Ya_r4)
278 if (allocated(ensObs%randPert_r4)) deallocate(ensObs%randPert_r4)
279 deallocate(ensObs%meanYb)
280 deallocate(ensObs%deterYb)
281 deallocate(ensObs%assFlag)
282
283 ensObs%allocated = .false.
284
285 end subroutine eob_deallocate
286
287 !--------------------------------------------------------------------------
288 ! eob_zero
289 !--------------------------------------------------------------------------
290 subroutine eob_zero(ensObs)
291 !
292 !:Purpose: Initialize an ensObs object to zero
293 !
294 implicit none
295
296 ! Arguments:
297 type(struct_eob), intent(inout) :: ensObs
298
299 if ( .not.ensObs%allocated ) then
300 call utl_abort('eob_zero: this object is not allocated')
301 end if
302
303 ensObs%lat(:) = 0.0d0
304 ensObs%lon(:) = 0.0d0
305 ensObs%vertLocation(:) = 0.0d0
306 ensObs%obsValue(:) = 0.0d0
307 ensObs%obsErrInv(:) = 0.0d0
308 if (allocated(ensObs%obsErrInv_sim)) ensObs%obsErrInv_sim(:) = 0.0
309 ensObs%Yb_r4(:,:) = 0.0
310 if (allocated(ensObs%Ya_r4)) ensObs%Ya_r4(:,:) = 0.0
311 if (allocated(ensObs%randPert_r4)) ensObs%randPert_r4(:,:) = 0.0
312 ensObs%meanYb(:) = 0.0d0
313 ensObs%deterYb(:) = 0.0d0
314 ensObs%assFlag(:) = 0
315
316 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
317
318 end subroutine eob_zero
319
320 !--------------------------------------------------------------------------
321 ! eob_setTypeVertCoord
322 !--------------------------------------------------------------------------
323 subroutine eob_setTypeVertCoord(ensObs, typeVertCoord)
324 !
325 !:Purpose: Set the type of vertical coordinate ('logPressure' or 'depth').
326 !
327 implicit none
328
329 ! Arguments:
330 type(struct_eob), intent(inout) :: ensObs
331 character(len=*), intent(in) :: typeVertCoord
332
333 if ( trim(typeVertCoord) /= 'logPressure' .and. &
334 trim(typeVertCoord) /= 'depth' ) then
335 write(*,*) 'eob_setTypeVertCoord: typeVertCoord = ', trim(typeVertCoord)
336 call utl_abort('eob_setTypeVertCoord: Unknown type of vertical coordinate')
337 end if
338
339 ensObs%typeVertCoord = typeVertCoord
340
341 end subroutine eob_setTypeVertCoord
342
343 !--------------------------------------------------------------------------
344 ! eob_clean (private routine)
345 !--------------------------------------------------------------------------
346 subroutine eob_clean(ensObs,ensObsClean)
347 !
348 !:Purpose: Remove all obs from the ensObs object that are not
349 ! flagged for assimilation. Put the cleaned result in the
350 ! locally created output object.
351 !
352 implicit none
353
354 ! Arguments:
355 type(struct_eob), intent(inout) :: ensObs
356 type(struct_eob), intent(out) :: ensObsClean
357
358 ! Locals:
359 integer :: obsIndex, obsCleanIndex, numObsClean
360
361 call eob_setAssFlag(ensObs)
362
363 numObsClean = 0
364 do obsIndex = 1, ensObs%numObs
365 if (ensObs%assFlag(obsIndex) == 1) numObsClean = numObsClean + 1
366 end do
367
368 write(*,*) 'eob_clean: reducing numObs from ', ensObs%numObs, ' to ', numObsClean
369 call eob_allocate(ensObsClean, ensObs%numMembers, numObsClean, ensObs%obsSpaceData)
370 if (allocated(ensObs%Ya_r4)) then
371 allocate(ensObsClean%Ya_r4(ensObs%numMembers,numObsClean))
372 end if
373 if (allocated(ensObs%randPert_r4)) then
374 allocate(ensObsClean%randPert_r4(ensObs%numMembers,numObsClean))
375 end if
376
377 obsCleanIndex = 0
378 do obsIndex = 1, ensObs%numObs
379 if (ensObs%assFlag(obsIndex) == 1) then
380 obsCleanIndex = obsCleanIndex + 1
381 ensObsClean%lat(obsCleanIndex) = ensObs%lat(obsIndex)
382 ensObsClean%lon(obsCleanIndex) = ensObs%lon(obsIndex)
383 ensObsClean%vertLocation(obsCleanIndex) = ensObs%vertLocation(obsIndex)
384 ensObsClean%obsErrInv(obsCleanIndex) = ensObs%obsErrInv(obsIndex)
385 if (allocated(ensObs%obsErrInv_sim)) then
386 ensObsClean%obsErrInv_sim(obsCleanIndex) = ensObs%obsErrInv_sim(obsIndex)
387 end if
388 ensObsClean%Yb_r4(:,obsCleanIndex) = ensObs%Yb_r4(:,obsIndex)
389 if (allocated(ensObs%Ya_r4)) then
390 ensObsClean%Ya_r4(:,obsCleanIndex) = ensObs%Ya_r4(:,obsIndex)
391 end if
392 if (allocated(ensObs%randPert_r4)) then
393 ensObsClean%randPert_r4(:,obsCleanIndex) = ensObs%randPert_r4(:,obsIndex)
394 end if
395 ensObsClean%meanYb(obsCleanIndex) = ensObs%meanYb(obsIndex)
396 ensObsClean%deterYb(obsCleanIndex) = ensObs%deterYb(obsIndex)
397 ensObsClean%obsValue(obsCleanIndex) = ensObs%obsValue(obsIndex)
398 ensObsClean%assFlag(obsCleanIndex) = ensObs%assFlag(obsIndex)
399 end if
400 end do
401
402 end subroutine eob_clean
403
404 !--------------------------------------------------------------------------
405 ! eob_copy
406 !--------------------------------------------------------------------------
407 subroutine eob_copy(ensObsIn,ensObsOut)
408 implicit none
409
410 ! Arguments:
411 type(struct_eob), intent(in) :: ensObsIn
412 type(struct_eob), intent(inout) :: ensObsOut
413
414 ensObsOut%lat(:) = ensObsIn%lat(:)
415 ensObsOut%lon(:) = ensObsIn%lon(:)
416 ensObsOut%vertLocation(:) = ensObsIn%vertLocation(:)
417 ensObsOut%obsErrInv(:) = ensObsIn%obsErrInv(:)
418 if (allocated(ensObsIn%obsErrInv_sim)) then
419 ensObsOut%obsErrInv(:) = ensObsIn%obsErrInv_sim(:)
420 end if
421 ensObsOut%Yb_r4(:,:) = ensObsIn%Yb_r4(:,:)
422 if (allocated(ensObsIn%Ya_r4)) then
423 allocate( ensObsOut%Ya_r4(ensObsIn%numMembers,ensObsIn%numObs))
424 ensObsOut%Ya_r4(:,:) = ensObsIn%Ya_r4(:,:)
425 end if
426 if (allocated(ensObsIn%randPert_r4)) then
427 allocate(ensObsOut%randPert_r4(ensObsIn%numMembers,ensObsIn%numObs))
428 ensObsOut%randPert_r4(:,:) = ensObsIn%randPert_r4(:,:)
429 end if
430 ensObsOut%meanYb(:) = ensObsIn%meanYb(:)
431 ensObsOut%deterYb(:) = ensObsIn%deterYb(:)
432 ensObsOut%obsValue(:) = ensObsIn%obsValue(:)
433 ensObsOut%assFlag(:) = ensObsIn%assFlag(:)
434 ensObsOut%typeVertCoord = ensObsIn%typeVertCoord
435
436 end subroutine eob_copy
437
438 !--------------------------------------------------------------------------
439 ! eob_allGather
440 !--------------------------------------------------------------------------
441 subroutine eob_allGather(ensObs,ensObs_mpiglobal)
442 !
443 !:Purpose: Collect obs information distributed over all mpi tasks and
444 ! make it available on all mpi tasks. The output ensObs object
445 ! will be allocated within this subroutine.
446 !
447 implicit none
448
449 ! Arguments:
450 type(struct_eob), intent(inout) :: ensObs
451 type(struct_eob), intent(out) :: ensObs_mpiglobal
452
453 ! Locals:
454 type(struct_eob) :: ensObsClean
455 integer :: ierr, procIndex, memberIndex, numObs_mpiglobal
456 integer :: allNumObs(mmpi_nprocs), displs(mmpi_nprocs)
457
458 write(*,*) 'eob_allGather: starting'
459 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
460
461 call utl_tmg_start(10,'--Observations')
462 call utl_tmg_start(24,'----Eob_AllGather')
463
464 ! refresh assimilation flag and then clean ensObs before communicating and writing
465 call eob_setAssFlag(ensObs)
466 call eob_clean(ensObs,ensObsClean)
467
468 call rpn_comm_allgather(ensObsClean%numObs, 1, 'mpi_integer', &
469 allNumObs, 1, 'mpi_integer', &
470 'GRID', ierr)
471 numObs_mpiglobal = sum(allNumObs(:))
472
473 if (ensObs_mpiglobal%allocated) then
474 call utl_abort('eob_allGather: output ensObs object must not be already allocated')
475 end if
476 call eob_allocate(ensObs_mpiglobal, ensObsClean%numMembers, numObs_mpiglobal, ensObsClean%obsSpaceData, &
477 fileMemberIndex1_opt=ensObs%fileMemberIndex1)
478 if (allocated(ensObsClean%Ya_r4)) then
479 allocate(ensObs_mpiglobal%Ya_r4(ensObsClean%numMembers,numObs_mpiglobal))
480 end if
481 if (allocated(ensObsClean%randPert_r4)) then
482 allocate(ensObs_mpiglobal%randPert_r4(ensObsClean%numMembers,numObs_mpiglobal))
483 end if
484 ensObs_mpiglobal%typeVertCoord = ensObsClean%typeVertCoord
485
486 if (mmpi_myid == 0) then
487 displs(1) = 0
488 do procIndex = 2, mmpi_nprocs
489 displs(procIndex) = displs(procIndex-1) + allNumObs(procIndex-1)
490 end do
491 else
492 displs(:) = 0
493 end if
494
495 call rpn_comm_gatherv(ensObsClean%lat, ensObsClean%numObs, 'mpi_real8', &
496 ensObs_mpiglobal%lat, allNumObs, displs, 'mpi_real8', &
497 0, 'GRID', ierr)
498 call rpn_comm_gatherv(ensObsClean%lon, ensObsClean%numObs, 'mpi_real8', &
499 ensObs_mpiglobal%lon, allNumObs, displs, 'mpi_real8', &
500 0, 'GRID', ierr)
501 call rpn_comm_gatherv(ensObsClean%vertLocation, ensObsClean%numObs, 'mpi_real8', &
502 ensObs_mpiglobal%vertLocation, allNumObs, displs, 'mpi_real8', &
503 0, 'GRID', ierr)
504 call rpn_comm_gatherv(ensObsClean%obsValue, ensObsClean%numObs, 'mpi_real8', &
505 ensObs_mpiglobal%obsValue, allNumObs, displs, 'mpi_real8', &
506 0, 'GRID', ierr)
507 call rpn_comm_gatherv(ensObsClean%obsErrInv, ensObsClean%numObs, 'mpi_real8', &
508 ensObs_mpiglobal%obsErrInv, allNumObs, displs, 'mpi_real8', &
509 0, 'GRID', ierr)
510 call rpn_comm_gatherv(ensObsClean%meanYb, ensObsClean%numObs, 'mpi_real8', &
511 ensObs_mpiglobal%meanYb, allNumObs, displs, 'mpi_real8', &
512 0, 'GRID', ierr)
513 call rpn_comm_gatherv(ensObsClean%deterYb, ensObsClean%numObs, 'mpi_real8', &
514 ensObs_mpiglobal%deterYb, allNumObs, displs, 'mpi_real8', &
515 0, 'GRID', ierr)
516 call rpn_comm_gatherv(ensObsClean%assFlag, ensObsClean%numObs, 'mpi_integer', &
517 ensObs_mpiglobal%assFlag, allNumObs, displs, 'mpi_integer', &
518 0, 'GRID', ierr)
519 if (allocated(ensObsClean%obsErrInv_sim)) then
520 call rpn_comm_gatherv(ensObsClean%obsErrInv_sim, ensObsClean%numObs, 'mpi_real8', &
521 ensObs_mpiglobal%obsErrInv_sim, allNumObs, displs, 'mpi_real8', &
522 0, 'GRID', ierr)
523 end if
524 do memberIndex = 1, ensObsClean%numMembers
525 call rpn_comm_gatherv(ensObsClean%Yb_r4(memberIndex,:), ensObsClean%numObs, 'mpi_real4', &
526 ensObs_mpiglobal%Yb_r4(memberIndex,:), allNumObs, displs, 'mpi_real4', &
527 0, 'GRID', ierr)
528 if (allocated(ensObsClean%Ya_r4)) then
529 call rpn_comm_gatherv(ensObsClean%Ya_r4(memberIndex,:), ensObsClean%numObs, 'mpi_real4', &
530 ensObs_mpiglobal%Ya_r4(memberIndex,:), allNumObs, displs, 'mpi_real4', &
531 0, 'GRID', ierr)
532 end if
533 if (allocated(ensObsClean%randPert_r4)) then
534 call rpn_comm_gatherv(ensObsClean%randPert_r4(memberIndex,:), ensObsClean%numObs, 'mpi_real4', &
535 ensObs_mpiglobal%randPert_r4(memberIndex,:), allNumObs, displs, 'mpi_real4', &
536 0, 'GRID', ierr)
537 end if
538 end do
539
540 call rpn_comm_bcast(ensObs_mpiglobal%lat, ensObs_mpiglobal%numObs, 'mpi_real8', &
541 0, 'GRID', ierr)
542 call rpn_comm_bcast(ensObs_mpiglobal%lon, ensObs_mpiglobal%numObs, 'mpi_real8', &
543 0, 'GRID', ierr)
544 call rpn_comm_bcast(ensObs_mpiglobal%vertLocation, ensObs_mpiglobal%numObs, 'mpi_real8', &
545 0, 'GRID', ierr)
546 call rpn_comm_bcast(ensObs_mpiglobal%obsValue, ensObs_mpiglobal%numObs, 'mpi_real8', &
547 0, 'GRID', ierr)
548 call rpn_comm_bcast(ensObs_mpiglobal%obsErrInv, ensObs_mpiglobal%numObs, 'mpi_real8', &
549 0, 'GRID', ierr)
550 if (allocated(ensObs_mpiglobal%obsErrInv_sim)) then
551 call rpn_comm_bcast(ensObs_mpiglobal%obsErrInv_sim, ensObs_mpiglobal%numObs, 'mpi_real8', &
552 0, 'GRID', ierr)
553 end if
554 call rpn_comm_bcast(ensObs_mpiglobal%meanYb, ensObs_mpiglobal%numObs, 'mpi_real8', &
555 0, 'GRID', ierr)
556 call rpn_comm_bcast(ensObs_mpiglobal%deterYb, ensObs_mpiglobal%numObs, 'mpi_real8', &
557 0, 'GRID', ierr)
558 call rpn_comm_bcast(ensObs_mpiglobal%assFlag, ensObs_mpiglobal%numObs, 'mpi_integer', &
559 0, 'GRID', ierr)
560 do memberIndex = 1, ensObsClean%numMembers
561 call rpn_comm_bcast(ensObs_mpiglobal%Yb_r4(memberIndex,:), ensObs_mpiglobal%numObs, 'mpi_real4', &
562 0, 'GRID', ierr)
563 if (allocated(ensObs_mpiglobal%Ya_r4)) then
564 call rpn_comm_bcast(ensObs_mpiglobal%Ya_r4(memberIndex,:), ensObs_mpiglobal%numObs, 'mpi_real4', &
565 0, 'GRID', ierr)
566 end if
567 if (allocated(ensObs_mpiglobal%randPert_r4)) then
568 call rpn_comm_bcast(ensObs_mpiglobal%randPert_r4(memberIndex,:), ensObs_mpiglobal%numObs, 'mpi_real4', &
569 0, 'GRID', ierr)
570 end if
571 end do
572
573 call eob_deallocate(ensObsClean)
574
575 write(*,*) 'eob_allGather: total number of obs to be assimilated =', sum(ensObs_mpiglobal%assFlag(:))
576
577 call utl_tmg_stop(24)
578 call utl_tmg_stop(10)
579
580 write(*,*) 'eob_allGather: finished'
581 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
582
583 end subroutine eob_allGather
584
585 !--------------------------------------------------------------------------
586 ! eob_writeToFiles
587 !--------------------------------------------------------------------------
588 subroutine eob_writeToFiles(ensObs, outputFilenamePrefix, writeObsInfo, &
589 numGroupsToDivideMembers_opt, &
590 maxNumMembersPerGroup_opt)
591 !
592 !:Purpose: Write the contents of an ensObs mpi local object to files
593 !
594 implicit none
595
596 ! Arguments:
597 type(struct_eob), intent(in) :: ensObs
598 character(len=*), intent(in) :: outputFilenamePrefix
599 logical, intent(in) :: writeObsInfo
600 integer, optional, intent(in) :: numGroupsToDivideMembers_opt
601 integer, optional, intent(in) :: maxNumMembersPerGroup_opt
602
603 ! Locals:
604 integer :: unitNum, ierr, obsIndex, memberIndex
605 integer :: obsVcoCode(ensObs%numObs), obsAssFlag(ensObs%numObs)
606 integer :: obsFlag(ensObs%numObs)
607 integer, allocatable :: memberIndexArray(:)
608 character(len=40) :: fileName
609 character(len=4) :: myidxStr, myidyStr
610 character(len=30) :: fileNameExtention
611 integer :: fnom, fclos
612 logical :: fileExists
613
614 if (.not. ensObs%allocated) then
615 call utl_abort('eob_writeToFiles: this object is not allocated')
616 end if
617
618 call obs_extractObsIntBodyColumn(obsVcoCode, ensObs%obsSpaceData, OBS_VCO)
619 call obs_extractObsIntBodyColumn(obsAssFlag, ensObs%obsSpaceData, OBS_ASS)
620 call obs_extractObsIntBodyColumn(obsFlag, ensObs%obsSpaceData, OBS_FLG)
621
622 write(myidxStr,'(I4.4)') (mmpi_myidx + 1)
623 write(myidyStr,'(I4.4)') (mmpi_myidy + 1)
624 fileNameExtention = trim(myidxStr) // '_' // trim(myidyStr)
625
626 ! write observation info to a file
627 if (writeObsInfo) then
628 fileName = 'eob_obsInfo_' // trim(fileNameExtention)
629 write(*,*) 'eob_writeToFiles: writing ',trim(filename)
630 inquire(file=trim(fileName),exist=fileExists)
631 if ( fileExists ) then
632 call utl_abort('eob_writeToFiles: file should not exist')
633 end if
634
635 unitNum = 0
636 ierr = fnom(unitNum, fileName, 'FTN+SEQ+UNF+R/W', 0)
637 write(unitNum) ensObs%numMembers, ensObs%numObs
638 write(unitNum) (ensObs%lat(obsIndex), obsIndex = 1, ensObs%numObs)
639 write(unitNum) (ensObs%lon(obsIndex), obsIndex = 1, ensObs%numObs)
640 write(unitNum) (obsVcoCode(obsIndex), obsIndex = 1, ensObs%numObs)
641 write(unitNum) (ensObs%obsValue(obsIndex), obsIndex = 1, ensObs%numObs)
642 write(unitNum) (obsAssFlag(obsIndex), obsIndex = 1, ensObs%numObs)
643 write(unitNum) (obsFlag(obsIndex), obsIndex = 1, ensObs%numObs)
644 ierr = fclos(unitNum)
645 end if
646
647 ! get memberIndex in the full ensemble set
648 allocate(memberIndexArray(ensObs%numMembers))
649 call getMemberIndexInFullEnsSet(ensObs, memberIndexArray, &
650 numGroupsToDivideMembers_opt=numGroupsToDivideMembers_opt, &
651 maxNumMembersPerGroup_opt=maxNumMembersPerGroup_opt)
652
653 ! Open file and write ensObs%Yb for all the members to one file
654 fileName = trim(outputFilenamePrefix) // '_' // trim(fileNameExtention)
655 write(*,*) 'eob_writeToFiles: writing ',trim(filename)
656 inquire(file=trim(fileName),exist=fileExists)
657 if (fileExists) then
658 call utl_abort('eob_writeToFiles: file should not exist')
659 end if
660
661 unitNum = 0
662 ierr = fnom(unitNum, fileName, 'FTN+SEQ+UNF+R/W', 0)
663 write(unitNum) ensObs%numMembers
664 write(unitNum) (memberIndexArray(memberIndex), memberIndex = 1, ensObs%numMembers)
665 do memberIndex = 1, ensObs%numMembers
666 if (mmpi_myid == 0) then
667 write(*,*) 'eob_writeToFiles: fileMemberIndex1=', ensObs%fileMemberIndex1, &
668 ', memberIndex=', memberIndex, &
669 ', memberIndex in full ensemble set=', memberIndexArray(memberIndex)
670 end if
671 write(unitNum) (ensObs%Yb_r4(memberIndex,obsIndex), obsIndex = 1, ensObs%numObs)
672 end do
673 ierr = fclos(unitNum)
674
675 deallocate(memberIndexArray)
676
677 end subroutine eob_writeToFiles
678
679 !--------------------------------------------------------------------------
680 ! eob_readFromFiles
681 !--------------------------------------------------------------------------
682 subroutine eob_readFromFiles(ensObs, numMembersToRead, inputFilenamePrefix, &
683 readObsInfo)
684 !
685 !:Purpose: Read mpi local ensObs%Yb object from file. Several files in separate subdirectories
686 ! can be read. Some examples of path+filename are:
687 ! ensObs_0001/eob_HX_0001_0001
688 ! ensObs_0002/eob_HX_0001_0001
689 !
690 implicit none
691
692 ! Arguments:
693 type(struct_eob), intent(inout) :: ensObs
694 integer , intent(in) :: numMembersToRead
695 character(len=*), intent(in) :: inputFilenamePrefix
696 logical, intent(in) :: readObsInfo
697
698 ! Locals:
699 real(8) :: latFromFile(ensObs%numObs), lonFromFile(ensObs%numObs)
700 real(8) :: obsValueFromFile(ensObs%numObs)
701 integer :: obsVcoCode(ensObs%numObs), obsVcoCodeFromFile(ensObs%numObs)
702 integer :: obsFlag(ensObs%numObs), assFlagFrom1File(ensObs%numObs)
703 integer :: assFlagFromAllFiles(ensObs%numObs)
704 integer :: unitNum, ierr, memberIndex, obsIndex, numObsFromFile
705 integer :: numMembersFromFile, fnom, fclos
706 integer :: fileIndex, numMembersAlreadyRead
707 integer, allocatable :: memberIndexFromFile(:)
708 logical :: fileExists
709 character(len=256) :: fileName
710 character(len=100) :: fileBaseName
711 character(len=4) :: myidxStr, myidyStr
712 character(len=3) :: fileIndexStr
713 character(len=30) :: fileNameExtention
714
715 if ( .not. ensObs%allocated ) then
716 call utl_abort('eob_readFromFiles: this object is not allocated')
717 end if
718
719 call obs_extractObsIntBodyColumn(obsVcoCode, ensObs%obsSpaceData, OBS_VCO)
720
721 write(myidxStr,'(I4.4)') (mmpi_myidx + 1)
722 write(myidyStr,'(I4.4)') (mmpi_myidy + 1)
723 fileNameExtention = trim(myidxStr) // '_' // trim(myidyStr)
724
725 if (readObsInfo) then
726 ! loop on all file from different directories, read obsInfo and check they match ensObs
727 fileBaseName = 'eob_obsInfo_' // trim(fileNameExtention)
728
729 fileIndex = 0
730 numMembersAlreadyRead = 0
731 do while (numMembersAlreadyRead < numMembersToRead)
732 fileIndex = fileIndex + 1
733 write(fileIndexStr,'(i3.3)') fileIndex
734 fileName = './ensObs_' // fileIndexStr // '/' // fileBaseName
735
736 write(*,*) 'eob_readFromFiles: reading ',trim(fileName)
737 inquire(file=trim(fileName),exist=fileExists)
738 if (.not. fileExists) then
739 write(*,*) 'fileName=', fileName
740 call utl_abort('eob_readFromFiles: file does not exist')
741 end if
742
743 unitNum = 0
744 ierr = fnom(unitNum,trim(fileName),'FTN+SEQ+UNF',0)
745 read(unitNum) numMembersFromFile, numObsFromFile
746 if (ensObs%numObs /= numObsFromFile) then
747 call utl_abort('eob_readFromFiles: ensObs%numObs does not match with that of file')
748 end if
749
750 read(unitNum) (latFromFile(obsIndex), obsIndex = 1, ensObs%numObs)
751 read(unitNum) (lonFromFile(obsIndex), obsIndex = 1, ensObs%numObs)
752 read(unitNum) (obsVcoCodeFromFile(obsIndex), obsIndex = 1, ensObs%numObs)
753 read(unitNum) (obsValueFromFile(obsIndex), obsIndex = 1, ensObs%numObs)
754
755 if (maxval(abs(latFromFile(:) - ensObs%lat(:))) > 1.0d-5 .or. &
756 maxval(abs(lonFromFile(:) - ensObs%lon(:))) > 1.0d-5 .or. &
757 maxval(abs(obsValueFromFile(:) - ensObs%obsValue(:))) > 1.0d-7 .or. &
758 .not. all(obsVcoCodeFromFile(:) == obsVcoCode(:))) then
759
760 call utl_abort('eob_readFromFiles: obsInfo file do not match ensObs')
761 end if
762
763 ! Read assimilation flag for of all files and apply a "logical or" to get the value
764 ! to put in obsSpaceData. Read obs flag only on the first file.
765 read(unitNum) (assFlagFrom1File(obsIndex), obsIndex = 1, ensObs%numObs)
766 if (numMembersAlreadyRead == 0) then
767 read(unitNum) (obsFlag(obsIndex), obsIndex = 1, ensObs%numObs)
768 end if
769
770 if (numMembersAlreadyRead == 0) assFlagFromAllFiles(:) = assFlagFrom1File(:)
771 where (assFlagFrom1File(:) == obs_notAssimilated .and. numMembersAlreadyRead > 0)
772 assFlagFromAllFiles(:) = obs_notAssimilated
773 end where
774
775 ierr = fclos(unitNum)
776
777 numMembersAlreadyRead = numMembersAlreadyRead + numMembersFromFile
778 end do
779
780 ! update assimilation flag in obsSpaceData
781 do obsIndex = 1, ensObs%numObs
782 ! skip this obs it is already set to be assimilated
783 if (assFlagFromAllFiles(obsIndex) == obs_assimilated) cycle
784
785 call obs_bodySet_i(ensObs%obsSpaceData, OBS_ASS, obsIndex, obs_notAssimilated)
786 call obs_bodySet_i(ensObs%obsSpaceData, OBS_FLG, obsIndex, obsFlag(obsIndex))
787 end do
788 end if
789
790 ! loop on all files from different directories to read ensObs%Yb for all members
791 fileBaseName = trim(inputFilenamePrefix) // '_' // trim(fileNameExtention)
792
793 fileIndex = 0
794 numMembersAlreadyRead = 0
795 do while (numMembersAlreadyRead < numMembersToRead)
796 fileIndex = fileIndex + 1
797 write(fileIndexStr,'(i3.3)') fileIndex
798 fileName = './ensObs_' // fileIndexStr // '/' // fileBaseName
799
800 write(*,*) 'eob_readFromFiles: reading ',trim(fileName)
801 inquire(file=trim(fileName),exist=fileExists)
802 if (.not. fileExists) then
803 write(*,*) 'fileName=', fileName
804 call utl_abort('eob_readFromFiles: file does not exist')
805 end if
806
807 unitNum = 0
808 ierr = fnom(unitNum,trim(fileName),'FTN+SEQ+UNF',0)
809 read(unitNum) numMembersFromFile
810 allocate(memberIndexFromFile(numMembersFromFile))
811 read(unitNum) (memberIndexFromFile(memberIndex), memberIndex = 1, numMembersFromFile)
812 do memberIndex = 1, numMembersFromFile
813 read(unitNum) (ensObs%Yb_r4(memberIndexFromFile(memberIndex),obsIndex), obsIndex = 1, ensObs%numObs)
814 end do
815 ierr = fclos(unitNum)
816
817 deallocate(memberIndexFromFile)
818
819 numMembersAlreadyRead = numMembersAlreadyRead + numMembersFromFile
820 end do
821
822 end subroutine eob_readFromFiles
823
824 !--------------------------------------------------------------------------
825 ! eob_getLocalBodyIndices
826 !--------------------------------------------------------------------------
827 function eob_getLocalBodyIndices(ensObs,localBodyIndices,distances,lat,lon,vertLocation, &
828 hLocalize,vLocalize,numLocalObsFound) result(numLocalObs)
829 !
830 !:Purpose: Return a list of values of bodyIndex for all observations within
831 ! the local volume around the specified lat/lon used for assimilation
832 ! (as defined by h/vLocalize). The kdtree2 module is used to efficiently
833 ! perform this task. The kdtree itself is constructed on the first call.
834 !
835 implicit none
836
837 ! Arguments:
838 type(struct_eob), intent(in) :: ensObs
839 integer , intent(out) :: localBodyIndices(:)
840 real(8) , intent(out) :: distances(:)
841 real(8) , intent(in) :: lat, lon, vertLocation, hLocalize, vLocalize
842 integer , intent(out) :: numLocalObsFound
843 ! Result:
844 integer :: numLocalObs ! function output
845
846 ! Locals:
847 integer :: bodyIndex, numLocalObsFoundSearch, maxNumLocalObs, localObsIndex
848 real(8) :: distance
849 real(kdkind), allocatable :: positionArray(:,:)
850 type(kdtree2_result), allocatable :: searchResults(:)
851 real(kdkind) :: maxRadius
852 real(kdkind) :: refPosition(3)
853
854 ! create the kdtree on the first call
855 if (.not. associated(tree)) then
856 write(*,*) 'eob_getLocalBodyIndices: start creating kdtree'
857 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
858 allocate(positionArray(3,ensObs%numObs))
859 do bodyIndex = 1, ensObs%numObs
860 positionArray(1,bodyIndex) = ec_ra * sin(ensObs%lon(bodyIndex)) * cos(ensObs%lat(bodyIndex))
861 positionArray(2,bodyIndex) = ec_ra * cos(ensObs%lon(bodyIndex)) * cos(ensObs%lat(bodyIndex))
862 positionArray(3,bodyIndex) = ec_ra * sin(ensObs%lat(bodyIndex))
863 end do
864 tree => kdtree2_create(positionArray, sort=.true., rearrange=.true.)
865 write(*,*) 'eob_getLocalBodyIndices: done creating kdtree'
866 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
867 end if
868
869 ! do the search
870 maxNumLocalObs = size(localBodyIndices)
871 maxRadius = hLocalize**2
872 refPosition(1) = ec_ra * sin(lon) * cos(lat)
873 refPosition(2) = ec_ra * cos(lon) * cos(lat)
874 refPosition(3) = ec_ra * sin(lat)
875 allocate(searchResults(maxNumLocalObsSearch))
876 call kdtree2_r_nearest(tp=tree, qv=refPosition, r2=maxRadius, nfound=numLocalObsFoundSearch,&
877 nalloc=maxNumLocalObsSearch, results=searchResults)
878 if (numLocalObsFoundSearch > maxNumLocalObsSearch) then
879 call utl_abort('eob_getLocalBodyIndices: the parameter maxNumLocalObsSearch must be increased')
880 end if
881
882 if ( vLocalize > 0.0d0 .and. vertLocation /= MPC_missingValue_R8 ) then
883 ! copy search results to output vectors, only those within vertical localization distance
884 numLocalObsFound = 0
885 numLocalObs = 0
886 do localObsIndex=1, numLocalObsFoundSearch
887 distance = abs( vertLocation - ensObs%vertLocation(searchResults(localObsIndex)%idx) )
888 if (distance <= vLocalize .and. ensObs%assFlag(searchResults(localObsIndex)%idx)==1) then
889 numLocalObsFound = numLocalObsFound + 1
890 if (numLocalObs < maxNumLocalObs) then
891 numLocalObs = numLocalObs + 1
892 localBodyIndices(numLocalObs) = searchResults(localObsIndex)%idx
893 distances(numLocalObs) = sqrt(searchResults(localObsIndex)%dis)
894 end if
895 end if
896 end do
897 else
898 ! no vertical location, so just copy results
899 numLocalObsFound = 0
900 numLocalObs = 0
901 do localObsIndex=1, numLocalObsFoundSearch
902 if (ensObs%assFlag(searchResults(localObsIndex)%idx)==1) then
903 numLocalObsFound = numLocalObsFound + 1
904 if (numLocalObs < maxNumLocalObs) then
905 numLocalObs = numLocalObs + 1
906 localBodyIndices(numLocalObs) = searchResults(localObsIndex)%idx
907 distances(numLocalObs) = sqrt(searchResults(localObsIndex)%dis)
908 end if
909 end if
910 end do
911 end if
912 deallocate(searchResults)
913
914 end function eob_getLocalBodyIndices
915
916 !--------------------------------------------------------------------------
917 ! eob_setLatLonObs
918 !--------------------------------------------------------------------------
919 subroutine eob_setLatLonObs(ensObs)
920 implicit none
921
922 ! Arguments:
923 type(struct_eob), intent(inout) :: ensObs
924
925 call obs_extractObsRealHeaderColumn(ensObs%lat, ensObs%obsSpaceData, OBS_LAT)
926 call obs_extractObsRealHeaderColumn(ensObs%lon, ensObs%obsSpaceData, OBS_LON)
927 call obs_extractObsRealBodyColumn(ensObs%obsValue, ensObs%obsSpaceData, OBS_VAR)
928
929 end subroutine eob_setLatLonObs
930
931 !--------------------------------------------------------------------------
932 ! eob_setobsErrInv
933 !--------------------------------------------------------------------------
934 subroutine eob_setobsErrInv(ensObs)
935 implicit none
936
937 ! Arguments:
938 type(struct_eob), intent(inout) :: ensObs
939
940 ! Locals:
941 integer :: obsIndex
942
943 call obs_extractObsRealBodyColumn(ensObs%obsErrInv, ensObs%obsSpaceData, OBS_OER)
944 do obsIndex = 1, ensObs%numObs
945 if(ensObs%obsErrInv(obsIndex) > 0.0d0) then
946 ensObs%obsErrInv(obsIndex) = 1.0d0/(ensObs%obsErrInv(obsIndex)**2)
947 else
948 ensObs%obsErrInv(obsIndex) = 0.0d0
949 end if
950 end do
951
952 ! read namelist if necessary and calculate obs error inverse for
953 ! passive and simulated observations
954 call eob_init()
955 if (psvObsAssim) call eob_setPsvObsErrInv(ensObs)
956 if (eob_simObsAssim) call eob_setSimObsErrInv(ensObs)
957
958 end subroutine eob_setobsErrInv
959
960 !--------------------------------------------------------------------------
961 ! eob_setPsvObsErrInv
962 !--------------------------------------------------------------------------
963 subroutine eob_setPsvObsErrInv(ensObs)
964 !
965 !:Purpose: Updates the inverse of the observation error variance
966 ! for passive osbervations and stores this in ensObs%obsErrInv.
967 ! This is done assuming that ensObs%obsErrInv was already set.
968 !
969 implicit none
970
971 ! Arguments:
972 type(struct_eob), intent(inout) :: ensObs
973
974 ! Locals:
975 integer :: obsIndex, headerIndex
976 integer :: codtyp, varnum, obsfamIndex
977 character(2) :: obsfamCurrent
978 logical :: psvFlag
979
980 do obsIndex = 1, ensObs%numObs
981 psvFlag = .false.
982 headerIndex = obs_bodyElem_i(ensObs%obsSpaceData, OBS_HIND, obsIndex)
983 obsfamCurrent = obs_getFamily(ensObs%obsSpaceData, headerIndex_opt=headerIndex)
984 ! update obs error inverse to 0 if current observation is passive
985 if (ANY(psvObsFamily == obsfamCurrent)) then
986 obsfamIndex = utl_findloc(psvObsFamily(:), obsfamCurrent)
987 if ((numPsvCodTyp(obsfamIndex) > 0) .and. (numPsvVarNum(obsfamIndex) > 0)) then
988 ! at least 1 codtyp AND varnum specified for current obs family so
989 ! see if current observation matches any of those codtypes AND
990 ! any of those varnums
991 codtyp = obs_headElem_i(ensObs%obsSpaceData, OBS_ITY, headerIndex)
992 varnum = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, obsIndex)
993 if (ANY(psvCodTyp(obsfamIndex,:) == codtyp) .and. ANY(psvVarNum(obsfamIndex,:) == varnum)) then
994 ensObs%obsErrInv(obsIndex) = 0.0d0
995 psvFlag = .true.
996 end if
997 else if (numPsvCodTyp(obsfamIndex) > 0) then
998 ! at least 1 codtype is specified for current obs family so
999 ! see if current observation matches any of those codtypes
1000 codtyp = obs_headElem_i(ensObs%obsSpaceData, OBS_ITY, headerIndex)
1001 if (ANY(psvCodTyp(obsfamIndex,:) == codtyp)) then
1002 ensObs%obsErrInv(obsIndex) = 0.0d0
1003 psvFlag = .true.
1004 end if
1005 else if (numPsvVarNum(obsfamIndex) > 0) then
1006 ! at least 1 varnum is specified for current obs family so
1007 ! see if current observation matches any of those varnums
1008 varnum = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, obsIndex)
1009 if (ANY(psvVarNum(obsfamIndex,:) == varnum)) then
1010 ensObs%obsErrInv(obsIndex) = 0.0d0
1011 psvFlag = .true.
1012 end if
1013 else
1014 ! passive observation family doesn't include any codtype or
1015 ! any varnum, so set error inverse to 0 irrespective of current
1016 ! observation codtype and varnum
1017 ensObs%obsErrInv(obsIndex) = 0.0d0
1018 psvFlag = .true.
1019 end if
1020 ! set OBS_FLG to indicate passive observation
1021 if (psvFlag) call obs_bodySet_i(ensObs%obsSpaceData,OBS_FLG,obsIndex, &
1022 ibset(obs_bodyElem_i(ensObs%obsSpaceData,OBS_FLG,obsIndex),25))
1023 end if
1024 end do
1025
1026 end subroutine eob_setPsvObsErrInv
1027
1028 !--------------------------------------------------------------------------
1029 ! eob_setSimObsErrInv
1030 !--------------------------------------------------------------------------
1031 subroutine eob_setSimObsErrInv(ensObs)
1032 !
1033 !:Purpose: Computes the inverse of the observation error variance if
1034 ! simulating any observations. Stores this in
1035 ! ensObs%obsErrInv_sim.
1036 !
1037 implicit none
1038
1039 ! Arguments:
1040 type(struct_eob), intent(inout) :: ensObs
1041
1042 ! Locals:
1043 integer :: obsIndex, headerIndex
1044 integer :: codtyp, varnum, obsfamIndex
1045 character(2) :: obsfamCurrent
1046 logical :: simFlag
1047
1048 write(*,*) 'eob_setSimObsErrInv: starting'
1049
1050 ! set to copy of regular obs error inverse
1051 ensObs%obsErrInv_sim(:) = ensObs%obsErrInv(:)
1052
1053 ! loop through all observations, and update the obs error inverse
1054 ! to a value of 0 for simulated observations
1055 do obsIndex = 1, ensObs%numObs
1056 simFlag = .false.
1057 headerIndex = obs_bodyElem_i(ensObs%obsSpaceData, OBS_HIND, obsIndex)
1058 obsfamCurrent = obs_getFamily(ensObs%obsSpaceData, headerIndex_opt=headerIndex)
1059 ! update obs error inverse for mean update to 0 if current observation is simulated
1060 if (ANY(simObsFamily == obsfamCurrent)) then
1061 obsfamIndex = utl_findloc(simObsFamily(:), obsfamCurrent)
1062 if ((numSimCodTyp(obsfamIndex) > 0) .and. (numSimVarNum(obsfamIndex) > 0)) then
1063 ! at least 1 codtyp AND varnum specified for current obs family so
1064 ! see if current observation matches any of those codtypes AND
1065 ! any of those varnums
1066 codtyp = obs_headElem_i(ensObs%obsSpaceData, OBS_ITY, headerIndex)
1067 varnum = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, obsIndex)
1068 if (ANY(simCodTyp(obsfamIndex,:) == codtyp) .and. ANY(simVarNum(obsfamIndex,:) == varnum)) then
1069 ensObs%obsErrInv_sim(obsIndex) = 0.0d0
1070 simFlag = .true.
1071 end if
1072 else if (numSimCodTyp(obsfamIndex) > 0) then
1073 ! at least 1 codtype is specified for current obs family so
1074 ! see if current observation matches any of those codtypes
1075 codtyp = obs_headElem_i(ensObs%obsSpaceData, OBS_ITY, headerIndex)
1076 if (ANY(simCodTyp(obsfamIndex,:) == codtyp)) then
1077 ensObs%obsErrInv_sim(obsIndex) = 0.0d0
1078 simFlag = .true.
1079 end if
1080 else if (numSimVarNum(obsfamIndex) > 0) then
1081 ! at least 1 varnum is specified for current obs family so
1082 ! see if current observation matches any of those varnums
1083 varnum = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, obsIndex)
1084 if (ANY(simVarNum(obsfamIndex,:) == varnum)) then
1085 ensObs%obsErrInv(obsIndex) = 0.0d0
1086 simFlag = .true.
1087 end if
1088 else
1089 ! simulated observation family doesn't include any codtype or
1090 ! any varnum, so set error inverse to 0 irrespective of current
1091 ! observation codtype and varnum
1092 ensObs%obsErrInv_sim(obsIndex) = 0.0d0
1093 simFlag = .true.
1094 end if
1095 ! set OBS_FLG to indicate simulated observation
1096 if (simFlag) call obs_bodySet_i(ensObs%obsSpaceData,OBS_FLG,obsIndex, &
1097 ibset(obs_bodyElem_i(ensObs%obsSpaceData,OBS_FLG,obsIndex),24))
1098 end if
1099 end do
1100
1101 end subroutine eob_setSimObsErrInv
1102
1103 !--------------------------------------------------------------------------
1104 ! eob_setVertLocation
1105 !--------------------------------------------------------------------------
1106 subroutine eob_setVertLocation(ensObs, columnMeanTrl)
1107 !
1108 !:Purpose: Set the vertical location value for each observation that
1109 ! will be used when doing vertical localization. For
1110 ! radiance observations, the level of the maximum value
1111 ! of the derivative of transmission is used. This value
1112 ! is also written in obsSpaceData in the OBS_ZHA column.
1113 !
1114 implicit none
1115
1116 ! Arguments:
1117 type(struct_eob) , intent(inout) :: ensObs
1118 type(struct_columnData), intent(in) :: columnMeanTrl
1119
1120 ! Locals:
1121 integer :: obsIndex, headerIndex, channelIndex, tovsIndex, numTovsLevels, nosensor
1122 integer :: levIndex, levIndexBelow, levIndexAbove, nLev_M
1123 integer :: varNumber(ensObs%numObs), obsVcoCode(ensObs%numObs), codType(ensObs%numObs)
1124 real(8) :: obsHeight, interpFactor, obsPPP(ensObs%numObs)
1125 real(8), pointer :: sfcPres_ptr(:,:), presM_ptr(:,:), heightM_ptr(:,:)
1126 type(rttov_profile), pointer :: profiles(:)
1127 logical :: verbose = .false.
1128
1129 call eob_setAssFlag(ensObs)
1130
1131 call obs_extractObsRealBodyColumn(obsPPP, ensObs%obsSpaceData, OBS_PPP)
1132 call obs_extractObsIntBodyColumn(varNumber, ensObs%obsSpaceData, OBS_VNM)
1133 call obs_extractObsIntBodyColumn(obsVcoCode, ensObs%obsSpaceData, OBS_VCO)
1134 call obs_extractObsIntHeaderColumn(codType, ensObs%obsSpaceData, OBS_ITY)
1135
1136 if (ensObs%typeVertCoord == 'logPressure') then
1137
1138 presM_ptr => col_getAllColumns(columnMeanTrl,'P_M')
1139 heightM_ptr => col_getAllColumns(columnMeanTrl,'Z_M')
1140 sfcPres_ptr => col_getAllColumns(columnMeanTrl,'P0')
1141 nLev_M = col_getNumLev(columnMeanTrl,'MM')
1142
1143 call tvs_getProfile(profiles,'nl')
1144
1145 end if
1146
1147 OBS_LOOP: do obsIndex = 1, ensObs%numObs
1148 headerIndex = obs_bodyElem_i(ensObs%obsSpaceData,OBS_HIND,obsIndex)
1149
1150 if( varNumber(obsIndex) == BUFR_NETS .or. varNumber(obsIndex) == BUFR_NEPS .or. &
1151 varNumber(obsIndex) == BUFR_NEUS .or. varNumber(obsIndex) == BUFR_NEVS .or. &
1152 varNumber(obsIndex) == BUFR_NESS .or. varNumber(obsIndex) == BUFR_NEPN .or. &
1153 varNumber(obsIndex) == BUFR_VIS .or. varNumber(obsIndex) == BUFR_LOGVIS .or. &
1154 varNumber(obsIndex) == BUFR_GUST .or. &
1155 varNumber(obsIndex) == BUFR_radarPrecip .or. varNumber(obsIndex) == BUFR_logRadarPrecip ) then
1156
1157 ! all surface observations
1158 if (ensObs%typeVertCoord == 'logPressure') then
1159 ensObs%vertLocation(obsIndex) = log(sfcPres_ptr(1,headerIndex))
1160 else if (ensObs%typeVertCoord == 'depth') then
1161 ensObs%vertLocation(obsIndex) = 0.0D0
1162 else
1163 call utl_abort('eob_setVertLocation: unknown typeVertCoord:' // trim(ensObs%typeVertCoord))
1164 end if
1165
1166 else if (varNumber(obsIndex) == BUFR_NEZD) then
1167
1168 ! ZTD observation, try 0.7*Psfc (i.e. ~700hPa when Psfc=1000hPa)
1169 if (ensObs%typeVertCoord == 'logPressure') then
1170 ensObs%vertLocation(obsIndex) = log(0.7D0 * sfcPres_ptr(1,headerIndex))
1171 else
1172 call utl_abort('eob_setVertLocation: ZTD obs only compatible with logPressure coordinate')
1173 end if
1174
1175 else if (obsPPP(obsIndex) > 0.0d0 .and. obsVcoCode(obsIndex)==2) then
1176
1177 ! all pressure level observations
1178 if (ensObs%typeVertCoord == 'logPressure') then
1179 ensObs%vertLocation(obsIndex) = log(obsPPP(obsIndex))
1180 else
1181 call utl_abort('eob_setVertLocation: pressure obs only compatible with logPressure coordinate')
1182 end if
1183
1184 else if(obsVcoCode(obsIndex) == 1 .and. .not.bufr_isOceanObs(varNumber(obsIndex))) then
1185
1186 if (ensObs%typeVertCoord /= 'logPressure') then
1187 ! skip this obs if it will not be assimilated
1188 if (ensObs%assFlag(obsIndex) == 0) cycle OBS_LOOP
1189 ! otherwise, abort
1190 write(*,*) 'eob_setVertLocation: varNum = ', varNumber(obsIndex)
1191 call utl_abort('eob_setVertLocation: height level obs only compatible with logPressure coordinate')
1192 end if
1193
1194 ! all height level observations (not including surface obs)
1195 obsHeight = obsPPP(obsIndex)
1196
1197 ! find level just below the observation
1198 levIndexBelow = 0
1199 LEV_LOOP: do levIndex = 1, nLev_M
1200 if (obsHeight > heightM_ptr(levIndex,headerIndex)) then
1201 levIndexBelow = levIndex
1202 exit LEV_LOOP
1203 end if
1204 end do LEV_LOOP
1205
1206 ! set the log pressure for observation
1207 if (levIndexBelow == 1) then
1208 ! above top level, use top level pressure
1209 ensObs%vertLocation(obsIndex) = log(presM_ptr(1,headerIndex))
1210 else if (levIndexBelow == 0) then
1211 ! below bottom level, use surface pressure
1212 ensObs%vertLocation(obsIndex) = log(sfcPres_ptr(1,headerIndex))
1213 else
1214 ! interpolate
1215 levIndexAbove = levIndexBelow - 1
1216 interpFactor = ( obsHeight - heightM_ptr(levIndexBelow,headerIndex) ) / &
1217 ( heightM_ptr(levIndexAbove,headerIndex) - heightM_ptr(levIndexBelow,headerIndex) )
1218 ensObs%vertLocation(obsIndex) = interpFactor * log(presM_ptr(levIndexAbove,headerIndex)) + &
1219 (1.0D0 - interpFactor) * log(presM_ptr(levIndexBelow,headerIndex))
1220 end if
1221
1222 else if(tvs_isIdBurpTovs(codType(obsIndex))) then
1223
1224 if (ensObs%typeVertCoord /= 'logPressure') then
1225 call utl_abort('eob_setVertLocation: radiance obs only compatible with logPressure coordinate')
1226 end if
1227
1228 tovsIndex = tvs_tovsIndex(headerIndex)
1229 nosensor = tvs_lsensor(tovsIndex)
1230 numTovsLevels = size(tvs_transmission(tovsIndex)%tau_levels,1)
1231 channelIndex = nint(obsPPP(obsIndex))
1232 channelIndex = max(0,min(channelIndex,tvs_maxChannelNumber+1))
1233 channelIndex = channelIndex - tvs_channelOffset(nosensor)
1234 channelIndex = utl_findloc(tvs_ichan(:,nosensor), channelIndex)
1235 if (channelIndex > 0 .and. ensObs%assFlag(obsIndex)==1) then
1236 call max_transmission(tvs_transmission(tovsIndex), numTovsLevels, &
1237 channelIndex, profiles(tovsIndex)%p, ensObs%vertLocation(obsIndex))
1238 if(mmpi_myid == 0 .and. verbose) then
1239 write(*,*) 'eob_setVertLocation for tovs: ', codType(obsIndex), &
1240 obsPPP(obsIndex), 0.01*exp(ensObs%vertLocation(obsIndex))
1241 end if
1242 else
1243 ensObs%vertLocation(obsIndex) = log(500.0D2)
1244 end if
1245
1246 else if (varNumber(obsIndex) == BUFR_SST) then
1247
1248 if (ensObs%typeVertCoord /= 'depth') then
1249 call utl_abort('eob_setVertLocation: SST obs only compatible with ocean depth coordinate')
1250 end if
1251
1252 ! SST observations
1253 ensObs%vertLocation(obsIndex) = minval(columnMeanTrl%vco%depths(:))
1254
1255 else if(ensObs%assFlag(obsIndex)==1) then
1256
1257 write(*,*) 'eob_setLatLonPresObs: ERROR! cannot compute pressure for this observation: ', &
1258 obsPPP(obsIndex), varNumber(obsIndex), obsVcoCode(obsIndex)
1259 call utl_abort('eob_setVertLocation')
1260
1261 end if
1262
1263 ! write the value into obsSpaceData for later diagnostics
1264 if (ensObs%assFlag(obsIndex)==1) then
1265 call obs_bodySet_r(ensObs%obsSpaceData, OBS_ZHA, obsIndex, &
1266 ensObs%vertLocation(obsIndex))
1267 end if
1268
1269 end do OBS_LOOP
1270
1271 nullify(profiles)
1272
1273 end subroutine eob_setVertLocation
1274
1275 !--------------------------------------------------------------------------
1276 ! eob_setAssFlag
1277 !--------------------------------------------------------------------------
1278 subroutine eob_setAssFlag(ensObs)
1279 implicit none
1280
1281 ! Arguments:
1282 type(struct_eob), intent(inout) :: ensObs
1283
1284 call obs_extractObsIntBodyColumn(ensObs%assFlag, ensObs%obsSpaceData, OBS_ASS)
1285
1286 end subroutine eob_setAssFlag
1287
1288 !--------------------------------------------------------------------------
1289 ! eob_setYb
1290 !--------------------------------------------------------------------------
1291 subroutine eob_setYb(ensObs, memberIndex)
1292 implicit none
1293
1294 ! Arguments:
1295 type(struct_eob), intent(inout) :: ensObs
1296 integer , intent(in) :: memberIndex
1297
1298 ! get the Y-HX value from obsSpaceData
1299 call obs_extractObsRealBodyColumn_r4(ensObs%Yb_r4(memberIndex,:), ensObs%obsSpaceData, OBS_OMP)
1300
1301 ! now compute HX = Y - (Y-HX)
1302 ensObs%Yb_r4(memberIndex,:) = ensObs%obsValue(:) - ensObs%Yb_r4(memberIndex,:)
1303
1304 end subroutine eob_setYb
1305
1306 !--------------------------------------------------------------------------
1307 ! eob_setYa (like eob_setYb but for the analysis)
1308 !--------------------------------------------------------------------------
1309 subroutine eob_setYa(ensObs, memberIndex, obsColumnName)
1310 implicit none
1311
1312 ! Arguments:
1313 type(struct_eob), intent(inout) :: ensObs
1314 integer , intent(in) :: memberIndex
1315 integer , intent(in) :: obsColumnName
1316
1317 if ( .not. allocated(ensObs%Ya_r4) ) then
1318 call utl_abort('eob_setYa: ensObs%Ya_r4 must be allocated and it is not')
1319 end if
1320
1321 ! get the Y-HX value from obsSpaceData
1322 call obs_extractObsRealBodyColumn_r4(ensObs%Ya_r4(memberIndex,:), ensObs%obsSpaceData, obsColumnName)
1323
1324 ! now compute HX = Y - (Y-HX)
1325 ensObs%Ya_r4(memberIndex,:) = ensObs%obsValue(:) - ensObs%Ya_r4(memberIndex,:)
1326
1327 end subroutine eob_setYa
1328
1329 !--------------------------------------------------------------------------
1330 ! eob_setSimObsVal
1331 !--------------------------------------------------------------------------
1332 subroutine eob_setSimObsVal(ensObs)
1333 !
1334 !:Purpose: Set the observed value for simulated observations to
1335 ! the background ensemble mean in observation space.
1336 !
1337 implicit none
1338
1339 ! Arguments:
1340 type(struct_eob) , intent(inout) :: ensObs
1341
1342 ! Locals:
1343 integer :: obsIndex, headerIndex
1344 integer :: codtyp, varnum, obsfamIndex
1345 character(2) :: obsfamCurrent
1346
1347 if (.not. eob_simObsAssim) return
1348 ! Loop through observations and set y to mean(H(x)) if y is in obs family of interest
1349 do obsIndex = 1, ensObs%numObs
1350 headerIndex = obs_bodyElem_i(ensObs%obsSpaceData, OBS_HIND, obsIndex)
1351 obsfamCurrent = obs_getFamily(ensObs%obsSpaceData, headerIndex_opt=headerIndex)
1352 if (ANY(simObsFamily == obsfamCurrent)) then
1353 obsfamIndex = utl_findloc(simObsFamily(:), obsfamCurrent)
1354 if ((numSimCodTyp(obsfamIndex) > 0) .and. (numSimVarNum(obsfamIndex) > 0)) then
1355 ! at least 1 codtyp AND varnum specified for current obs family so
1356 ! see if current observation matches any of those codtypes AND
1357 ! any of those varnums
1358 codtyp = obs_headElem_i(ensObs%obsSpaceData, OBS_ITY, headerIndex)
1359 varnum = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, obsIndex)
1360 if (ANY(simCodTyp(obsfamIndex,:) == codtyp) .and. ANY(simVarNum(obsfamIndex,:) == varnum)) then
1361 ensObs%obsvalue(obsIndex) = ensObs%meanYb(obsIndex)
1362 end if
1363 else if (numSimCodTyp(obsfamIndex) > 0) then
1364 ! at least 1 codtype is specified for current obs family so
1365 ! see if current observation matches any of those codtypes
1366 codtyp = obs_headElem_i(ensObs%obsSpaceData, OBS_ITY, headerIndex)
1367 if (ANY(simCodTyp(obsfamIndex,:) == codtyp)) then
1368 ensObs%obsvalue(obsIndex) = ensObs%meanYb(obsIndex)
1369 end if
1370 else if (numSimVarNum(obsfamIndex) > 0) then
1371 ! at least 1 varnum is specified for current obs family so
1372 ! see if current observation matches any of those varnums
1373 varnum = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, obsIndex)
1374 if (ANY(simVarNum(obsfamIndex,:) == varnum)) then
1375 ensObs%obsvalue(obsIndex) = ensObs%meanYb(obsIndex)
1376 end if
1377 else
1378 ! simulated observation family doesn't include any codtype
1379 ! so set irrespective of current observation's codtype
1380 ensObs%obsvalue(obsIndex) = ensObs%meanYb(obsIndex)
1381 end if
1382 end if
1383 end do
1384
1385 end subroutine eob_setSimObsVal
1386
1387 !--------------------------------------------------------------------------
1388 ! eob_setDeterYb
1389 !--------------------------------------------------------------------------
1390 subroutine eob_setDeterYb(ensObs)
1391 implicit none
1392
1393 ! Arguments:
1394 type(struct_eob), intent(inout) :: ensObs
1395
1396 ! get the Y-HX value from obsSpaceData
1397 call obs_extractObsRealBodyColumn(ensObs%DeterYb(:), ensObs%obsSpaceData, OBS_OMP)
1398
1399 ! now compute HX = Y - (Y-HX)
1400 ensObs%DeterYb(:) = ensObs%obsValue(:) - ensObs%DeterYb(:)
1401
1402 end subroutine eob_setDeterYb
1403
1404 !--------------------------------------------------------------------------
1405 ! eob_calcAndRemoveMeanYb
1406 !--------------------------------------------------------------------------
1407 subroutine eob_calcAndRemoveMeanYb(ensObs)
1408 implicit none
1409
1410 ! Arguments:
1411 type(struct_eob), intent(inout) :: ensObs
1412
1413 ! Locals:
1414 integer :: obsIndex
1415
1416 do obsIndex = 1, ensObs%numObs
1417 ensObs%meanYb(obsIndex) = sum(ensObs%Yb_r4(:,obsIndex)) / ensObs%numMembers
1418 ensObs%Yb_r4(:,obsIndex) = ensObs%Yb_r4(:,obsIndex) - ensObs%meanYb(obsIndex)
1419 end do
1420
1421 ensObs%meanRemoved = .true.
1422
1423 end subroutine eob_calcAndRemoveMeanYb
1424
1425 !--------------------------------------------------------------------------
1426 ! eob_calcRandPert
1427 !--------------------------------------------------------------------------
1428 subroutine eob_calcRandPert(ensObs, randomSeed)
1429 implicit none
1430
1431 ! Arguments:
1432 type(struct_eob), intent(inout) :: ensObs
1433 integer , intent(in) :: randomSeed
1434
1435 ! Locals:
1436 integer :: obsIndex, memberIndex
1437 real(4) :: meanRandPert, sigObs
1438
1439 call rng_setup(abs(randomSeed))
1440 if ( allocated(ensObs%randPert_r4) ) then
1441 call utl_abort('eob_calcRandPert: ensObs%randPert_r4 must not be already allocated')
1442 end if
1443
1444 allocate(ensObs%randPert_r4(ensObs%numMembers,ensObs%numObs))
1445
1446 do obsIndex = 1, ensObs%numObs
1447 sigObs = obs_bodyElem_r(ensObs%obsSpaceData, OBS_OER, obsIndex)
1448 do memberIndex = 1, ensObs%numMembers
1449 ensObs%randPert_r4(memberIndex,obsIndex) = sigObs * rng_gaussian()
1450 end do
1451
1452 meanRandPert = sum(ensObs%randPert_r4(:,obsIndex)) / real(ensObs%numMembers,4)
1453 ensObs%randPert_r4(:,obsIndex) = ensObs%randPert_r4(:,obsIndex) - meanRandPert
1454 end do
1455
1456 end subroutine eob_calcRandPert
1457
1458 !--------------------------------------------------------------------------
1459 ! eob_setMeanOMP
1460 !--------------------------------------------------------------------------
1461 subroutine eob_setMeanOMP(ensObs)
1462 implicit none
1463
1464 ! Arguments:
1465 type(struct_eob), intent(inout) :: ensObs
1466
1467 ! Locals:
1468 integer :: obsIndex
1469
1470 do obsIndex = 1, ensObs%numObs
1471 if (obs_bodyElem_i(ensObs%obsSpaceData, OBS_ASS, obsIndex) == obs_notAssimilated) cycle
1472 call obs_bodySet_r(ensObs%obsSpaceData, OBS_OMP, obsIndex, &
1473 ensObs%obsValue(obsIndex)-ensObs%meanYb(obsIndex))
1474 end do
1475
1476 end subroutine eob_setMeanOMP
1477
1478 !--------------------------------------------------------------------------
1479 ! eob_setHPHT
1480 !--------------------------------------------------------------------------
1481 subroutine eob_setHPHT(ensObs)
1482 implicit none
1483
1484 ! Arguments:
1485 type(struct_eob), intent(inout) :: ensObs
1486
1487 ! Locals:
1488 integer :: obsIndex, memberIndex
1489 real(8) :: hpht
1490
1491 do obsIndex = 1, ensObs%numObs
1492 if (obs_bodyElem_i(ensObs%obsSpaceData, OBS_ASS, obsIndex) == obs_notAssimilated) cycle
1493 hpht = 0.0d0
1494 do memberIndex = 1, ensObs%numMembers
1495 hpht = hpht + ensObs%Yb_r4(memberIndex,obsIndex)**2 / ensObs%numMembers
1496 end do
1497 if (hpht > 0.0D0) then
1498 hpht = sqrt(hpht)
1499 else
1500 hpht = 0.0D0
1501 end if
1502 call obs_bodySet_r(ensObs%obsSpaceData, OBS_HPHT, obsIndex, hpht)
1503
1504 end do
1505
1506 end subroutine eob_setHPHT
1507
1508 !--------------------------------------------------------------------------
1509 ! eob_backgroundCheck
1510 !--------------------------------------------------------------------------
1511 subroutine eob_backgroundCheck(ensObs)
1512 !
1513 !:Purpose: Apply additional background using the ensemble spread.
1514 !
1515 implicit none
1516
1517 ! Arguments::
1518 type(struct_eob), intent(inout) :: ensObs
1519
1520 ! Locals:
1521 integer :: bodyIndexBeg, bodyIndexEnd, headerIndex, ivar, bodyIndex
1522 integer :: numRejected, numRejectedMpiGlobal, ierr, windCount
1523 real :: sigo, sigb, omp, sig, reject_limit
1524 logical :: reject_wind
1525
1526 numRejected = 0
1527 reject_limit = 5.0
1528 do headerIndex = 1, obs_numheader(ensObs%obsSpaceData)
1529 bodyIndexBeg = obs_headElem_i(ensObs%obsSpaceData, OBS_RLN, headerIndex)
1530 bodyIndexEnd = obs_headElem_i(ensObs%obsSpaceData, OBS_NLV, headerIndex) + bodyIndexBeg - 1
1531 reject_wind = .false.
1532
1533 do bodyIndex = bodyIndexBeg, bodyIndexEnd
1534 if (obs_bodyElem_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex) == obs_notAssimilated) cycle
1535 sigo = obs_bodyElem_r(ensObs%obsSpaceData, OBS_OER, bodyIndex)
1536 sigb = obs_bodyElem_r(ensObs%obsSpaceData, OBS_HPHT, bodyIndex)
1537 ! cut off at reject_limit standard deviations
1538 sig = reject_limit*(sigo**2 + sigb**2)**0.5
1539 omp = abs(obs_bodyElem_r(ensObs%obsSpaceData, OBS_OMP, bodyIndex))
1540 if (omp > sig) then
1541 call obs_bodySet_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
1542 call obs_bodySet_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex, &
1543 IBSET(obs_bodyElem_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex),9))
1544 ivar = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, bodyIndex)
1545 reject_wind = bufr_isWindComponent(ivar)
1546 write(*,*) 'eob_backgroundCheck: headerIndex, omp, sig, ivar, reject_wind', &
1547 headerIndex, omp, sig, ivar, reject_wind
1548 numRejected = numRejected + 1
1549 end if
1550 end do
1551
1552 ! take the same reject decision for both components of the
1553 ! wind vector (and any other winds for that station).
1554 ! N.B.: This seems to assume only one level per header, this is generally
1555 ! the case for wind observations currently, since radiosondes are 4D
1556 if (reject_wind) then
1557 ! first count how many wind observations we have for this station
1558 windCount = 0
1559 do bodyIndex = bodyIndexBeg, bodyIndexEnd
1560 if (bufr_isWindComponent(obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, bodyIndex))) then
1561 windCount = windCount + 1
1562 end if
1563 end do
1564 if (windCount > 2) then
1565 write(*,*) 'eob_backgroundCheck: WARNING'
1566 write(*,*) 'Station ',headerIndex,' has ',windCount,' wind observations '
1567 write(*,*) 'Perhaps old radiosonde format - std dev not changed for other wind component'
1568 else
1569 do bodyIndex = bodyIndexBeg, bodyIndexEnd
1570 if (obs_bodyElem_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex) == obs_notassimilated) cycle
1571 ivar = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, bodyIndex)
1572 if (bufr_isWindComponent(ivar) .and. &
1573 obs_bodyElem_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) then
1574 call obs_bodySet_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
1575 call obs_bodySet_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex, &
1576 IBSET(obs_bodyElem_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex),9))
1577 write(*,*) 'eob_backgroundCheck: other wind component headerIndex, ivar', &
1578 headerIndex, ivar
1579 numRejected = numRejected + 1
1580 end if
1581 end do
1582 end if
1583 end if
1584 end do
1585
1586 call rpn_comm_allreduce(numRejected,numRejectedMpiGlobal,1,'mpi_integer','mpi_sum','GRID',ierr)
1587 write(*,*)
1588 write(*,*) 'eob_backgroundCheck: number of observations rejected (local) =', numRejected
1589 write(*,*) 'eob_backgroundCheck: number of observations rejected (global)=', numRejectedMpiGlobal
1590
1591 end subroutine eob_backgroundCheck
1592
1593 !--------------------------------------------------------------------------
1594 ! eob_removeObsNearLand
1595 !--------------------------------------------------------------------------
1596 subroutine eob_removeObsNearLand(ensObs, oceanMask, minDistanceToLand)
1597 !
1598 !:Purpose: Reject observations that are close to land as determined by
1599 ! the argument "minDistanceToLand". In the case of a depth-
1600 ! varying land mask, the first level (should be the surface) is
1601 ! used.
1602 !
1603 implicit none
1604
1605 ! Arguments:
1606 type(struct_eob), intent(inout) :: ensObs
1607 type(struct_ocm), intent(in) :: oceanMask
1608 real(8), intent(in) :: minDistanceToLand
1609
1610 ! Locals:
1611 integer :: headerIndex, bodyIndex, bodyIndexBeg, bodyIndexEnd, levIndex
1612 integer :: numRejected, numRejectedMpiGlobal, ierr
1613 real(8) :: obsLon, obsLat
1614
1615 write(*,*) 'eob_removeObsNearLand: starting'
1616
1617 numRejected = 0
1618
1619 HEADER_LOOP: do headerIndex = 1, obs_numheader(ensObs%obsSpaceData)
1620 bodyIndexBeg = obs_headElem_i(ensObs%obsSpaceData, OBS_RLN, headerIndex)
1621 bodyIndexEnd = obs_headElem_i(ensObs%obsSpaceData, OBS_NLV, headerIndex) + bodyIndexBeg - 1
1622
1623 levIndex = 1
1624 obsLat = obs_headElem_r(ensObs%obsSpaceData, OBS_LAT, headerIndex)
1625 obsLon = obs_headElem_r(ensObs%obsSpaceData, OBS_LON, headerIndex)
1626
1627 ! skip this obs if it is far from land
1628 if (ocm_farFromLand(oceanMask, levIndex, obsLon, obsLat, minDistanceToLand)) then
1629 cycle HEADER_LOOP
1630 end if
1631
1632 ! otherwise it is rejected
1633 BODY_LOOP: do bodyIndex = bodyIndexBeg, bodyIndexEnd
1634 if (obs_bodyElem_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex) == obs_notAssimilated) cycle BODY_LOOP
1635
1636 call obs_bodySet_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
1637 call obs_bodySet_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex, &
1638 IBSET(obs_bodyElem_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex),9))
1639 numRejected = numRejected + 1
1640 end do BODY_LOOP
1641 end do HEADER_LOOP
1642
1643 call rpn_comm_allreduce(numRejected,numRejectedMpiGlobal,1,'mpi_integer','mpi_sum','GRID',ierr)
1644 write(*,*)
1645 write(*,*) 'eob_removeObsNearLand: number of observations rejected (local) =', numRejected
1646 write(*,*) 'eob_removeObsNearLand: number of observations rejected (global)=', numRejectedMpiGlobal
1647
1648 write(*,*) 'eob_removeObsNearLand: finished'
1649
1650 end subroutine eob_removeObsNearLand
1651
1652 !--------------------------------------------------------------------------
1653 ! eob_setSigiSigo
1654 !--------------------------------------------------------------------------
1655 subroutine eob_setSigiSigo(ensObs)
1656 !
1657 !:Purpose: Apply huber norm quality control procedure. This modifies
1658 ! the OBS_OER value, but before that its value is copied into
1659 ! OBS_SIGO and also OBS_SIGI computed
1660 !
1661 implicit none
1662
1663 ! Arguments:
1664 type(struct_eob), intent(inout) :: ensObs
1665
1666 ! Locals:
1667 integer :: bodyIndex
1668 real(pre_obsReal) :: sigo, sigb, sigi
1669
1670 ! Set 'sigi' and 'sigo' before oer is modified by Huber norm
1671 do bodyIndex = 1, obs_numbody(ensObs%obsSpaceData)
1672 sigb = obs_bodyElem_r(ensObs%obsSpaceData, OBS_HPHT, bodyIndex)
1673 sigo = obs_bodyElem_r(ensObs%obsSpaceData, OBS_OER, bodyIndex)
1674 if (obs_bodyElem_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex) == obs_assimilated) then
1675 sigi = (sigo**2 + sigb**2)**0.5
1676 call obs_bodySet_r(ensObs%obsSpaceData, OBS_SIGI, bodyIndex, sigi)
1677 call obs_bodySet_r(ensObs%obsSpaceData, OBS_SIGO, bodyIndex, sigo)
1678 end if
1679 end do
1680
1681 end subroutine eob_setSigiSigo
1682
1683 !--------------------------------------------------------------------------
1684 ! eob_huberNorm
1685 !--------------------------------------------------------------------------
1686 subroutine eob_huberNorm(ensObs)
1687 !
1688 !:Purpose: Apply huber norm quality control procedure. This modifies
1689 ! the OBS_OER value.
1690 !
1691 implicit none
1692
1693 ! Arguments:
1694 type(struct_eob), intent(inout) :: ensObs
1695
1696 ! Locals:
1697 integer :: huberCount, huberCountMpiGlobal, ivar, windCount, ierr
1698 integer :: bodyIndex, bodyIndexBeg, bodyIndexEnd, headerIndex
1699 real(pre_obsReal) :: c_limit, sig, sigo, sigb, omp, sigo_hub, sigo_hub_wind
1700 logical :: reject_wind
1701
1702 c_limit = 2.0
1703 huberCount = 0
1704 do headerIndex = 1, obs_numheader(ensObs%obsSpaceData)
1705 bodyIndexBeg = obs_headElem_i(ensObs%obsSpaceData, OBS_RLN, headerIndex)
1706 bodyIndexEnd = obs_headElem_i(ensObs%obsSpaceData, OBS_NLV, headerIndex) + bodyIndexBeg - 1
1707 reject_wind = .false.
1708 sigo_hub_wind = 0.0
1709 do bodyIndex = bodyIndexBeg, bodyIndexEnd
1710 sigo = obs_bodyElem_r(ensObs%obsSpaceData, OBS_OER, bodyIndex)
1711 sigb = obs_bodyElem_r(ensObs%obsSpaceData, OBS_HPHT, bodyIndex)
1712 ! cut off at reject_limit standard deviations
1713 sig = c_limit*(sigo**2 + sigb**2)**0.5
1714 omp = abs(obs_bodyElem_r(ensObs%obsSpaceData, OBS_OMP, bodyIndex))
1715 if (omp > sig) then
1716 ! redefining the observational error such that the innovation
1717 ! is at exactly c_limit standard deviations.
1718 huberCount = huberCount + 1
1719 sigo_hub = ((omp/c_limit)**2-sigb**2)**0.5
1720 ivar = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, bodyIndex)
1721 if (bufr_isWindComponent(ivar)) then
1722 ! the other wind components will be changed at the same time (next if)
1723 reject_wind = .true.
1724 call obs_bodySet_r(ensObs%obsSpaceData, OBS_OER, bodyIndex, sigo_hub)
1725 ! this is for the special case both components of the wind innovation
1726 ! suggest using the Huber norm.
1727 if (sigo_hub > sigo_hub_wind) then
1728 sigo_hub_wind = sigo_hub
1729 end if
1730 else
1731 call obs_bodySet_r(ensObs%obsSpaceData, OBS_OER, bodyIndex, sigo_hub)
1732 end if
1733 end if
1734 end do
1735
1736 ! Give the same inflated observation error to both wind components.
1737 ! this part assumes that modern sondes are used (at most two
1738 ! wind observations per station.
1739 if (reject_wind) then
1740 ! first count how many wind observations we have for this station
1741 windCount = 0
1742 do bodyIndex = bodyIndexBeg, bodyIndexEnd
1743 if (bufr_isWindComponent(obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, bodyIndex))) then
1744 windCount = windCount + 1
1745 end if
1746 end do
1747 if (windCount > 2) then
1748 write(*,*) 'Warning Hubernorm'
1749 write(*,*) 'Station ',headerIndex,' has ',windCount,' wind observations '
1750 write(*,*) 'Perhaps old radiosonde format - std dev not changed for other wind component'
1751 else
1752 huberCount = huberCount + 1
1753 do bodyIndex = bodyIndexBeg, bodyIndexEnd
1754 if (bufr_isWindComponent(obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, bodyIndex))) then
1755 call obs_bodySet_r(ensObs%obsSpaceData, OBS_OER, bodyIndex, sigo_hub_wind)
1756 end if
1757 end do
1758 end if
1759 end if
1760 end do
1761
1762 call rpn_comm_allreduce(huberCount, huberCountMpiGlobal, 1, 'mpi_integer', 'mpi_sum', 'GRID', ierr)
1763 write(*,*)
1764 write(*,*) 'eob_huberNorm: number of obs with increased error stddev (local) = ', huberCount
1765 write(*,*) 'eob_huberNorm: number of obs with increased error stddev (global)= ', huberCountMpiGlobal
1766
1767 end subroutine eob_huberNorm
1768
1769 !--------------------------------------------------------------------------
1770 ! eob_rejectRadNearSfc
1771 !--------------------------------------------------------------------------
1772 subroutine eob_rejectRadNearSfc(ensObs)
1773 !
1774 !:Purpose: Reject all radiance observations with peak sensitivity
1775 ! too close to the surface.
1776 !
1777 implicit none
1778
1779 ! Arguments:
1780 type(struct_eob), intent(inout) :: ensObs
1781
1782 ! Locals:
1783 integer :: acceptCount, rejectCount, acceptCountMpiGlobal, rejectCountMpiGlobal
1784 integer :: varNumber, bodyIndex, ierr
1785 ! reject lower than 975 hPa
1786 real, parameter :: logPresRadianceLimit = 11.4876E0
1787
1788 acceptCount = 0
1789 rejectCount = 0
1790 do bodyIndex = 1, obs_numbody(ensObs%obsSpaceData)
1791 varNumber = obs_bodyElem_i(ensObs%obsSpaceData, OBS_VNM, bodyIndex)
1792 if (varNumber == bufr_nbt1 .or. &
1793 varNumber == bufr_nbt2 .or. &
1794 varNumber == bufr_nbt3) then
1795 if (ensObs%vertLocation(bodyIndex) < logPresRadianceLimit) then
1796 acceptCount = acceptCount + 1
1797 else
1798 rejectCount = rejectCount + 1
1799 call obs_bodySet_i(ensObs%obsSpaceData, OBS_ASS, bodyIndex, obs_notAssimilated)
1800 call obs_bodySet_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex, &
1801 IBSET(obs_bodyElem_i(ensObs%obsSpaceData, OBS_FLG, bodyIndex),9))
1802 end if
1803 end if
1804 end do
1805
1806 call rpn_comm_allreduce(acceptCount, acceptCountMpiGlobal, 1, 'mpi_integer', 'mpi_sum', 'GRID', ierr)
1807 call rpn_comm_allreduce(rejectCount, rejectCountMpiGlobal, 1, 'mpi_integer', 'mpi_sum', 'GRID', ierr)
1808 write(*,*)
1809 write(*,*) 'eob_rejectRadNearSfc: Number of accepted, rejected observations (local) : ', &
1810 acceptCount, rejectCount
1811 write(*,*) 'eob_rejectRadNearSfc: Number of accepted, rejected observations (global): ', &
1812 acceptCountMpiGlobal, rejectCountMpiGlobal
1813
1814 end subroutine eob_rejectRadNearSfc
1815
1816 !--------------------------------------------------------------------------
1817 ! getMemberIndexInFullEnsSet (private routine)
1818 !--------------------------------------------------------------------------
1819 subroutine getMemberIndexInFullEnsSet(ensObs, memberIndexArray, &
1820 numGroupsToDivideMembers_opt, &
1821 maxNumMembersPerGroup_opt)
1822 !
1823 !:Purpose: get memberIndex array corresponding to the full ensemble set. This
1824 ! is useful when ensObs is a subset of full ensemble members.
1825 ! If first member in ensObs is member 6, to get the full ensemble set equivalent of ensObs members:
1826 ! a) When members are not grouped (numGroupsToDivideMembers=1), all members are offset by
1827 ! memberIndexOffset (e.g. 6, 7, ..., 6+ensObs%numMermbers)
1828 ! b) When members are grouped (numGroupsToDivideMembers/=1), members within each group
1829 ! are offset by memberIndexOffset but there is increment of maxNumMembersPerGroup_opt
1830 ! to jump to the next group (e.g. if maxNumMembersPerGroup_opt=10, for first group 6, 7, 8, 9,
1831 ! for second group 6+10, 7+10, 8+10, 9+10, and so on)
1832 !
1833 implicit none
1834
1835 ! Arguments:
1836 type(struct_eob), intent(in) :: ensObs
1837 integer, intent(inout) :: memberIndexArray(:)
1838 integer, optional, intent(in) :: numGroupsToDivideMembers_opt
1839 integer, optional, intent(in) :: maxNumMembersPerGroup_opt
1840
1841 ! Locals:
1842 integer :: memberIndex, groupIndex, memberIndexOffset, memberIndexInGroup
1843 integer :: numGroupsToDivideMembers, numMembersPerGroup
1844
1845 if (present(numGroupsToDivideMembers_opt)) then
1846 numGroupsToDivideMembers = numGroupsToDivideMembers_opt
1847 else
1848 numGroupsToDivideMembers = 1
1849 end if
1850
1851 memberIndexOffset = ensObs%fileMemberIndex1
1852
1853 if (numGroupsToDivideMembers == 1) then
1854 do memberIndex = 1, ensObs%numMembers
1855 memberIndexArray(memberIndex) = memberIndex + memberIndexOffset - 1
1856 end do
1857 else
1858 if (.not. present(maxNumMembersPerGroup_opt)) then
1859 call utl_abort('getMemberIndexInFullEnsSet: maxNumMembersPerGroup_opt input argument missing')
1860 end if
1861
1862 ! divide members into groups
1863 numMembersPerGroup = ensObs%numMembers / numGroupsToDivideMembers
1864 if (numMembersPerGroup > maxNumMembersPerGroup_opt) then
1865 call utl_abort('getMemberIndexInFullEnsSet: numMembersPerGroup > maxNumMembersPerGroup_opt')
1866 end if
1867
1868 memberIndex = 0
1869 do groupIndex = 1, numGroupsToDivideMembers
1870 do memberIndexInGroup = 1, numMembersPerGroup
1871 memberIndex = memberIndex + 1
1872 memberIndexArray(memberIndex) = (groupIndex - 1) * maxNumMembersPerGroup_opt + &
1873 memberIndexInGroup + memberIndexOffset - 1
1874 end do
1875 end do
1876
1877 end if
1878
1879 end subroutine getMemberIndexInFullEnsSet
1880
1881 !--------------------------------------------------------------------------
1882 ! max_transmission (private routine)
1883 !--------------------------------------------------------------------------
1884 subroutine max_transmission(transmission, numLevels, transIndex, rttovPres, maxLnP)
1885 !
1886 !:Purpose: Determine the height in log pressure where we find the maximum
1887 ! value of the first derivative of transmission with respect to
1888 ! log pressure
1889 !
1890 implicit none
1891
1892 ! Arguments:
1893 type(rttov_transmission), intent(in) :: transmission ! transmission (rttov type)
1894 integer(kind=jpim) , intent(in) :: numLevels ! number of RTTOV levels
1895 integer , intent(in) :: transIndex ! index of transmission%tau_levels
1896 real(kind=jprb), pointer, intent(in) :: rttovPres(:) ! pressure of RTTOV levels
1897 real(8) , intent(out) :: maxLnP ! log pressure of maximum
1898
1899 ! Locals:
1900 integer :: levIndex
1901 real(8) :: lnPres(numLevels), avgPres(numLevels-1)
1902 real(8) :: diffTau, derivTau(numLevels), maxDeriv
1903 integer :: nAvgLev, maxIndex
1904
1905 nAvgLev = numLevels - 1
1906 lnPres(:) = log(rttovPres(:)*MPC_PA_PER_MBAR_R8)
1907 ! calculate the first derivative of transmission with respect to log pressure
1908 ! and find the level index for its maximum
1909 maxDeriv = -0.1d0
1910 derivTau(1) = 0.0d0
1911 maxIndex = numLevels
1912 do levIndex = 2, numLevels
1913 avgPres(levIndex-1) = 0.5d0*(lnPres(levIndex)+lnPres(levIndex-1))
1914 diffTau = transmission%tau_levels(levIndex-1,transIndex) - transmission%tau_levels(levIndex,transIndex)
1915 derivTau(levIndex) = diffTau / (lnPres(levIndex)-lnPres(levIndex-1))
1916 if (derivTau(levIndex)>maxDeriv) then
1917 maxDeriv = derivTau(levIndex)
1918 maxIndex = levIndex
1919 end if
1920 end do
1921
1922 ! get the height in log pressure for the level index (maxIndex) found above
1923 if (maxIndex==1) maxIndex = maxIndex + 1
1924 if ((maxIndex==2).or.(maxIndex==numLevels)) then
1925 maxLnP = avgPres(maxIndex-1)
1926 else
1927 call get_peak(maxIndex,nAvgLev,avgPres,derivTau,maxLnP)
1928 end if
1929
1930 end subroutine max_transmission
1931
1932 !--------------------------------------------------------------------------
1933 ! get_peak (private routine)
1934 !--------------------------------------------------------------------------
1935 subroutine get_peak(maxIndex,nlev,lnp,deriv,maxLnP)
1936 !
1937 !:Purpose: Do quadratic interpolation to find pressure of peak transmission.
1938 !
1939 implicit none
1940
1941 ! Arguments:
1942 integer, intent(in) :: maxIndex, nlev
1943 real(8), intent(in) :: lnp(nlev), deriv(nlev+1)
1944 real(8), intent(inout) :: maxLnP
1945
1946 ! Locals:
1947 external :: dgesv
1948 integer, parameter :: N=3
1949 integer :: info
1950 integer, parameter :: lda=N, ldb=N, nrhs=1
1951 integer :: ipiv(N)
1952 real(8) :: A(lda,N),B(ldb,nrhs)
1953 integer :: index1, index2
1954
1955 index2 = 0
1956 do index1=maxIndex-1,maxIndex+1
1957 index2 = index2 + 1
1958 A(index2,1) = lnp(index1-1)*lnp(index1-1)
1959 A(index2,2) = lnp(index1-1)
1960 A(index2,3) = 1.0d0
1961 B(index2,1) = deriv(index1)
1962 end do
1963
1964 call dgesv(N,nrhs,A,lda,ipiv,B,ldb,info)
1965
1966 if (info==0) then
1967 maxLnP = -0.5*(B(2,1)/B(1,1))
1968 else
1969 maxLnP = lnp(maxIndex-1)
1970 end if
1971
1972 end subroutine get_peak
1973
1974end module ensembleObservations_mod
