1module innovation_mod
2 ! MODULE innovation_mod (prefix='inn' category='1. High-level functionality')
3 !
4 !:Purpose: Several high-level subroutines used to compute the innovations:
5 ! that is, the observation-minus-background values. This includes
6 ! the subroutine that reads in the gridded high-res background state
7 ! from standard files.
8 !
9 use midasMpi_mod
10 use obsSpaceData_mod
11 use columnData_mod
12 use timeCoord_mod
13 use obsTimeInterp_mod
14 use obsOperators_mod
15 use mathPhysConstants_mod
16 use horizontalCoord_mod
17 use varNameList_mod
18 use verticalCoord_mod
19 use gridStateVector_mod
20 use utilities_mod
21 use message_mod
22 use obsFilter_mod
23 use gps_mod
24 use tovsNL_mod
25 use multiIRbgck_mod
26 use obsFiles_mod
27 use randomNumber_mod
28 use obsErrors_mod
29 use bufr_mod
30 use statetocolumn_mod
31 use biascorrectionSat_mod
32 use rmatrix_mod
33 use costFunction_mod
34 use varqc_mod
35 use humidityLimits_mod
36 use interpolation_mod
37 implicit none
38 save
39 private
40
41 ! public procedures
42 public :: inn_setupObs, inn_computeInnovation
43 public :: inn_perturbObs, inn_setupColumnsOnTrlLev, inn_setupColumnsOnAnlIncLev
44 public :: inn_getHcoVcoFromTrlmFile
45
46 character(len=48) :: innovationMode
47
48contains
49
50 !--------------------------------------------------------------------------
51 ! inn_setupObs
52 !--------------------------------------------------------------------------
53 subroutine inn_setupobs(obsSpaceData, hco_anl, obsColumnMode, obsMpiStrategy, &
54 innovationMode_in, obsClean_opt )
55 !
56 !:Purpose: To initialize the observation parameters and constants
57 !
58 implicit none
59
60 ! Arguments:
61 type(struct_obs), intent(out) :: obsSpaceData
62 type(struct_hco), pointer, intent(in) :: hco_anl
63 character(len=*), intent(in) :: obsMpiStrategy
64 character(len=*), intent(in) :: obsColumnMode
65 character(len=*), intent(in) :: innovationMode_in
66 logical, optional, intent(in) :: obsClean_opt
67
68 ! Locals:
69 character(len=20) :: nameDimFile
70 integer :: get_max_rss, ierr, fnom, fclos, unitDimFile, mxstn, mxobs
71 logical :: obsDimFileExists
72
73 write(*,FMT=9000)
749000 FORMAT(/,1x,' INN_SETUPOBS - Initialisation of observations',/,1x,3('- -----------'))
75
76 call utl_tmg_start(10,'--Observations')
77
78 !
79 !- Setup de the mode
80 !
81 innovationMode = innovationMode_in
82
83 !
84 !- Specify the active observation-array columns
85 !
86 call obs_class_initialize(obsColumnMode)
87 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
88 !
89 !- Allocate memory for observation arrays
90 !
91 nameDimFile = './obs/cmadim'
92 inquire( file=trim(nameDimFile), exist=obsDimFileExists )
93 if ( obsDimFileExists ) then
94 unitDimFile = 0
95 ierr = fnom( unitDimFile, trim(nameDimFile), 'FTN+SEQ+R/O', 0 )
96 read(unitDimFile,*) mxstn
97 read(unitDimFile,*) mxobs
98 ierr=fclos(unitDimFile)
99 call obs_initialize(obsSpaceData, numHeader_max_opt=mxstn, numBody_max_opt=mxobs, &
100 mpi_local_opt=obsf_filesSplit())
101 else
102 call obs_initialize(obsSpaceData, mpi_local_opt=obsf_filesSplit())
103 end if
104 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
105
106 !
107 !- Set up the list of elements to be assimilated and flags for rejection
108 !
109 call filt_setup(innovationMode) ! IN
110
111 !
112 !- Initialize TOVS processing
113 !
114 call tvs_setup
115
116 !
117 !- Read the observations from files
118 !
119 call utl_tmg_start(11,'----ReadObsFiles')
120 call obsf_readFiles( obsSpaceData )
121 call utl_tmg_stop(11)
122
123 !
124 !- Initialize GPS processing
125 !
126 if (obs_famExist(obsSpaceData,'RO')) call gps_setupro
127 if (obs_famExist(obsSpaceData,'GP')) call gps_setupgb
128
129 !
130 !- Filter out data from the obs data base
131 !
132 call filt_suprep(obsSpaceData)
133
134 ! Additional filtering for bias correction if requested
135 call bcs_setup()
136 call bcs_filterObs(obsSpaceData)
137
138 if ( present(obsClean_opt) ) then
139 if ( obsClean_opt ) then
140 write(*,*) ''
141 write(*,*) 'inn_setupObs: !!WARNING!! Performing a cleanup of obsSpaceData'
142 write(*,*) ' This may make it impossible to update burp files'
143 write(*,*) ''
144 call obs_clean2(obsSpaceData)
145 end if
146 end if
147
148 !
149 !- Set (OBS_IPC and OBS_IPT) or OBS_IP columns according to the chosen strategy
150 !
151 write(*,*)
152 write(*,*) 'inn_setupObs - Using obsMpiStrategy = ', trim(obsMpiStrategy)
153 call setObsMpiStrategy(obsSpaceData,hco_anl,obsMpiStrategy)
154
155 !
156 !- Check if burp files already split
157 !
158 if ( obsf_filesSplit() ) then
159 ! local observations files, so just do reallocation to reduce memory used
160 call obs_squeeze(obsSpaceData)
161 if ( obs_columnActive_IH(obsSpaceData,OBS_IPC) ) then
162 call obs_MpiRedistribute(obsSpaceData,OBS_IPC)
163 end if
164 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
165 else
166 ! complete set of obs on each MPI process, only keep subset according to OBS_IP
167 call obs_reduceToMpiLocal(obsSpaceData)
168 end if
169
170 !
171 !- Initialization and memory allocation for TOVS processing
172 !
173 if (obs_famExist(obsSpaceData,'TO') .and. (trim(innovationMode) /= 'thinning')) then
174 call tvs_setupAlloc(obsSpaceData)
175 if (trim(innovationMode) == 'bgck' ) call irbg_setup()
176 ! Initialize non diagonal observation error matrices
177 if ( trim(innovationMode) == 'analysis' .or. trim(innovationMode) == 'FSO') call oer_setInterchanCorr()
178 end if
179
180 call utl_tmg_stop(10)
181
182 end subroutine inn_setupobs
183
184 !--------------------------------------------------------------------------
185 ! inn_setupColumnsOnTrlLev
186 !--------------------------------------------------------------------------
187 subroutine inn_setupColumnsOnTrlLev( columnTrlOnTrlLev, obsSpaceData, hco_core, &
188 stateVectorUpdateHighRes, &
189 deallocInterpInfoNL_opt )
190 !
191 !:Purpose: To compute vertical (and potentially slanted) columns of trial data interpolated to obs location
192 !
193 implicit none
194
195 ! Arguments:
196 type(struct_columnData), intent(out) :: columnTrlOnTrlLev
197 type(struct_obs), intent(inout) :: obsSpaceData
198 type(struct_hco), pointer, intent(in) :: hco_core
199 type(struct_gsv), intent(inout) :: stateVectorUpdateHighRes
200 logical , optional, intent(in) :: deallocInterpInfoNL_opt
201
202 ! Locals:
203 type(struct_vco), pointer :: vco_trl => null()
204 integer :: ierr, nulnam, fnom, fclos
205 logical :: deallocInterpInfoNL
206 real(8), pointer :: onecolumn(:)
207 character(len=20) :: timeInterpType_nl ! 'NEAREST' or 'LINEAR'
208 integer :: numObsBatches ! number of batches for calling interp setup
209
210 NAMELIST /NAMINN/timeInterpType_nl, numObsBatches
211
212 write(*,*)
213 write(*,*) 'inn_setupColumnsOnTrlLev: START'
214
215 timeInterpType_nl = 'NEAREST'
216 numObsBatches = 20
217
218 if (utl_isNamelistPresent('naminn','./flnml')) then
219 nulnam = 0
220 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
221 if (ierr /= 0) call utl_abort('inn_setupColumnsOnTrlLev: Error opening file flnml')
222 read(nulnam,nml=naminn,iostat=ierr)
223 if (ierr /= 0) call utl_abort('inn_setupColumnsOnTrlLev: Error reading namelist')
224 if (mmpi_myid == 0) write(*,nml=naminn)
225 ierr = fclos(nulnam)
226 else
227 write(*,*)
228 write(*,*) 'inn_setupColumnsOnTrlLev: Namelist block NAMINN is missing in the namelist.'
229 write(*,*) ' The default values will be taken.'
230 if (mmpi_myid == 0) write(*,nml=naminn)
231 end if
232
233 if ( present(deallocInterpInfoNL_opt) ) then
234 deallocInterpInfoNL = deallocInterpInfoNL_opt
235 else
236 deallocInterpInfoNL = .true.
237 end if
238
239 nullify(vco_trl)
240 vco_trl => gsv_getVco(stateVectorUpdateHighRes)
241
242 call col_setVco(columnTrlOnTrlLev,vco_trl)
243 call col_allocate(columnTrlOnTrlLev,obs_numHeader(obsSpaceData),mpiLocal_opt=.true.)
244
245 ! copy latitude from obsSpaceData
246 if ( obs_numHeader(obsSpaceData) > 0 ) then
247 call obs_extractObsRealHeaderColumn(columnTrlOnTrlLev%lat(:), obsSpaceData, OBS_LAT)
248 end if
249
250 call s2c_nl( stateVectorUpdateHighRes, obsSpaceData, columnTrlOnTrlLev, hco_core, &
251 timeInterpType=timeInterpType_nl, &
252 moveObsAtPole_opt=.true., numObsBatches_opt=numObsBatches, &
253 dealloc_opt=deallocInterpInfoNL )
254
255 if ( col_getNumCol(columnTrlOnTrlLev) > 0 .and. col_varExist(columnTrlOnTrlLev,'Z_T ') ) then
256 write(*,*) 'inn_setupBackgroundColumns, statevector->Column 1:'
257 write(*,*) 'Z_T:'
258 onecolumn => col_getColumn(columnTrlOnTrlLev,1,'Z_T ')
259 write(*,*) onecolumn(:)
260 write(*,*) 'Z_M:'
261 onecolumn => col_getColumn(columnTrlOnTrlLev,1,'Z_M ')
262 write(*,*) onecolumn(:)
263
264 nullify(onecolumn)
265 end if
266 if ( col_getNumCol(columnTrlOnTrlLev) > 0 .and. col_varExist(columnTrlOnTrlLev,'P_T ') ) then
267 write(*,*) 'inn_setupBackgroundColumns, statevector->Column 1:'
268 write(*,*) 'P_T:'
269 onecolumn => col_getColumn(columnTrlOnTrlLev,1,'P_T ')
270 write(*,*) onecolumn(:)
271 write(*,*) 'P_M:'
272 onecolumn => col_getColumn(columnTrlOnTrlLev,1,'P_M ')
273 write(*,*) onecolumn(:)
274
275 nullify(onecolumn)
276 end if
277
278 write(*,*) 'inn_setupColumnsOnTrlLev: END'
279
280 end subroutine inn_setupColumnsOnTrlLev
281
282 !--------------------------------------------------------------------------
283 ! inn_setupColumnsOnAnlIncLev
284 !--------------------------------------------------------------------------
285 subroutine inn_setupColumnsOnAnlIncLev(columnTrlOnTrlLev,columnTrlOnAnlIncLev)
286 !
287 !:Purpose: To create trial data columns on analysis increment levels
288 !
289 implicit none
290
291 ! Arguments:
292 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev
293 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
294
295 ! Locals:
296 integer :: jvar, jlev, columnIndex
297 real(8), pointer :: columnTrlOnAnlIncLev_ptr(:), columnTrlOnTrlLev_ptr(:)
298
299 if (col_getNumCol(columnTrlOnAnlIncLev) == 0) return
300
301 call msg('inn_setupColumnsOnAnlIncLev','START',verb_opt=2)
302
303 !
304 !- Data copying from columnh to columnTrlOnAnlIncLev
305 !
306
307 ! copy latitude
308 columnTrlOnAnlIncLev%lat(:) = columnTrlOnTrlLev%lat(:)
309
310 ! copy 2D surface variables
311 do jvar = 1, vnl_numvarmax2D
312 if ( .not. col_varExist(columnTrlOnAnlIncLev,vnl_varNameList2D(jvar)) ) cycle
313 do columnIndex = 1, col_getNumCol(columnTrlOnAnlIncLev)
314 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev , columnIndex, vnl_varNameList2D(jvar))
315 columnTrlOnTrlLev_ptr => col_getColumn(columnTrlOnTrlLev, columnIndex, vnl_varNameList2D(jvar))
316 columnTrlOnAnlIncLev_ptr(:) = columnTrlOnTrlLev_ptr(:)
317 end do
318 end do
319
320 !
321 !- Vertical interpolation of 3D variables from trials levels to analysis increment levels
322 !
323 do jvar = 1, vnl_numvarmax3D
324
325 if ( .not. col_varExist(columnTrlOnAnlIncLev,vnl_varNameList3D(jvar)) ) cycle
326
327 call int_vInterp_col( columnTrlOnTrlLev, columnTrlOnAnlIncLev, &
328 vnl_varNameList3D(jvar), useColumnPressure_opt=.false.)
329
330 if ( vnl_varNameList3D(jvar) == 'HU ') then
331 ! Imposing a minimum value for HU
332 do columnIndex = 1, col_getNumCol(columnTrlOnAnlIncLev)
333 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,columnIndex,'HU')
334 do jlev = 1, col_getNumLev(columnTrlOnAnlIncLev,'TH')
335 columnTrlOnAnlIncLev_ptr(jlev) = max(columnTrlOnAnlIncLev_ptr(jlev),col_rhumin)
336 end do
337 end do
338
339 ! Imposing minimum/maximum value for ALL cloud variables at all levels
340 else if (vnl_isCloudVar(vnl_varNameList3D(jvar))) then
341 do columnIndex = 1, col_getNumCol(columnTrlOnAnlIncLev)
342 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,columnIndex,vnl_varNameList3D(jvar))
343 do jlev = 1, col_getNumLev(columnTrlOnAnlIncLev,'TH')
344 columnTrlOnAnlIncLev_ptr(jlev) = max(columnTrlOnAnlIncLev_ptr(jlev), &
345 qlim_getMinValueCloud(vnl_varNameList3D(jvar)))
346 columnTrlOnAnlIncLev_ptr(jlev) = min(columnTrlOnAnlIncLev_ptr(jlev), &
347 qlim_getMaxValueCloud(vnl_varNameList3D(jvar)))
348 end do
349 end do
350
351 else if (trim(vnl_varKindFromVarname(vnl_varNameList3D(jvar))) == 'CH') then
352 ! Imposing boundary values for CH kind variables. This is to prevent
353 ! undesired values usually from vertical extrapolation.
354 do columnIndex = 1, col_getNumCol(columnTrlOnAnlIncLev)
355 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,columnIndex,trim(vnl_varNameList3D(jvar)))
356 if ( col_minValVarKindCH(vnl_varListIndex(vnl_varNameList3D(jvar))) > 1.01*MPC_missingValue_R8 ) then
357 do jlev=1,col_getNumLev(columnTrlOnAnlIncLev,'TH')
358 columnTrlOnAnlIncLev_ptr(jlev) = max(columnTrlOnAnlIncLev_ptr(jlev),col_minValVarKindCH(vnl_varListIndex(vnl_varNameList3D(jvar))))
359 end do
360 end if
361 end do
362 end if
363
364 end do
365
366 ! Print pressure on thermo levels for the first column
367 if (col_varExist(columnTrlOnAnlIncLev,'P_T')) then
368 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,1,'P_T')
369 call msg('inn_setupColumnsOnAnlIncLev', ' after vintprof:'&
370 //new_line('')//'P_T (columnTrlOnAnlIncLev(1)) = '//str(columnTrlOnAnlIncLev_ptr(:)), mpiAll_opt=.false.)
371 end if
372
373 !
374 !- Height adjustments
375 !
376
377 ! Set surface height
378 do columnIndex = 1, col_getNumCol(columnTrlOnAnlIncLev)
379 columnTrlOnAnlIncLev%heightSfc(columnIndex) = columnTrlOnTrlLev%heightSfc(columnIndex)
380 end do
381
382 ! Remove the height offset for the diagnostic levels for backward compatibility only
383 if ( col_varExist(columnTrlOnAnlIncLev,'Z_T') .and. .not.columnTrlOnAnlIncLev%addHeightSfcOffset ) then
384 do columnIndex = 1, col_getNumCol(columnTrlOnAnlIncLev)
385 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,columnIndex,'Z_T')
386 columnTrlOnAnlIncLev_ptr(col_getNumLev(columnTrlOnAnlIncLev,'TH')) = columnTrlOnAnlIncLev%heightSfc(columnIndex)
387 end do
388 end if
389 if ( col_varExist(columnTrlOnAnlIncLev,'Z_M') .and. .not.columnTrlOnAnlIncLev%addHeightSfcOffset ) then
390 do columnIndex = 1, col_getNumCol(columnTrlOnAnlIncLev)
391 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,columnIndex,'Z_M')
392 columnTrlOnAnlIncLev_ptr(col_getNumLev(columnTrlOnAnlIncLev,'MM')) = columnTrlOnAnlIncLev%heightSfc(columnIndex)
393 end do
394 end if
395
396 ! Print height info of the first original and interpolated columns
397 if ( col_getNumLev(columnTrlOnAnlIncLev,'TH') > 0 .and. col_varExist(columnTrlOnAnlIncLev,'Z_T') ) then
398 columnTrlOnTrlLev_ptr => col_getColumn(columnTrlOnTrlLev,1,'Z_T')
399 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,1,'Z_T')
400 call msg('inn_setupColumnsOnAnlIncLev','vIntProf output:' &
401 //new_line('')//'Z_T (columnTrlOnTrlLev) = '//str(columnTrlOnTrlLev_ptr(:)) &
402 //new_line('')//'Z_T (columnTrlOnAnlIncLev) = '//str(columnTrlOnAnlIncLev_ptr(:)))
403 end if
404
405 if ( col_getNumLev(columnTrlOnAnlIncLev,'MM') > 0 .and. col_varExist(columnTrlOnAnlIncLev,'Z_M') ) then
406 columnTrlOnTrlLev_ptr => col_getColumn(columnTrlOnTrlLev,1,'Z_M')
407 columnTrlOnAnlIncLev_ptr => col_getColumn(columnTrlOnAnlIncLev,1,'Z_M')
408 call msg('inn_setupColumnsOnAnlIncLev','vIntProf output:'&
409 //new_line('')//'Z_M (columnTrlOnTrlLev) = '//str(columnTrlOnTrlLev_ptr(:)) &
410 //new_line('')//'Z_M (columnTrlOnAnlIncLev) = '//str(columnTrlOnAnlIncLev_ptr(:)))
411 end if
412
413 call msg('inn_setupColumnsOnAnlIncLev','HeightSfc:'//str(columnTrlOnAnlIncLev%heightSfc(1)))
414
415 call msg('inn_setupColumnsOnAnlIncLev','END',verb_opt=2)
416
417 end subroutine inn_setupColumnsOnAnlIncLev
418
419 !--------------------------------------------------------------------------
420 ! inn_computeInnovation
421 !--------------------------------------------------------------------------
422 subroutine inn_computeInnovation( columnTrlOnTrlLev, obsSpaceData, filterObsAndInitOer_opt, &
423 applyVarqcOnNlJo_opt, destObsColumn_opt, &
424 beSilent_opt, callFiltTopo_opt, callSetErrGpsgb_opt, &
425 analysisMode_opt )
426 !
427 !:Purpose: To initialize observation innovations using the nonlinear H
428 !
429 implicit none
430
431 ! Arguments:
432 type(struct_columnData), intent(in) :: columnTrlOnTrlLev
433 type(struct_obs) , intent(inout) :: obsSpaceData
434 logical, optional , intent(in) :: filterObsAndInitOer_opt
435 logical, optional , intent(in) :: applyVarqcOnNlJo_opt
436 integer, optional , intent(in) :: destObsColumn_opt ! column where result stored, default is OBS_OMP
437 logical, optional , intent(in) :: beSilent_opt
438 logical, optional , intent(in) :: callFiltTopo_opt ! whether to make call to FiltTopo
439 logical, optional , intent(in) :: callSetErrGpsgb_opt ! whether to make call to oer_SETERRGPSGB
440 logical, optional , intent(in) :: analysisMode_opt ! analysisMode argument for oer_SETERRGPSGB and oop_gpsgb_nl
441
442 ! Locals:
443 real(8) :: Jo
444 integer :: destObsColumn, get_max_rss
445 logical :: applyVarqcOnNlJo, filterObsAndInitOer, beSilent, callFiltTopo, callSetErrGpsgb, analysisMode
446 logical, save :: lgpdata = .false.
447
448 call utl_tmg_start(10,'--Observations')
449
450 if ( present(beSilent_opt) ) then
451 beSilent = beSilent_opt
452 else
453 beSilent = .false.
454 end if
455
456 if ( .not. beSilent ) then
457 write(*,*)
458 write(*,*) '--Starting subroutine inn_computeInnovation--'
459 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
460 end if
461
462 if ( present(filterObsAndInitOer_opt) ) then
463 filterObsAndInitOer = filterObsAndInitOer_opt
464 else
465 filterObsAndInitOer = .true.
466 end if
467
468 if ( present(applyVarqcOnNlJo_opt) ) then
469 applyVarqcOnNlJo = applyVarqcOnNlJo_opt
470 else
471 applyVarqcOnNlJo = .false.
472 end if
473
474 if ( present(destObsColumn_opt) ) then
475 destObsColumn = destObsColumn_opt
476 else
477 destObsColumn = obs_omp
478 end if
479
480 if ( present(callFiltTopo_opt) ) then
481 callFiltTopo = callFiltTopo_opt
482 else
483 callFiltTopo = .true.
484 end if
485
486 if ( present(callSetErrGpsgb_opt ) ) then
487 callSetErrGpsgb = callSetErrGpsgb_opt
488 else
489 callSetErrGpsgb = .true.
490 end if
491
492 if ( present(analysisMode_opt ) ) then
493 analysisMode = analysisMode_opt
494 else
495 analysisMode = .true.
496 end if
497
498 if ( .not.beSilent ) write(*,*) 'oti_timeBinning: Before filtering done in inn_computeInnovation'
499 if ( .not.beSilent ) call oti_timeBinning(obsSpaceData,tim_nstepobs)
500
501 ! Reject observed elements too far below the surface. Pressure values
502 ! for elements slightly below the surface are replaced by the surface
503 ! pressure values of the trial field.
504 !
505 ! GB-GPS (met and ZTD) observations are processed in s/r filt_topoSFC (in obsFilter_mod.ftn90)
506 !
507 if ( filterObsAndInitOer .and. callFiltTopo ) then
508 call filt_topo(columnTrlOnTrlLev,obsSpaceData,beSilent)
509 else
510 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,*) 'inn_computeInnovation: skip filt_topo'
511 end if
512
513 ! Remove surface station wind observations
514 if ( trim(innovationMode) == 'analysis' .or. trim(innovationMode) == 'FSO' ) then
515 if ( filterObsAndInitOer ) then
516 call filt_surfaceWind(obsSpaceData, beSilent)
517 else
518 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,*) 'inn_computeInnovation: skip filt_surfaceWind'
519 end if
520 end if
521
522 ! Find interpolation layer in model profiles
523 if ( col_getNumLev(columnTrlOnTrlLev,'MM') > 1 ) call oop_vobslyrs(columnTrlOnTrlLev, obsSpaceData, beSilent)
524
525 !
526 !- Calculate the innovations [Y - H(Xb)] and place the result in obsSpaceData in destObsColumn column
527 !
528 call utl_tmg_start(17,'----ObsOper_NL')
529
530 ! Radiosondes
531 call oop_ppp_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'UA', destObsColumn)
532
533 ! Aircrafts
534 call oop_ppp_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'AI', destObsColumn)
535
536 ! SatWinds
537 if ( filterObsAndInitOer ) then
538 call oer_sw(columnTrlOnTrlLev,obsSpaceData)
539 else
540 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,*) 'inn_computeInnovation: skip oer_sw'
541 end if
542
543 call oop_ppp_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'SW', destObsColumn)
544
545 ! Surface (SF, UA, SC, GP and RA families)
546 call oop_sfc_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'SF', destObsColumn)
547 call oop_sfc_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'UA', destObsColumn)
548 call oop_sfc_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'SC', destObsColumn)
549 call oop_sfc_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'GP', destObsColumn)
550 call oop_sfc_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'RA', destObsColumn)
551
552 ! RADAR Doppler velocity
553 if ( filterObsAndInitOer ) then
554 ! Filter Radar for Doppler velocity
555 call filt_radvel(columnTrlOnTrlLev, obsSpaceData, beSilent)
556 else
557 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,*) 'inn_computeInnovation: skip filt_radvel'
558 end if
559
560 call oop_raDvel_nl(columnTrlOnTrlLev,obsSpaceData, beSilent, 'RA', destObsColumn)
561
562 ! Sea surface temperature
563 call oop_sst_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'TM', destObsColumn)
564
565 ! Sea ice concentration
566 if ( filterObsAndInitOer ) then
567 call filt_iceConcentration(obsSpaceData, beSilent)
568 call filt_backScatAnisIce(obsSpaceData, beSilent)
569 call oer_setErrBackScatAnisIce(obsSpaceData, beSilent, columnTrlOnTrlLev_opt=columnTrlOnTrlLev)
570 else
571 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,*) 'inn_computeInnovation: skip filt_iceConcentration, filt_backScatAnisIce, and oer_setErrBackScatAnisIce'
572 end if
573
574 call oop_ice_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'GL', destObsColumn)
575
576 ! Hydrology
577 call oop_hydro_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'HY', destObsColumn)
578
579 ! TOVS / Radiances
580 if (trim(innovationMode) == 'bgck' ) then
581 call oop_tovs_nl(columnTrlOnTrlLev, obsSpaceData, tim_getDatestamp(), &
582 beSilent, bgckMode_opt=.true., destObs_opt=destObsColumn)
583 else
584 call oop_tovs_nl(columnTrlOnTrlLev, obsSpaceData, tim_getDatestamp(), &
585 beSilent, bgckMode_opt=.false., destObs_opt=destObsColumn)
586 end if
587
588 ! Profilers
589 call oop_zzz_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'PR', destObsColumn)
590
591 ! Aladin winds
592 call oop_zzz_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, 'AL', destObsColumn)
593
594 ! GPS radio occultation
595 if (obs_famExist(obsSpaceData,'RO', localMPI_opt = .true. )) then
596 if ( filterObsAndInitOer ) then
597 call filt_gpsro(columnTrlOnTrlLev, obsSpaceData, beSilent)
598 call oer_SETERRGPSRO(columnTrlOnTrlLev, obsSpaceData, beSilent)
599 else
600 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,*) 'inn_computeInnovation: skip filt_gpsro, and oer_SETERRGPSRO'
601 end if
602 call oop_gpsro_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, destObsColumn)
603 end if
604
605 ! Chemical constituents
606 call oop_chm_nl(columnTrlOnTrlLev, obsSpaceData, destObsColumn)
607
608 ! GPS ground-based zenith delay
609 if (obs_famExist(obsSpaceData,'GP', localMPI_opt = .true. )) then
610 if ( CallSetErrGpsgb ) then
611 call oer_SETERRGPSGB(columnTrlOnTrlLev, obsSpaceData, beSilent, lgpdata, analysisMode)
612 else
613 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,*) 'inn_computeInnovation: skip oer_SETERRGPSGB'
614 end if
615 if (lgpdata) call oop_gpsgb_nl(columnTrlOnTrlLev, obsSpaceData, beSilent, &
616 destObsColumn, analysisMode_opt=analysisMode)
617 end if
618
619 call utl_tmg_stop(17)
620
621 ! Save as OBS_WORK : R**-1/2 (d)
622 call rmat_RsqrtInverseAllObs(obsSpaceData,OBS_WORK,destObsColumn)
623
624 ! Store J-obs in OBS_JOBS : 1/2 * R**-1 (d)**2
625 call cfn_calcJo(obsSpaceData)
626
627 ! applying varqc, if asked for
628 if ( applyVarqcOnNlJo ) call vqc_NlTl(obsSpaceData)
629
630 ! Compute Jo components and print
631 call cfn_sumJo(obsSpaceData, Jo, beSilent_opt=beSilent)
632 if ( mmpi_myid == 0 .and. .not. beSilent ) write(*,'(a15,f25.17)') 'Total Jo = ',Jo
633
634 if ( .not.beSilent ) write(*,*) 'oti_timeBinning: After filtering done in inn_computeInnovation'
635 if ( .not.beSilent ) call oti_timeBinning(obsSpaceData,tim_nstepobs)
636
637 if ( .not. beSilent ) then
638 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
639 write(*,*) '--Done subroutine inn_computeInnovation--'
640 end if
641
642 call utl_tmg_stop(10)
643
644 end subroutine inn_computeInnovation
645
646 !--------------------------------------------------------------------------
647 ! setObsMpiStrategy
648 !--------------------------------------------------------------------------
649 subroutine setObsMpiStrategy(obsSpaceData, hco_anl, mpiStrategy)
650 !
651 !:Purpose: To distribute header indices following the chosen strategy,
652 ! current options: "LIKESPLITFILES", "ROUNDROBIN", "LATLONTILES".
653 !
654 implicit none
655
656 ! Arguments:
657 type(struct_obs), intent(inout) :: obsSpaceData
658 type(struct_hco), pointer, intent(in) :: hco_anl
659 character(len=*), intent(in) :: mpiStrategy
660
661 ! Locals:
662 real(8) :: lat_r8, lon_r8
663 real :: lat_r4, lon_r4
664 real :: xpos_r4, ypos_r4
665 integer :: numHeaderFile, headerIndex, latIndex, lonIndex, ierr
666 integer :: IP, IP_x, IP_y
667 integer :: gdxyfll
668
669 !
670 !- Determine obs_ipc (column) and obs_ipt (tile) according to distribution strategy
671 !
672 write(*,*)
673 numHeaderFile = obs_numheader(obsSpaceData)
674 write(*,*) 'setObsMpiStrategy: numHeader after reading files = ',numHeaderFile
675 write(*,*) 'setObsMpiStrategy: strategy = ',trim(mpiStrategy)
676
677 ! make sure old PE index is not used (set to -1)
678 do headerIndex = 1, numHeaderFile
679 call obs_headSet_i(obsSpaceData,OBS_IP,headerIndex,-1)
680 end do
681
682 ! set PE index for the file (where reading and updating was/will be done)
683 if ( obs_columnActive_IH(obsSpaceData,OBS_IPF) ) then
684 do headerIndex = 1, numHeaderFile
685 call obs_headSet_i(obsSpaceData,OBS_IPF,headerIndex,mmpi_myid)
686 end do
687 end if
688
689 select case (trim(mpiStrategy))
690 case ('LIKESPLITFILES')
691 !- Keep distribution exactly as it is in the split files:
692 do headerIndex = 1, numHeaderFile
693 if ( obs_columnActive_IH(obsSpaceData,OBS_IPC) ) &
694 call obs_headSet_i(obsSpaceData,OBS_IPC,headerIndex, mmpi_myid)
695 if ( obs_columnActive_IH(obsSpaceData,OBS_IPT) ) &
696 call obs_headSet_i(obsSpaceData,OBS_IPT,headerIndex, mmpi_myid)
697 if ( obs_columnActive_IH(obsSpaceData,OBS_IP) ) &
698 call obs_headSet_i(obsSpaceData,OBS_IP,headerIndex, mmpi_myid)
699 end do
700 case ('ROUNDROBIN')
701 !- Distribute by a round-robin strategy for both obs_ipc and obs_ipt:
702 ! (Only use if files already not split by round robin)
703 do headerIndex = 1, numHeaderFile
704 IP = mod((headerIndex-1),mmpi_nprocs)
705 if ( obs_columnActive_IH(obsSpaceData,OBS_IPC) ) &
706 call obs_headSet_i(obsSpaceData,OBS_IPC,headerIndex, IP)
707 if ( obs_columnActive_IH(obsSpaceData,OBS_IPT) ) &
708 call obs_headSet_i(obsSpaceData,OBS_IPT,headerIndex, IP)
709 if ( obs_columnActive_IH(obsSpaceData,OBS_IP) ) &
710 call obs_headSet_i(obsSpaceData,OBS_IP ,headerIndex, IP)
711 end do
712 case ('LATLONTILES')
713 !- Distribute by latitude/longitude tiles for both obs_ipc and obs_ipt:
714 do headerIndex = 1, numHeaderFile
715 ! compute grid index for each observation header
716 lat_r8 = obs_headElem_r(obsSpaceData,OBS_LAT,headerIndex)
717 lon_r8 = obs_headElem_r(obsSpaceData,OBS_LON,headerIndex)
718 lat_r4 = real(lat_r8) * MPC_DEGREES_PER_RADIAN_R4
719 lon_r4 = real(lon_r8) * MPC_DEGREES_PER_RADIAN_R4
720 ierr = gdxyfll( hco_anl%EZscintID, & ! IN
721 xpos_r4, ypos_r4, & ! OUT
722 lat_r4, lon_r4, 1 ) ! IN
723
724 ! compute correponding mpi task id for each observation
725 latIndex = floor(ypos_r4)
726 lonIndex = floor(xpos_r4)
727 IP_y = mmpi_myidYfromLat(latIndex, hco_anl%nj)
728 IP_x = mmpi_myidXfromLon(lonIndex, hco_anl%ni)
729 IP = IP_x + IP_y*mmpi_npex
730
731 call obs_headSet_i(obsSpaceData,OBS_IPC,headerIndex, IP)
732 call obs_headSet_i(obsSpaceData,OBS_IPT,headerIndex, IP)
733 end do
734 case ('LATLONTILESBALANCED')
735 !- Distribute by latitude/longitude tiles, but with simple & cheap balancing for obs_ipc (1 send or recv):
736 call utl_abort('setObsMpiStrategy: Sorry, LATLONTILESBALANCED no longer available')
737 case default
738 write(*,*)
739 write(*,*) 'ERROR unknown mpiStrategy: ', trim(mpiStrategy)
740 call utl_abort('setObsMpiStrategy')
741 end select
742
743 end subroutine setObsMpiStrategy
744
745 !--------------------------------------------------------------------------
746 ! inn_perturbObs
747 !--------------------------------------------------------------------------
748 subroutine inn_perturbObs(obsSpaceData,numAnalyses,indexAnalysis, &
749 indexBatch,obs_column_index_src,obs_column_index_dest)
750 !
751 !:Purpose: To perturb the innovation vector to simulate effect of
752 ! observation uncertainty
753 !
754 !.. WARNING:: perturbations are not the same when MPI topology changes!!!
755 !
756 implicit none
757
758 ! Arguments:
759 type(struct_obs), intent(inout) :: obsSpaceData
760 integer, intent(in) :: numAnalyses
761 integer, intent(in) :: indexAnalysis
762 integer, intent(in) :: indexBatch
763 integer, intent(in) :: obs_column_index_src
764 integer, intent(in) :: obs_column_index_dest
765
766 ! Locals:
767 integer :: numPerturbations
768 integer :: nrandseed,iseed,indexAnalysis2,indexBody,indexFamily
769 integer, parameter :: numFamily=10
770 real*8 :: zmean,originalOmp
771 real*8 :: scaleFactor(numFamily)
772 character(len=2) :: familyList(numFamily)
773 real*8, save, pointer :: obsPerturbations(:,:) => NULL()
774 logical, save :: firstTime = .true.
775
776 familyList(1)='UA' ; scaleFactor(1)=1.00d0
777 familyList(2)='AI' ; scaleFactor(2)=1.00d0
778 familyList(3)='SF' ; scaleFactor(3)=1.00d0
779 familyList(4)='TO' ; scaleFactor(4)=1.00d0
780 familyList(5)='SW' ; scaleFactor(5)=1.00d0
781 familyList(6)='SC' ; scaleFactor(6)=1.00d0
782 familyList(7)='PR' ; scaleFactor(7)=1.00d0
783 familyList(8)='RO' ; scaleFactor(8)=1.00d0
784 familyList(9)='GP' ; scaleFactor(9)=1.00d0
785 familyList(10)='RA'; scaleFactor(10)=1.00d0
786
787 numPerturbations = numAnalyses
788
789 if (firstTime) then
790
791 if (.not.associated(obsPerturbations)) then
792 write(*,*) 'perturbObs: allocating space for all perturbations'
793 allocate(obsPerturbations(obs_numBody(obsSpaceData),numPerturbations))
794 end if
795
796 write(*,*) 'perturbObs: computing random numbers'
797
798 ! compute random perturbations
799 do indexAnalysis2 = 1,numPerturbations
800 nrandseed = indexAnalysis2 + (indexBatch-1)*numAnalyses
801 iseed = ABS(nrandseed)
802 write(*,*) 'perturbobs: indexAnalysis, iseed=',indexAnalysis2, iseed
803 call rng_setup(iseed) ! JFC: should be called only once, no???
804 do indexBody = 1,obs_numBody(obsSpaceData)
805 obsPerturbations(indexBody,indexAnalysis2)=rng_gaussian()*obs_bodyElem_r(obsSpaceData,OBS_OER,indexBody)
806 end do
807 end do
808
809 ! apply scale factor
810 do indexFamily = 1,numFamily
811 do indexBody = 1,obs_numBody(obsSpaceData)
812 if (obs_getFamily(obsSpaceData,bodyIndex_opt=indexBody).eq.familyList(indexFamily)) then
813 do indexAnalysis2 = 1,numPerturbations
814 obsPerturbations(indexBody,indexAnalysis2)=obsPerturbations(indexBody,indexAnalysis2)*scaleFactor(indexFamily)
815 end do
816 end if
817 end do
818 end do
819
820 ! remove ensemble mean
821 do indexBody = 1,obs_numBody(obsSpaceData)
822 zmean=0.0d0
823 do indexAnalysis2 = 1,numPerturbations
824 zmean=zmean+obsPerturbations(indexBody,indexAnalysis2)
825 end do
826 zmean=zmean/numPerturbations
827 do indexAnalysis2 = 1,numPerturbations
828 obsPerturbations(indexBody,indexAnalysis2)=obsPerturbations(indexBody,indexAnalysis2)-zmean
829 end do
830 end do
831
832 firstTime=.false.
833
834 end if
835
836 ! apply perturbation for current analysis
837 write(*,*) 'perturbObs: applying perturbation for analysis number: ',indexAnalysis
838 do indexBody = 1,obs_numBody(obsSpaceData)
839 if (obs_bodyElem_i(obsSpaceData,OBS_ASS,indexBody) == obs_assimilated) then
840 originalOmp = obs_bodyElem_r(obsSpaceData,obs_column_index_src,indexBody)
841 call obs_bodySet_r(obsSpaceData,obs_column_index_dest,indexBody,originalOmp+obsPerturbations(indexBody,indexAnalysis))
842 end if
843 end do
844
845 end subroutine inn_perturbObs
846
847 !--------------------------------------------------------------------------
848 ! inn_getHcoVcoFromTrlmFile
849 !--------------------------------------------------------------------------
850 subroutine inn_getHcoVcoFromTrlmFile( hco_trl, vco_trl )
851 !
852 !:Purpose: Get hco/vco of the trials
853 !
854 implicit none
855
856 ! Arguments:
857 type(struct_hco), pointer, intent(inout) :: hco_trl
858 type(struct_vco), pointer, intent(inout) :: vco_trl
859
860 ! Locals:
861 character(len=4), pointer :: anlVar(:)
862
863 write(*,*) 'inn_getHcoVcoFromTrlmFile: START'
864 nullify(hco_trl,vco_trl)
865
866 ! check if gsv is initialized.
867 if ( .not. gsv_isInitialized() ) then
868 write(*,*)
869 write(*,*) 'inn_getHcoVcoFromTrlmFile: gsv_setup must be called first. Call it now'
870 call gsv_setup
871 end if
872
873 nullify(anlVar)
874 call gsv_varNamesList(anlVar)
875 call hco_SetupFromFile(hco_trl, './trlm_01', ' ', 'Trial', varName_opt=anlVar(1))
876
877 call vco_SetupFromFile(vco_trl, './trlm_01')
878
879 write(*,*) 'inn_getHcoVcoFromTrlmFile: END'
880
881 end subroutine inn_getHcoVcoFromTrlmFile
882
883end module innovation_mod
