1program midas_diagHBHt
2 !
3 !:Purpose: Main program for computing a single random realization of
4 ! background error in observation space, stored in ``obs_hbht``. This
5 ! can then be used by python or other scripts to compute the
6 ! background error variance (consistent with the specified B matrix)
7 ! in observation space for comparison with the innovation variance
8 ! and observation error variance.
9 !
10 ! --
11 !
12 !:Algorithm: The random realization of background error in observation space
13 ! is computed following these steps:
14 !
15 ! 1. Compute random values for the control vector with each element
16 ! drawn from independent Gaussian distribution with variance of one
17 ! and bias of zero.
18 !
19 ! 2. Multiply random vector by sqrt of B matrix.
20 !
21 ! 3. Apply observation operator to obtain random perturbation in
22 ! observation space.
23 !
24 ! --
25 !
26 !:File I/O: The required input files and produced output files can vary
27 ! according to the application. Below are tables of files for
28 ! typical NWP 4D-EnVar (e.g. GDPS) and sea ice or SST 3D-Var
29 ! applications.
30 !
31 ! --
32 !
33 !============================================== ==============================================================
34 ! Input and Output Files (for NWP application) Description of file
35 !============================================== ==============================================================
36 ! ``flnml`` In - Main namelist file with parameters user may modify
37 ! ``flnml_static`` In - The "static" namelist that should not be modified
38 ! ``trlm_$NN`` (e.g. ``trlm_01``) In - Background state (a.k.a. trial) files for each timestep
39 ! ``analysisgrid`` In - File defining grid for computing the random gridded perturbation
40 ! ``bgcov`` In - Static (i.e. NMC) B matrix file for NWP fields
41 ! ``bgchemcov`` In - Static B matrix file for chemistry fields
42 ! ``ensemble/$YYYYMMDDHH_006_$NNNN`` In - Ensemble member files defining ensemble B matrix
43 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
44 ! ``obserr`` In - Observation error statistics
45 ! ``obsinfo_chm`` In - Something needed for chemistry assimilation?
46 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - Updated obs file for each "family" and MPI task
47 ! Remainder are files related to radiance obs:
48 ! ``stats_$SENSOR_assim`` In - Satellite radiance observation errors of different sensors
49 ! ``stats_tovs`` In - Satellite radiance observation errors
50 ! ``stats_tovs_symmetricObsErr`` In - User-defined symmetric TOVS errors for all sky
51 ! ``ceres_global.std`` In - High-res surface type and water fraction for radiance obs
52 ! ``rtcoef_$PLATFORM_$SENSOR.dat`` In - RTTOV coefficient files
53 ! ``ozoneclim98`` In - Ozone climatology
54 !============================================== ==============================================================
55 !
56 ! --
57 !
58 !:Synopsis: Below is a summary of the ``diagHBHt`` program calling sequence:
59 !
60 ! - **Initial setups:**
61 !
62 ! - Setup horizontal and vertical grid objects for "analysis
63 ! grid" from ``analysisgrid`` file and for "trial grid" from
64 ! first trial file: ``trlm_01``.
65 !
66 ! - Setup ``obsSpaceData`` object and read observations from
67 ! files: ``inn_setupObs``.
68 !
69 ! - Setup ``columnData`` and ``gridStateVector`` modules (read
70 ! list of analysis variables from namelist) and allocate column
71 ! object for storing trial on analysis levels.
72 !
73 ! - Setup the observation error statistics in ``obsSpaceData``
74 ! object: ``oer_setObsErrors``.
75 !
76 ! - Allocate a stateVector object on the trial grid and then
77 ! read the trials: ``gio_readTrials``.
78 !
79 ! - Setup the B matrices: ``bmat_setup``.
80 !
81 ! - Setup the ``gridVariableTransforms``.
82 !
83 ! - **Calculation**
84 !
85 ! - Compute ``columnTrlOnTrlLev`` and ``columnTrlOnAnlIncLev``
86 ! from trials: ``inn_setupColumnsOnTrlLev``,
87 ! ``inn_setupColumnsOnAnlIncLev``
88 !
89 ! - Compute innovation from updated state:
90 ! ``inn_computeInnovation``.
91 !
92 ! - Compute an MPI global random vector, then extract only
93 ! portion needed for this MPI task (to reduce sensitivity of
94 ! results to MPI topology).
95 !
96 ! - Multiply random vector by sqrt of B matrix with resulting
97 ! gridded state random perturbation in ``statevector``.
98 !
99 ! - Apply linearized observation operators to the random gridded
100 ! state: ``s2c_tl`` and ``oop_Htl`` with final result in
101 ! observation space: ``obs_work`` column of ``obsSpaceData``.
102 !
103 ! - Copy result from ``obs_work`` to ``obs_hbht`` column.
104 !
105 ! - **Final steps**, after the outer loop:
106 !
107 ! - Various final steps, including: update the observation files
108 ! (``obsf_writeFiles``).
109 !
110 ! --
111 !
112 !:Options: `List of namelist blocks <../namelists_in_each_program.html#diaghbht>`_
113 ! that can affect the ``diagHBHt`` program.
114 !
115 ! * The choice of what B matrix is used for the calculation is
116 ! controlled for each individual B matrix component through it's
117 ! own namelist block. The weights for all B matrix components are
118 ! zero be default and can be set to a nonzero value through the
119 ! namelist variable ``SCALEFACTOR`` in the namelist block for each
120 ! corresponding fortran module.
121 !
122 ! * All other namelist blocks related to observations are relevant
123 ! for the diagHBHt calculation, including ``NAMFILT`` and
124 ! ``NAMTOV``.
125 !
126 ! * Some of the other relevant namelist blocks used to configure the
127 ! diagHBHt calculation are listed in the following table:
128 !
129 !======================== ============ ==============================================================
130 ! Module Namelist Description of what is controlled
131 !======================== ============ ==============================================================
132 ! ``timeCoord_mod`` ``NAMTIME`` assimilation time window length, temporal resolution of
133 ! the background state and increment (i.e. perturbation)
134 ! ``bMatrixEnsemble_mod`` ``NAMBEN`` weight and other parameters for ensemble-based B matrix
135 ! component
136 ! ``bMatrixHI_mod`` ``NAMBHI`` weight and other parameters for the climatological B matrix
137 ! component based on homogeneous-isotropic covariances
138 ! represented in spectral space
139 ! Other B matrix modules various weight and other parameters for each type of B matrix
140 !======================== ============ ==============================================================
141 !
142 use version_mod
143 use codePrecision_mod
144 use ramDisk_mod
145 use utilities_mod
146 use midasMpi_mod
147 use MathPhysConstants_mod
148 use horizontalCoord_mod
149 use verticalCoord_mod
150 use timeCoord_mod
151 use obsSpaceData_mod
152 use columnData_mod
153 use gridStateVector_mod
154 use gridStateVectorFileIO_mod
155 use controlVector_mod
156 use obsFiles_mod
157 use randomNumber_mod
158 use obsTimeInterp_mod
159 use stateToColumn_mod
160 use innovation_mod
161 use bMatrix_mod
162 use obsErrors_mod
163 use gridVariableTransforms_mod
164 use obsOperators_mod
165 implicit none
166
167 integer :: istamp,exdb,exfin
168 integer :: ierr
169
170 type(struct_obs), target :: obsSpaceData
171 type(struct_columnData), target :: columnTrlOnAnlIncLev
172 type(struct_columnData), target :: columnTrlOnTrlLev
173 type(struct_gsv) :: stateVectorTrialHighRes
174 type(struct_hco), pointer :: hco_trl => null()
175 type(struct_vco), pointer :: vco_trl => null()
176
177 logical :: allocHeightSfc
178
179 character(len=48) :: obsMpiStrategy, varMode
180
181 type(struct_hco), pointer :: hco_anl => null()
182 type(struct_hco), pointer :: hco_core => null()
183
184 istamp = exdb('diagHBHt','DEBUT','NON')
185
186 call ver_printNameAndVersion('diagHBHt','RANDOMIZED DIAGNOSTIC of HBHt')
187
188 ! MPI initilization
189 call mmpi_initialize
190
191 call tmg_init(mmpi_myid, 'TMG_INFO')
192
193 call utl_tmg_start(0,'Main')
194
195 varMode='analysis'
196
197 call ram_setup
198
199 ! Do initial set up
200
201 obsMpiStrategy = 'LIKESPLITFILES'
202
203 call var_setup('VAR') ! obsColumnMode
204
205 ! Reading trials
206 call inn_getHcoVcoFromTrlmFile( hco_trl, vco_trl )
207 allocHeightSfc = ( vco_trl%Vcode /= 0 )
208
209 call gsv_allocate( stateVectorTrialHighRes, tim_nstepobs, hco_trl, vco_trl, &
210 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
211 mpi_distribution_opt='Tiles', dataKind_opt=4, &
212 allocHeightSfc_opt=allocHeightSfc, hInterpolateDegree_opt='LINEAR', &
213 beSilent_opt=.false. )
214 call gsv_zero( stateVectorTrialHighRes )
215 call gio_readTrials( stateVectorTrialHighRes )
216
217 ! Horizontally interpolate trials to trial columns
218 call inn_setupColumnsOnTrlLev( columnTrlOnTrlLev, obsSpaceData, hco_core, &
219 stateVectorTrialHighRes )
220
221 ! Interpolate trial columns to analysis levels and setup for linearized H
222 call inn_setupColumnsOnAnlIncLev(columnTrlOnTrlLev,columnTrlOnAnlIncLev)
223
224 ! Compute observation innovations and prepare obsSpaceData for minimization
225 call inn_computeInnovation(columnTrlOnTrlLev,obsSpaceData)
226
227 ! Compute perturbed
228 call diagHBHt(columnTrlOnAnlIncLev,obsSpaceData)
229
230 ! Deallocate memory related to B matrices
231 call bmat_finalize()
232
233 ! Now write out the observation data files
234 if ( .not. obsf_filesSplit() ) then
235 write(*,*) 'We read/write global observation files'
236 call obs_expandToMpiGlobal(obsSpaceData)
237 if (mmpi_myid == 0) call obsf_writeFiles(obsSpaceData)
238 else
239 ! redistribute obs data to how it was just after reading the files
240 call obs_MpiRedistribute(obsSpaceData,OBS_IPF)
241 call obsf_writeFiles(obsSpaceData)
242 end if
243
244
245 ! Deallocate copied obsSpaceData
246 call obs_finalize(obsSpaceData)
247
248 ! 3. Job termination
249
250 istamp = exfin('diagHBHt','FIN','NON')
251
252 call utl_tmg_stop(0)
253
254 call tmg_terminate(mmpi_myid, 'TMG_INFO')
255
256 call rpn_comm_finalize(ierr)
257
258contains
259
260 !--------------------------------------------------------------------------
261 ! var_setup
262 !--------------------------------------------------------------------------
263 subroutine var_setup(obsColumnMode)
264 !
265 !:Purpose: Control of the preprocessing of the variational assimilation
266 !
267 implicit none
268
269 ! Arguments:
270 character (len=*), intent(in) :: obsColumnMode
271
272 ! Locals:
273 integer :: dateStampFromObs
274 type(struct_vco),pointer :: vco_anl => null()
275 integer :: get_max_rss
276
277 write(*,*)
278 write(*,*) '-----------------------------------'
279 write(*,*) '-- Starting subroutine var_setup --'
280 write(*,*) '-----------------------------------'
281
282 !
283 !- Initialize the Temporal grid and set dateStamp from env variable
284 !
285 call tim_setup()
286
287 !
288 !- Initialize burp file names and set datestamp if not already
289 !
290 call obsf_setup(dateStampFromObs, 'analysis')
291 if (tim_getDateStamp() == 0) then
292 if (dateStampFromObs > 0) then
293 call tim_setDatestamp(dateStampFromObs)
294 else
295 call utl_abort('var_setup: dateStamp was not set')
296 end if
297 end if
298
299 !
300 !- Initialize constants
301 !
302 if ( mmpi_myid == 0 ) then
303 call mpc_printConstants(6)
304 call pre_printPrecisions
305 end if
306
307 !
308 !- Initialize variables of the model states
309 !
310 call gsv_setup
311 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
312
313 !
314 !- Initialize the Analysis grid
315 !
316 if(mmpi_myid.eq.0) write(*,*)''
317 if(mmpi_myid.eq.0) write(*,*)' preproc: Set hco parameters for analysis grid'
318 call hco_SetupFromFile(hco_anl, './analysisgrid', 'ANALYSIS', 'Analysis' ) ! IN
319
320 if ( hco_anl % global ) then
321 hco_core => hco_anl
322 else
323 !- Iniatilized the core (Non-Exteded) analysis grid
324 call hco_SetupFromFile( hco_core, './analysisgrid', 'COREGRID', 'AnalysisCore' ) ! IN
325 end if
326
327 !
328 !- Initialisation of the analysis grid vertical coordinate from analysisgrid file
329 !
330 call vco_SetupFromFile( vco_anl, & ! OUT
331 './analysisgrid') ! IN
332
333 call col_setVco(columnTrlOnAnlIncLev,vco_anl)
334 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
335
336 !
337 !- Setup and read observations
338 !
339 call inn_setupObs(obsSpaceData, hco_anl, obsColumnMode, obsMpiStrategy, varMode) ! IN
340 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
341
342 !
343 !- Basic setup of columnData module
344 !
345 call col_setup
346 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
347
348 !
349 !- Memory allocation for background column data
350 !
351 call col_allocate(columnTrlOnAnlIncLev,obs_numheader(obsSpaceData),mpiLocal_opt=.true.)
352
353 !
354 !- Initialize the observation error covariances
355 !
356 call oer_setObsErrors(obsSpaceData, varMode) ! IN
357 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
358
359 !
360 !- Initialize the background-error covariance, also sets up control vector module (cvm)
361 !
362 call bmat_setup(hco_anl,hco_core,vco_anl)
363 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
364
365 !
366 ! - Initialize the gridded variable transform module
367 !
368
369 call gvt_setup(hco_anl,hco_core,vco_anl)
370 call gvt_setupRefFromTrialFiles('HU')
371 call gvt_setupRefFromTrialFiles('height')
372
373 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
374
375 end subroutine var_setup
376
377 !--------------------------------------------------------------------------
378 ! diagHBHt
379 !--------------------------------------------------------------------------
380 subroutine diagHBHt(columnTrlOnAnlIncLev, obsSpaceData)
381 !
382 !:Purpose: Main calculations of HBHt
383 !
384 implicit none
385
386 ! Arguments:
387 type(struct_obs), intent(inout) :: obsSpaceData ! Observation-related data
388 type(struct_columnData), intent(inout) :: columnTrlOnAnlIncLev ! Background columns interpolated to anl levels and obs horiz locations
389
390 ! Locals:
391 type(struct_columnData) :: columnAnlInc
392 type(struct_gsv) :: statevector
393 type(struct_vco), pointer :: vco_anl
394 real(8) ,allocatable :: random_vector(:)
395 real(8) ,allocatable :: local_random_vector(:)
396 integer :: index_body, local_dimension, jj, ierr, dateprnt, timeprnt, nrandseed, istat
397 integer ,external :: newdate, get_max_rss
398 real(8) ,external :: gasdev
399
400 !
401 !- 1. Initialization
402
403 write(*,*)
404 write(*,*) 'Computing perturbations for randomized HBHT evaluation START'
405
406 vco_anl => col_getVco(columnTrlOnAnlIncLev)
407 !- 1.3 Create a gridstatevector to store the perturbations
408 call gsv_allocate(statevector,tim_nstepobsinc,hco_anl,vco_anl, &
409 dataKind_opt=pre_incrReal,mpi_local_opt=.true.)
410
411 !- 1.4 Create column vectors to store the perturbation interpolated to obs horizontal locations
412 call col_setVco(columnAnlInc,vco_anl)
413 call col_allocate(columnAnlInc,col_getNumCol(columnTrlOnAnlIncLev),mpiLocal_opt=.true.)
414
415 !- 1.6
416 call oti_timeBinning(obsSpaceData,tim_nstepobsinc)
417
418 !
419 !- 2. Compute the perturbations
420 !
421
422 !- 2.1 Random perturbations
423 write(*,*)
424 write(*,*) 'Generating random perturbation:'
425
426 !- Global vector (same for each processors)
427 allocate(random_vector(cvm_nvadim_mpiglobal),stat =istat )
428 allocate(local_random_vector(cvm_nvadim),stat =istat )
429
430 !- Initialize random number generator
431 ierr = newdate(tim_getDatestamp(), dateprnt, timeprnt, -3)
432 nrandseed=100*dateprnt + int(timeprnt/100.0)
433 write(*,*) 'diagHBHt: Random seed set to ',nrandseed
434 call rng_setup(nrandseed)
435 ! Generate a random vector from N(0,1)
436 do jj = 1, cvm_nvadim_mpiglobal
437 random_vector(jj) = rng_gaussian()
438 enddo
439 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
440 !- Extract only the subvector for this processor
441 call bmat_reduceToMPILocal(local_random_vector, & ! OUT
442 random_vector) ! IN
443 local_dimension = size(local_random_vector)
444 !- Transform to control variables in physical space
445 call bmat_sqrtB(local_random_vector,local_dimension,statevector)
446 !- 2.2 Interpolation to the observation horizontal locations
447
448 call s2c_tl( statevector, & ! IN
449 columnAnlInc, & ! OUT (H_horiz EnsPert)
450 columnTrlOnAnlIncLev, obsSpaceData ) ! IN
451 !- 2.3 Interpolation to observation space
452 call oop_Htl(columnAnlInc,columnTrlOnAnlIncLev,obsSpaceData,min_nsim=1)
453
454 !- Copy from OBS_WORK to OBS_HPHT
455
456 do index_body = 1, obs_numBody(obsSpaceData)
457 call obs_bodySet_r(obsSpaceData,OBS_HPHT,index_body, &
458 obs_bodyElem_r(obsSpaceData,OBS_WORK,index_body) )
459 end do
460 !
461 !- 3. Ending/Deallocation
462 !
463 deallocate(random_vector, stat=istat)
464 deallocate(local_random_vector, stat=istat)
465 call col_deallocate(columnAnlInc)
466 call gsv_deallocate(statevector)
467 write(*,*)
468 write(*,*) 'Computing perturbations for randomized HBHT evaluation END'
469
470 end subroutine diagHBHt
471
472end program midas_diagHBHt
