1program midas_var
2 !
3 !:Purpose: Main program for performing data assimilation using one of the
4 ! following algorithms based on the incremental variational
5 ! approach:
6 !
7 ! - 3D-Var
8 ! - 3D- or 4D-EnVar
9 !
10 ! ---
11 !
12 !:Algorithm: The incremental variational data assimilation approach uses
13 ! observations to compute a correction to the background state
14 ! (i.e. the analysis increment) while taking into account the
15 ! specified uncertainties for both the observations and the
16 ! background state (i.e. the R and B covariance matrices,
17 ! respectively). The analysis increment is computed by finding the
18 ! minimum value of a cost function. This cost function is a scalar
19 ! measure of the weighted departure of the corrected state from
20 ! both the background state and the observations. The minimization
21 ! is perfomed with a standard minimization algorithm (currently
22 ! the quasi-Newton algorithm) that employs the gradient of the
23 ! cost function to enable the minimum to be found (or at least
24 ! approximated with sufficient accuracy) with relatively few
25 ! iterations. The analysis increment is not assumed to be on the
26 ! same horizontal and vertical grid as the background state. In
27 ! addition, the temporal resolution of the background state and
28 ! analysis increment are not assumed to be equal.
29 !
30 ! --
31 !
32 ! In practical terms, the background state is used to compute the
33 ! "innovation" vector, that is, the difference between the
34 ! observations and the background state in observation space:
35 ! ``y-H(xb)``. The function ``H()`` represents the non-linear
36 ! observation operators that map a gridded state vector into
37 ! observation space. The innovation vector is used in the cost
38 ! function calculation. Versions of the observation operators that
39 ! are linearized with respect to the background state (or updated
40 ! background state when using the outer loop) are also used in the
41 ! cost function calculation to transform the analysis increment
42 ! into observation space. The gradient of the cost function is
43 ! computed using the adjoint of these linearized operators.
44 !
45 ! --
46 !
47 ! After the minimization, the analysis increment is spatially
48 ! interpolated to the same grid as the background state and then
49 ! added to this state to obtain the analysis. For some variables a
50 ! physical constraint is imposed on the analysis
51 ! (e.g. non-negative humidity or sea-ice concentration between 0
52 ! and 1) and then the analysis increment is recomputed.
53 !
54 ! --
55 !
56 !:File I/O: The required input files and produced output files can vary
57 ! according to the application. Below are tables of files for
58 ! typical NWP 4D-EnVar (e.g. GDPS) and sea ice or SST 3D-Var
59 ! applications.
60 !
61 ! --
62 !
63 !============================================== ==============================================================
64 ! Input and Output Files (NWP applicaton) Description of file
65 !============================================== ==============================================================
66 ! ``flnml`` In - Main namelist file with parameters user may modify
67 ! ``flnml_static`` In - The "static" namelist that should not be modified
68 ! ``trlm_$NN`` (e.g. ``trlm_01``) In - Background state (a.k.a. trial) files for each timestep
69 ! ``analysisgrid`` In - File defining grid for computing the analysis increment
70 ! ``bgcov`` In - Static (i.e. NMC) B matrix file for NWP fields
71 ! ``bgchemcov`` In - Static B matrix file for chemistry fields
72 ! ``ensemble/$YYYYMMDDHH_006_$NNNN`` In - Ensemble member files defining ensemble B matrix
73 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
74 ! ``obserr`` In - Observation error statistics
75 ! ``obsinfo_chm`` In - Something needed for chemistry assimilation?
76 ! ``preconin`` In - Preconditioning file (Hessian of the cost function)
77 ! ``rebm_$MMMm`` (e.g. ``rebm_180m``) Out - Analysis increment on the (low-res) analysis grid
78 ! ``rehm_$MMMm`` Out - Analysis increment on the (high-res) trial grid
79 ! ``anlm_$MMMm`` Out - Analysis on the trial grid
80 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - Updated obs file for each "family" and MPI task
81 ! Remainder are files related to radiance obs:
82 ! ``stats_$SENSOR_assim`` In - Satellite radiance observation errors of difference sensors
83 ! ``stats_tovs`` In - Satellite radiance observation errors
84 ! ``stats_tovs_symmetricObsErr`` In - User-defined symmetric TOVS errors for all sky
85 ! ``Cmat_$PLATFORM_$SENSOR.dat`` In - Inter-channel observation-error correlations
86 ! ``ceres_global.std`` In - High-res surface type and water fraction for radiance obs
87 ! ``rtcoef_$PLATFORM_$SENSOR.dat`` In - RTTOV coefficient files
88 ! ``rttov_h2o_limits.dat`` In - Min/max humidity limits applied to analysis
89 ! ``ozoneclim98`` In - Ozone climatology
90 !============================================== ==============================================================
91 !
92 ! --
93 !
94 !============================================== ==============================================================
95 ! Input and Output Files (SST or Sea ice) Description of file
96 !============================================== ==============================================================
97 ! ``flnml`` In - Main namelist file with parameters user may modify
98 ! ``trlm_$NN`` (e.g. ``trlm_01``) In - Background state (a.k.a. trial) files for each timestep
99 ! ``analysisgrid`` In - File defining grid for computing the analysis increment
100 ! ``sea_ice_obs-err`` In - Observation error statistics (sea ice only)
101 ! ``bgstddev`` In - Background error stddev
102 ! ``diffusmod.std`` In - Normalization factor for diffusion operator correlations
103 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
104 ! ``rebm_$MMMm`` (e.g. ``rebm_180m``) Out - Analysis increment on the (low-res) analysis grid
105 ! ``rehm_$MMMm`` Out - Analysis increment on the (high-res) trial grid
106 ! ``anlm_$MMMm`` Out - Analysis on the trial grid
107 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - Updated obs file for each "family" and MPI task
108 !============================================== ==============================================================
109 !
110 ! --
111 !
112 !:Synopsis: Below is a summary of the ``var`` program calling sequence:
113 !
114 ! - **Initial setups:**
115 !
116 ! - Setup horizontal and vertical grid objects for "analysis
117 ! grid" from ``analysisgrid`` file and for "trial grid" from
118 ! first trial file: ``trlm_01``.
119 !
120 ! - Setup ``obsSpaceData`` object and read observations from
121 ! files: ``inn_setupObs``.
122 !
123 ! - Setup ``columnData`` and ``gridStateVector`` modules (read
124 ! list of analysis variables from namelist) and allocate column
125 ! object for storing trial on analysis levels.
126 !
127 ! - Setup the observation error statistics in ``obsSpaceData``
128 ! object: ``oer_setObsErrors``.
129 !
130 ! - Allocate a stateVector object on the trial grid and then
131 ! read the trials: ``gio_readTrials``.
132 !
133 ! - Setup the B matrices: ``bmat_setup``.
134 !
135 ! - Setup the ``gridVariableTransforms`` and ``minimization``
136 ! modules.
137 !
138 ! - **Outer loop** (only 1 iteration when no outer loop is used)
139 ! during which the analysis increment is computed to correct the
140 ! current "updated state" which is used to compute the
141 ! innovation and serve as the starting point for the
142 ! minimization of the following outer loop iteration:
143 !
144 ! - Impose RTTOV humidity limits on updated state (initially the
145 ! trial) on trial grid: ``qlim_rttovLimit``.
146 !
147 ! - Use the updated state on trial grid to setup the reference
148 ! state used by ``gridVariableTransforms`` module for height
149 ! calculations.
150 !
151 ! - Compute ``columnTrlOnTrlLev`` and ``columnTrlOnAnlIncLev`` from
152 ! updated state: ``inn_setupColumnsOnTrlLev``,
153 ! ``inn_setupColumnsOnAnlIncLev``
154 !
155 ! - Compute innovation from updated state:
156 ! ``inn_computeInnovation``.
157 !
158 ! - Use the updated state on trial grid to setup the reference
159 ! state used by ``gridVariableTransforms`` module for ``LQ``
160 ! to ``HU`` calculations.
161 !
162 ! - Do the minimization for this outer loop iteration:
163 ! ``min_minimize`` to obtain ``controlVectorIncr``.
164 !
165 ! - Update sum of all computed increment control vectors:
166 ! ``controlVectorIncrSum(:)`` which is needed in the
167 ! background cost function for outer loop.
168 !
169 ! - Compute ``stateVectorIncr`` (on analysis grid) from
170 ! ``controlVectorIncr``: ``inc_getIncrement``.
171 !
172 ! - Interpolate and add ``stateVectorIncr`` to the updated
173 ! state: ``inc_computeHighResAnalysis``.
174 !
175 ! - If requested, impose saturation and RTTOV humidity limits on
176 ! the updated state.
177 !
178 ! - Write increment (or sum of increments when outer loop used)
179 ! to file: ``inc_writeIncrement``.
180 !
181 ! - End of outer loop.
182 !
183 ! - **Final steps**, after the outer loop:
184 !
185 ! - If requested, compute final cost function value using
186 ! non-linear observation operators.
187 !
188 ! - Write the final analysis and recomputed complete increment on
189 ! the trial grid: ``inc_writeincAndAnalHighRes``.
190 !
191 ! - Various final steps, including: write the Hessian to binary
192 ! file (``min_writeHessian``), update the observation files
193 ! (``obsf_writeFiles``).
194 !
195 ! --
196 !
197 !:Options: `List of namelist blocks <../namelists_in_each_program.html#var>`_
198 ! that can affect the ``var`` program.
199 !
200 ! * The use of an outer loop is controlled by the namelist block
201 ! ``&NAMVAR`` read by the ``var`` program.
202 !
203 ! * The choice of 3D-Var vs. EnVar is controlled by the weights given
204 ! to the climatological (i.e. static) and ensemble-based background
205 ! error covariance (i.e. B matrix) components. The weights for all
206 ! B matrix components are zero be default and can be set to a
207 ! nonzero value through the namelist variable ``SCALEFACTOR`` in
208 ! the namelist block for each corresponding fortran module.
209 !
210 ! * Either algorithm can be used in combination with "First guess at
211 ! appropriate time" (i.e. FGAT) which is controlled by the namelist
212 ! variable ``DSTEPOBS`` (in ``NAMTIME``) that specifies the length
213 ! of time (in hours) between times when the background state is
214 ! used to compute the innovation (i.e. O-B).
215 !
216 ! * Similarly, the choice between 3D-EnVar and 4D-EnVar is controlled
217 ! by the namelist variable ``DSTEPOBSINC`` (also in ``&NAMTIME``)
218 ! which specifies the length of time (in hours) between times when
219 ! the analysis increment is computed. In the context of EnVar, if
220 ! this is less than the assimilation time window, then the
221 ! ensembles used to construct the ensemble-based B matrix component
222 ! will be used at multiple times within the assimilation window to
223 ! obtain 4D covariances (i.e. 4D-EnVar).
224 !
225 ! * Some of the other relevant namelist blocks used to configure the
226 ! variational analysis are listed in the following table:
227 !
228 !======================== ============ ==============================================================
229 ! Module Namelist Description of what is controlled
230 !======================== ============ ==============================================================
231 ! ``minimization_mod`` ``NAMMIN`` maximum number of iterations, convergence criteria and many
232 ! additional parameters for controlling the minimization
233 ! ``timeCoord_mod`` ``NAMTIME`` assimilation time window length, temporal resolution of
234 ! the background state and increment
235 ! ``bMatrixEnsemble_mod`` ``NAMBEN`` weight and other parameters for ensemble-based B matrix
236 ! component
237 ! ``bMatrixHI_mod`` ``NAMBHI`` weight and other parameters for the climatological B matrix
238 ! component based on homogeneous-isotropic covariances
239 ! represented in spectral space
240 ! Other B matrix modules various weight and other parameters for each type of B matrix
241 !======================== ============ ==============================================================
242 !
243 use version_mod
244 use codePrecision_mod
245 use ramDisk_mod
246 use utilities_mod
247 use midasMpi_mod
248 use message_mod
249 use mathPhysConstants_mod
250 use horizontalCoord_mod
251 use verticalCoord_mod
252 use humidityLimits_mod
253 use timeCoord_mod
254 use obsSpaceData_mod
255 use columnData_mod
256 use gridStateVector_mod
257 use gridStateVectorFileIO_mod
258 use obsSpaceDiag_mod
259 use controlVector_mod
260 use obsFiles_mod
261 use minimization_mod
262 use innovation_mod
263 use bMatrix_mod
264 use rMatrix_mod
265 use obsErrors_mod
266 use gridVariableTransforms_mod
267 use increment_mod
268 use biasCorrectionSat_mod
269 use varQC_mod
270 use tovsNL_mod
271 use stateToColumn_mod
272
273 implicit none
274
275 integer :: istamp, exdb, exfin
276 integer :: ierr, dateStampFromObs, nulnam
277 integer :: fclos, fnom
278 character(len=9) :: clmsg
279 character(len=48) :: obsMpiStrategy, varMode
280 real(8), allocatable :: controlVectorIncr(:)
281 real(8), allocatable :: controlVectorIncrSum(:)
282
283 type(struct_obs) , target :: obsSpaceData
284 type(struct_columnData), target :: columnTrlOnAnlIncLev
285 type(struct_columnData), target :: columnTrlOnTrlLev
286 type(struct_gsv) :: stateVectorIncr
287 type(struct_gsv) :: stateVectorIncrSum
288 type(struct_gsv) :: stateVectorUpdateHighRes
289 type(struct_gsv) :: stateVectorTrial
290 type(struct_gsv) :: stateVectorPsfcHighRes
291 type(struct_gsv) :: stateVectorPsfc
292 type(struct_gsv) :: stateVectorAnal
293 type(struct_hco) , pointer :: hco_anl => null()
294 type(struct_vco) , pointer :: vco_anl => null()
295 type(struct_hco) , pointer :: hco_trl => null()
296 type(struct_vco) , pointer :: vco_trl => null()
297 type(struct_hco) , pointer :: hco_core => null()
298
299 integer :: outerLoopIndex, numIterMaxInnerLoopUsed
300 integer :: numIterWithoutVarqc, numInnerLoopIterDone
301
302 logical :: allocHeightSfc, applyLimitOnHU
303 logical :: deallocHessian, isMinimizationFinalCall
304 logical :: varqcActive, applyVarqcOnNlJo
305 logical :: filterObsAndInitOer
306 logical :: deallocInterpInfoNL
307
308 integer, parameter :: maxNumOuterLoopIter = 15
309
310 ! namelist variables
311 integer :: numOuterLoopIterations ! number of outer loop iterations (default=1)
312 integer :: numIterMaxInnerLoop(maxNumOuterLoopIter) ! number of each inner loop iterations
313 logical :: limitHuInOuterLoop ! impose humidity limits on each outer loop iteration
314 logical :: computeFinalNlJo ! compute final cost function using non-linear H()
315 NAMELIST /NAMVAR/ numOuterLoopIterations, numIterMaxInnerLoop, limitHuInOuterLoop
316 NAMELIST /NAMVAR/ computeFinalNlJo
317
318 istamp = exdb('VAR','DEBUT','NON')
319
320 call ver_printNameAndVersion('var','Variational Assimilation')
321
322 ! MPI initialization
323 call mmpi_initialize
324
325 call tmg_init(mmpi_myid, 'TMG_INFO')
326
327 call utl_tmg_start(0,'Main')
328
329 if (mmpi_myid == 0) then
330 clmsg = 'VAR3D_BEG'
331 call utl_writeStatus(clmsg)
332 end if
333
334 varMode='analysis'
335
336 ! Setup the ram disk
337 call ram_setup
338
339 ! Do initial set up
340
341 ! Set/Read values for the namelist NAMVAR
342 ! Setting default namelist variable values
343 numOuterLoopIterations = 1
344 limitHuInOuterLoop = .false.
345 numIterMaxInnerLoop(:) = 0
346 computeFinalNlJo = .false.
347
348 if ( .not. utl_isNamelistPresent('NAMVAR','./flnml') ) then
349 call msg('midas-var','namvar is missing in the namelist. '&
350 //'The default values will be taken.', mpiAll_opt=.false.)
351
352 else
353 ! read in the namelist NAMVAR
354 nulnam = 0
355 ierr = fnom(nulnam, './flnml', 'FTN+SEQ+R/O', 0)
356 read(nulnam, nml=namvar, iostat=ierr)
357 if( ierr /= 0) call utl_abort('midas-var: Error reading namelist')
358 ierr = fclos(nulnam)
359 end if
360 if ( mmpi_myid == 0 ) write(*,nml=namvar)
361
362 if ( numOuterLoopIterations > maxNumOuterLoopIter ) then
363 call utl_abort('midas-var: numOuterLoopIterations is greater than max value')
364 end if
365
366 if ( numOuterLoopIterations > 1 .and. &
367 .not. all(numIterMaxInnerLoop(1:numOuterLoopIterations) > 0) ) then
368 call utl_abort('midas-var: some numIterMaxInnerLoop(:) in namelist are negative or zero')
369 end if
370
371 obsMpiStrategy = 'LIKESPLITFILES'
372
373 ! Initialize the Temporal grid and set dateStamp from env variable
374 call tim_setup()
375
376 ! Initialize observation file names and set datestamp if not already
377 call obsf_setup(dateStampFromObs, varMode)
378 if (tim_getDateStamp() == 0) then
379 if (dateStampFromObs > 0) then
380 call tim_setDatestamp(datestampFromObs)
381 else
382 call utl_abort('var_setup: DateStamp was not set')
383 end if
384 end if
385
386 ! Initialize constants
387 if ( mmpi_myid == 0 ) then
388 call mpc_printConstants(6)
389 call pre_printPrecisions
390 end if
391
392 ! Initialize the Analysis grid
393 call msg('midas-var','Set hco parameters for analysis grid', mpiAll_opt=.false.)
394 call hco_SetupFromFile(hco_anl, './analysisgrid', 'ANALYSIS', 'Analysis' ) ! IN
395
396 if ( hco_anl % global ) then
397 hco_core => hco_anl
398 else
399 ! Initialize the core (Non-Extended) analysis grid
400 call msg('midas-var','Set hco parameters for core grid', mpiAll_opt=.false.)
401 call hco_SetupFromFile( hco_core, './analysisgrid', 'COREGRID', 'AnalysisCore' ) ! IN
402 end if
403
404 ! Initialisation of the analysis grid vertical coordinate from analysisgrid file
405 call vco_SetupFromFile( vco_anl, & ! OUT
406 './analysisgrid') ! IN
407
408 call col_setVco(columnTrlOnAnlIncLev,vco_anl)
409 call msg_memUsage('var')
410
411 ! Setup and read observations
412 call inn_setupObs(obsSpaceData, hco_anl, 'VAR', obsMpiStrategy, varMode) ! IN
413 call msg_memUsage('var')
414
415 ! Basic setup of columnData module
416 call col_setup
417 call msg_memUsage('var')
418
419 !- Memory allocation for background column data
420 call col_allocate(columnTrlOnAnlIncLev,obs_numheader(obsSpaceData),mpiLocal_opt=.true.)
421
422 ! Initialize the observation error covariances
423 call oer_setObsErrors(obsSpaceData, varMode) ! IN
424 call msg_memUsage('var')
425
426 ! Initialize list of analyzed variables.
427 call gsv_setup
428 call msg_memUsage('var')
429
430 ! Reading trials
431 call inn_getHcoVcoFromTrlmFile( hco_trl, vco_trl )
432 allocHeightSfc = ( vco_trl%Vcode /= 0 )
433
434 call gsv_allocate( stateVectorUpdateHighRes, tim_nstepobs, hco_trl, vco_trl, &
435 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
436 mpi_distribution_opt='Tiles', dataKind_opt=pre_incrReal, &
437 allocHeightSfc_opt=allocHeightSfc, hInterpolateDegree_opt='LINEAR')
438 call gsv_zero( stateVectorUpdateHighRes )
439 call gio_readTrials( stateVectorUpdateHighRes )
440 call msg_memUsage('var')
441
442 ! Initialize the background-error covariance, also sets up control vector module (cvm)
443 call bmat_setup(hco_anl,hco_core,vco_anl)
444 call msg_memUsage('var')
445
446 ! Initialize the gridded variable transform module
447 call gvt_setup(hco_anl,hco_core,vco_anl)
448
449 ! Set up the minimization module, now that the required parameters are known
450 ! NOTE: some global variables remain in minimization_mod that must be initialized before
451 ! inn_setupColumnsOnTrlLev
452 call min_setup( cvm_nvadim, hco_anl, & ! IN
453 varqc_opt=varqcActive, nwoqcv_opt=numIterWithoutVarqc ) ! OUT
454 allocate(controlVectorIncr(cvm_nvadim),stat=ierr)
455 if (ierr /= 0) then
456 call msg('var','Problem allocating memory for controlVectorIncr'//str(ierr))
457 call utl_abort('aborting in VAR')
458 end if
459 call utl_reallocate(controlVectorIncrSum,cvm_nvadim)
460 call msg_memUsage('var')
461
462 numInnerLoopIterDone = 0
463
464 ! Enter outer-loop
465 outer_loop: do outerLoopIndex = 1, numOuterLoopIterations
466 call msg('var','start of outer-loop index='//str(outerLoopIndex))
467
468 ! Impose limits on ALL cloud variables
469 call qlim_rttovLimit(stateVectorUpdateHighRes, applyLimitToCloud_opt=.true.)
470
471 ! Initialize stateVectorRefHeight for transforming TT/HU/P0 increments to
472 ! height/pressure increments.
473 if ( (gsv_varExist(stateVectorUpdateHighRes,'P_T') .and. &
474 gsv_varExist(stateVectorUpdateHighRes,'P_M')) .or. &
475 (gsv_varExist(stateVectorUpdateHighRes,'Z_T') .and. &
476 gsv_varExist(stateVectorUpdateHighRes,'Z_M')) ) then
477
478 call gvt_setupRefFromStateVector( stateVectorUpdateHighRes, 'height' )
479
480 call msg_memUsage('var')
481 end if
482
483 ! Horizontally interpolate high-resolution stateVectorUpdate to trial columns
484 deallocInterpInfoNL = ( numOuterLoopIterations <= 1 )
485 call inn_setupColumnsOnTrlLev( columnTrlOnTrlLev, obsSpaceData, hco_core, &
486 stateVectorUpdateHighRes, &
487 deallocInterpInfoNL_opt=deallocInterpInfoNL )
488 call msg_memUsage('var')
489
490 ! Interpolate trial columns to analysis levels and setup for linearized H
491 call inn_setupColumnsOnAnlIncLev( columnTrlOnTrlLev, columnTrlOnAnlIncLev )
492 call msg_memUsage('var')
493
494 ! Determine if to apply varqc to Jo of non-linear operator
495 applyVarqcOnNlJo = ( varqcActive .and. numInnerLoopIterDone > numIterWithoutVarqc )
496 if ( applyVarqcOnNlJo ) call msg('var','applying varqc to non-linear Jo', mpiAll_opt=.false.)
497
498 ! Compute observation innovations and prepare obsSpaceData for minimization
499 filterObsAndInitOer = ( outerLoopIndex == 1 )
500 call inn_computeInnovation( columnTrlOnTrlLev, obsSpaceData, &
501 filterObsAndInitOer_opt=filterObsAndInitOer, &
502 applyVarqcOnNlJo_opt=applyVarqcOnNlJo, &
503 callSetErrGpsgb_opt=filterObsAndInitOer )
504 call msg_memUsage('var')
505
506 ! Initialize stateVectorRefHU for doing variable transformation of the increments.
507 if ( gsv_varExist(stateVectorUpdateHighRes,'HU') ) then
508 applyLimitOnHU = ( limitHuInOuterLoop .and. outerLoopIndex > 1 )
509
510 call gvt_setupRefFromStateVector( stateVectorUpdateHighRes, 'HU', &
511 applyLimitOnHU_opt=applyLimitOnHU )
512
513 call msg_memUsage('var')
514 end if
515
516 ! Do minimization of cost function. Use numIterMaxInnerLoop from NAMVAR, instead of
517 ! nitermax from NAMMIN, when numOuterLoopIterations > 1
518 controlVectorIncr(:) = 0.0d0
519 deallocHessian = ( numOuterLoopIterations == 1 )
520 isMinimizationFinalCall = ( outerLoopIndex == numOuterLoopIterations )
521 call min_minimize( outerLoopIndex, columnTrlOnAnlIncLev, obsSpaceData, controlVectorIncrSum, &
522 controlVectorIncr, numIterMaxInnerLoop(outerLoopIndex), &
523 deallocHessian_opt=deallocHessian, &
524 isMinimizationFinalCall_opt=isMinimizationFinalCall, &
525 numIterMaxInnerLoopUsed_opt=numIterMaxInnerLoopUsed )
526 numInnerLoopIterDone = numInnerLoopIterDone + numIterMaxInnerLoopUsed
527 call msg_memUsage('var')
528
529 ! Accumulate control vector increments of all the previous iterations
530 controlVectorIncrSum(:) = controlVectorIncrSum(:) + controlVectorIncr(:)
531
532 ! Compute satellite bias correction increment and write to file on last outer-loop
533 ! iteration
534 if ( outerLoopIndex == numOuterLoopIterations ) then
535 call bcs_writebias(controlVectorIncr)
536 end if
537
538 call tvs_deallocateProfilesNlTlAd
539
540 call gsv_allocate(stateVectorIncr, tim_nstepobsinc, hco_anl, vco_anl, &
541 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
542 dataKind_opt=pre_incrReal, allocHeight_opt=.false., allocPressure_opt=.false.)
543
544 ! get final increment with mask if it exists
545 call inc_getIncrement( controlVectorIncr, stateVectorIncr, cvm_nvadim )
546 call gio_readMaskFromFile( stateVectorIncr, './analysisgrid' )
547 call msg_memUsage('var')
548
549 ! Compute high-resolution analysis on trial grid
550 call inc_computeHighResAnalysis( stateVectorIncr, & ! IN
551 stateVectorUpdateHighRes, stateVectorPsfcHighRes ) ! OUT
552 call msg_memUsage('var')
553
554 ! Impose limits on stateVectorUpdateHighRes only when outer loop is used.
555 if ( limitHuInOuterLoop ) then
556 call msg('var','impose limits on stateVectorUpdateHighRes')
557 call qlim_saturationLimit( stateVectorUpdateHighRes )
558 call qlim_rttovLimit( stateVectorUpdateHighRes )
559 end if
560
561 ! prepare to write incremnt when no outer-loop, or sum of increments at last
562 ! outer-loop iteration.
563 if ( numOuterLoopIterations > 1 .and. &
564 outerLoopIndex == numOuterLoopIterations ) then
565
566 call gsv_allocate(stateVectorIncrSum, tim_nstepobsinc, hco_anl, vco_anl, &
567 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
568 dataKind_opt=pre_incrReal, allocHeight_opt=.false., allocPressure_opt=.false.)
569 call inc_getIncrement( controlVectorIncrSum, stateVectorIncrSum, cvm_nvadim )
570 call gio_readMaskFromFile( stateVectorIncrSum, './analysisgrid' )
571 call msg_memUsage('var')
572
573 call inc_writeIncrement( stateVectorIncrSum, & ! IN
574 ip3ForWriteToFile_opt=0 ) ! IN
575 call gsv_deallocate( stateVectorIncrSum )
576 else if ( numOuterLoopIterations == 1 ) then
577 call inc_writeIncrement( stateVectorIncr, & ! IN
578 ip3ForWriteToFile_opt=0 ) ! IN
579 end if
580 call msg_memUsage('var')
581
582 call gsv_deallocate( stateVectorIncr )
583
584 call msg('var','end of outer-loop index='//str(outerLoopIndex))
585 end do outer_loop
586
587 ! Set the QC flags to be consistent with VAR-QC if control analysis
588 if ( varqcActive ) call vqc_listrej(obsSpaceData)
589
590 if ( computeFinalNlJo ) then
591 ! Horizontally interpolate high-resolution stateVectorUpdate to trial columns
592 call inn_setupColumnsOnTrlLev( columnTrlOnTrlLev, obsSpaceData, hco_core, &
593 stateVectorUpdateHighRes )
594 call msg_memUsage('var')
595
596 ! Determine if to apply varqc to Jo of non-linear operator
597 applyVarqcOnNlJo = ( varqcActive .and. numInnerLoopIterDone > numIterWithoutVarqc )
598 if ( applyVarqcOnNlJo ) call msg('var','applying varqc to non-linear Jo', mpiAll_opt=.false.)
599
600 ! Compute observation innovations and prepare obsSpaceData for minimization
601 filterObsAndInitOer = .false.
602 call inn_computeInnovation( columnTrlOnTrlLev, obsSpaceData, &
603 destObsColumn_opt=OBS_OMA, &
604 filterObsAndInitOer_opt=filterObsAndInitOer, &
605 applyVarqcOnNlJo_opt=applyVarqcOnNlJo , &
606 callSetErrGpsgb_opt=filterObsAndInitOer )
607 call msg_memUsage('var')
608 end if
609
610 ! Memory deallocations for non diagonal R matrices for radiances
611 call rmat_cleanup()
612
613 ! Conduct obs-space post-processing diagnostic tasks (some diagnostic
614 ! computations controlled by NAMOSD namelist in flnml)
615 call osd_ObsSpaceDiag( obsSpaceData, columnTrlOnAnlIncLev, hco_anl )
616
617 ! Deallocate memory related to B matrices and update stateVector
618 call bmat_finalize()
619
620 ! Deallocate structures needed for interpolation
621 call s2c_deallocInterpInfo( inputStateVectorType='nl' )
622 call s2c_deallocInterpInfo( inputStateVectorType='tlad' )
623
624 ! Post processing of analysis before writing (variable transform+humidity clipping)
625 call inc_analPostProcessing( stateVectorPsfcHighRes, stateVectorUpdateHighRes, & ! IN
626 stateVectorTrial, stateVectorPsfc, stateVectorAnal ) ! OUT
627 call gsv_deallocate( stateVectorUpdateHighRes )
628
629 ! compute and write the analysis (as well as the increment on the trial grid)
630 call inc_writeIncAndAnalHighRes( stateVectorTrial, stateVectorPsfc, &
631 stateVectorAnal )
632 call msg_memUsage('var')
633
634 if (mmpi_myid == 0) then
635 clmsg = 'REBM_DONE'
636 call utl_writeStatus(clmsg)
637 end if
638
639 ! write the Hessian
640 call min_writeHessian(controlVectorIncr)
641 deallocate(controlVectorIncr)
642
643 ! Deallocate memory related to variational bias correction
644 call bcs_finalize()
645
646 ! Now write out the observation data files
647 if (min_niter > 0) then
648 if ( .not. obsf_filesSplit() ) then
649 call msg('var','reading/writing global observation files')
650 call obs_expandToMpiGlobal(obsSpaceData)
651 if (mmpi_myid == 0) call obsf_writeFiles(obsSpaceData)
652 else
653 ! redistribute obs data to how it was just after reading the files
654 call obs_MpiRedistribute(obsSpaceData,OBS_IPF)
655 call obsf_writeFiles(obsSpaceData)
656 end if
657 end if
658
659 ! Deallocate copied obsSpaceData
660 call obs_finalize(obsSpaceData)
661
662 ! Job termination
663 istamp = exfin('VAR','FIN','NON')
664
665 if (mmpi_myid == 0) then
666 clmsg = 'VAR3D_END'
667 call utl_writeStatus(clmsg)
668 end if
669
670 call utl_tmg_stop(0)
671
672 call tmg_terminate(mmpi_myid, 'TMG_INFO')
673
674 call rpn_comm_finalize(ierr)
675
676end program midas_var
