1program midas_var1D
2 !
3 !:Purpose: Main program for one dimensional variational minimization
4 !
5 ! ---
6 !
7 !:Algorithm: This program performs a similar data assimilation procedure as the var program,
8 ! except without taking into account the horizontal or temporal dimensions.
9 ! The assimilation is performed separately at each horizontal location and time where
10 ! observations are present. The B matrix can be either an explicit representation of
11 ! the covariances produced by the program extractBmatrixFor1Dvar or from an ensemble
12 ! (controlled by the namelists). The resulting analysis and analysis increment at the
13 ! observation locations/times are output in standard files on a "Y" grid.
14 !
15 ! --
16 !
17 !:File I/O: The required input files and produced output files are listed as follows.
18 !
19 !============================================== ==============================================================
20 ! Input and Output Files (NWP application) Description of file
21 !============================================== ==============================================================
22 ! ``flnml`` In - Main namelist file with parameters user may modify
23 ! ``flnml_static`` In - The "static" namelist that should not be modified
24 ! ``trlm_$NN`` (e.g. ``trlm_01``) In - Background state (a.k.a. trial) files for each timestep
25 ! ``analysisgrid`` In - File defining grid for computing the analysis increment
26 ! ``Bmatrix_sea.bin`` In - 1DVar Bmatrix file over sea (output from extractBmatrixFor1DVar)
27 ! ``Bmatrix_land.bin`` In - 1DVar Bmatrix file over land (output from extractBmatrixFor1DVar)
28 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family" and MPI task
29 ! ``obserr`` In - Observation error statistics
30 ! ``rttov_h2o_limits.dat`` In - minimum and maximum humidity profile to clip analysis
31 ! ``pm1q`` In/Out - Preconditioning file (Hessian of the cost function)
32 ! ``rebm_$MMMm`` (e.g. ``rebm_180m``) Out - Analysis increment on Y grid
33 ! ``anlm_$MMMm`` Out - Analysis on Y grid
34 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - Updated obs file for each "family" and MPI task
35 ! Remainder are files related to radiance obs:
36 ! ``stats_tovs`` In - Observation error file for radiances
37 ! ``stats_tovs_symmetricObsErr`` In - user-defined symmetric TOVS errors for all sky
38 ! ``Cmat_$PLATFORM_$SENSOR.dat`` In - Inter-channel observation-error correlations
39 ! ``rtcoef_$PLATFORM_$SENSOR.H5`` In - RTTOV coefficient file HDF-5 format
40 ! ``rtcoef_$PLATFORM_$SENSOR.dat`` In - RTTOV coefficient file ASCII format
41 ! ``ozoneclim98`` In - ozone climatology standard file (Fortuin and Kelder)
42 !============================================== ==============================================================
43 !
44 ! --
45 !
46 !:Synopsis: Below is a summary of the ``var1D`` program calling sequence:
47 !
48 ! - **Initial setups:**
49 !
50 ! - Initialize temporal grid
51 !
52 ! - Initialize observation file names and set date stamp
53 !
54 ! - Setup horizontal and vertical grid objects for "analysis
55 ! grid" from ``analysisgrid`` file
56 !
57 ! - Setup ``obsSpaceData`` object and read observations from
58 ! files: ``inn_setupObs``
59 !
60 ! - Setup ``columnData`` module
61 !
62 ! - Setup the observation error statistics in ``obsSpaceData``
63 ! object: ``oer_setObsErrors``
64 !
65 ! - Setup the gridStateVector module (initialize list of analyzed variables)
66 !
67 ! - Get horizontal and vertical grid descriptors from trial fields: ``inn_getHcoVcoFromTrlmFile``,
68 ! and allocate a gridStateVector objects
69 !
70 ! - Read the trials: ``gio_readTrials``
71 !
72 ! - Setup the 1DVar B matrix: ``bmat1D_bsetup``
73 !
74 ! - Setup the ``gridVariableTransforms`` and ``minimization`` modules
75 !
76 ! - Horizontally interpolate high-resolution stateVectorUpdate to trial columns: ``inn_setupColumnsOnTrlLev``
77 !
78 ! - Interpolate trial columns to analysis levels and setup for linearized H: ``inn_setupColumnsOnAnlIncLev``
79 !
80 ! - Compute observation innovations and prepare obsSpaceData for minimization: ``inn_computeInnovation``
81 !
82 ! - **Minimization:**
83 !
84 ! - Do the minimization ``min_minimize`` to obtain ``controlVectorIncr``
85 !
86 ! - Get 1DVar increment from ``controlVectorIncr`` to ``columnAnlInc``: ``bmat1D_get1DvarIncrement``
87 ! - Transfer increment to ``stateVectorIncr``: ``var1D_transferColumnToYGrid``
88 !
89 ! - Write increment to file: ``inc_writeIncrement``
90 !
91 ! - **Final step:**
92 !
93 ! - Release resources
94 !
95 ! - Write hessian if requested to: ``min_writeHessian``
96 !
97 ! - Write updated observation files: ``obsf_writeFiles``
98 !
99 ! --
100 !
101 !:Options: `List of namelist blocks <../namelists_in_each_program.html#var1d>`_
102 ! that can affect the ``var1D`` program.
103 !
104 ! - The B matrix used by the ``var1D`` program is controlled by the namelist block
105 ! ``&nambmat1D`` of module bmatrix1DVar_mod
106 ! * scaleFactorHI scaling factors for HI variances
107 ! * scaleFactorHIHumidity scaling factors for HI humidity variances
108 ! * scaleFactorENs scaling factors for Ens variances
109 ! * scaleFactorEnsHumidity scaling factors for Ens humidity variances
110 ! * nEns number of ensemble members for the ensemble part of B matrix (negative means no ensemble)
111 ! * vLocalize vertical localization length scale for the ensemble part of B matrix
112 ! * includeAnlVar list of variable for the B matrix
113 ! * numIncludeAnlVar number of variables in includeAnlVar (to be removed later)
114 !
115 ! --
116 !
117 use version_mod
118 use codePrecision_mod
119 use ramDisk_mod
120 use utilities_mod
121 use midasMpi_mod
122 use mathPhysConstants_mod
123 use horizontalCoord_mod
124 use verticalCoord_mod
125 use timeCoord_mod
126 use obsSpaceData_mod
127 use columnData_mod
128 use gridStateVector_mod
129 use gridStateVectorFileIO_mod
130 use controlVector_mod
131 use obsFiles_mod
132 use minimization_mod
133 use innovation_mod
134 use obsErrors_mod
135 use gridVariableTransforms_mod
136 use increment_mod
137 use biasCorrectionSat_mod
138 use var1D_mod
139 use bMatrix1Dvar_mod
140
141 implicit none
142
143 integer :: istamp, exdb, exfin
144 integer :: ierr, dateStampFromObs
145 integer :: get_max_rss
146 character(len=48) :: obsMpiStrategy, varMode
147 real(8), allocatable :: controlVectorIncr(:)
148 real(8), allocatable :: controlVectorIncrSum(:)
149 type(struct_obs), target :: obsSpaceData
150 type(struct_columnData), target :: columnTrlOnAnlIncLev
151 type(struct_columnData), target :: columnTrlOnTrlLev
152 type(struct_columnData), target :: columnAnlInc
153 type(struct_gsv) :: stateVectorIncr
154 type(struct_gsv) :: stateVectorTrialHighRes
155 type(struct_gsv) :: stateVectorAnalysis
156 type(struct_hco), pointer :: hco_anl => null()
157 type(struct_vco), pointer :: vco_anl => null()
158 type(struct_hco), pointer :: hco_core => null()
159 type(struct_hco), pointer :: hco_trl => null()
160 type(struct_vco), pointer :: vco_trl => null()
161
162 integer :: outerLoopIndex, numIterMaxInnerLoop
163 logical :: allocHeightSfc
164
165 istamp = exdb('VAR1D', 'DEBUT', 'NON')
166
167 call ver_printNameAndVersion('var1D', '1D Variational Assimilation')
168
169 ! MPI initialization
170 call mmpi_initialize
171
172 call tmg_init(mmpi_myid, 'TMG_INFO')
173
174 call utl_tmg_start(0,'Main')
175
176 write(*,*)
177 write(*,*) 'Real Kind used for computing the increment =', pre_incrReal
178 write(*,*)
179
180 varMode='analysis'
181
182 ! Setup the ram disk
183 call ram_setup
184
185 ! Do initial set up
186
187 obsMpiStrategy = 'LIKESPLITFILES'
188
189 ! Initialize the Temporal grid and set dateStamp from env variable
190 call tim_setup
191
192 ! Initialize observation file names and set datestamp if not already
193 call obsf_setup(dateStampFromObs, varMode)
194 if (tim_getDateStamp() == 0) then
195 if (dateStampFromObs > 0) then
196 call tim_setDatestamp(dateStampFromObs)
197 else
198 call utl_abort('midas-var1D: DateStamp was not set')
199 end if
200 end if
201
202 ! Initialize constants
203 if (mmpi_myid == 0) call mpc_printConstants(6)
204
205 ! Initialize the Analysis grid
206 if (mmpi_myid == 0) write(*,*)
207 if (mmpi_myid == 0) write(*,*) 'midas-var1D: Set hco parameters for analysis grid'
208 call hco_SetupFromFile(hco_anl, './analysisgrid', 'ANALYSIS', 'Analysis' ) ! IN
209
210 if ( hco_anl % global ) then
211 hco_core => hco_anl
212 else
213 !- Initialize the core (Non-Extended) analysis grid
214 if (mmpi_myid == 0) write(*,*) 'midas-var1D: Set hco parameters for core grid'
215 call hco_SetupFromFile(hco_core, './analysisgrid', 'COREGRID', 'AnalysisCore' ) ! IN
216 end if
217
218 ! Initialisation of the analysis grid vertical coordinate from analysisgrid file
219 call vco_SetupFromFile(vco_anl, & ! OUT
220 './analysisgrid') ! IN
221
222 call col_setVco(columnTrlOnAnlIncLev, vco_anl)
223 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
224
225 ! Setup and read observations
226 call inn_setupObs(obsSpaceData, hco_anl, 'VAR', obsMpiStrategy, varMode) ! IN
227 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
228
229 ! Basic setup of columnData module
230 call col_setup
231 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
232
233 ! Memory allocation for background column data
234 call col_allocate(columnTrlOnAnlIncLev, obs_numheader(obsSpaceData), mpiLocal_opt=.true.)
235
236 ! Initialize the observation error covariances
237 call oer_setObsErrors(obsSpaceData, varMode) ! IN
238 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
239
240 ! Initialize list of analyzed variables.
241 call gsv_setup
242 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
243
244 ! Reading trials
245 call inn_getHcoVcoFromTrlmFile(hco_trl, vco_trl)
246 allocHeightSfc = ( vco_trl%Vcode /= 0 )
247
248 call gsv_allocate(stateVectorTrialHighRes, tim_nstepobs, hco_trl, vco_trl, &
249 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
250 mpi_distribution_opt='Tiles', dataKind_opt=pre_incrReal, &
251 allocHeightSfc_opt=allocHeightSfc, hInterpolateDegree_opt='LINEAR', &
252 beSilent_opt=.false.)
253 call gsv_zero(stateVectorTrialHighRes)
254 call gio_readTrials(stateVectorTrialHighRes)
255 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
256
257 ! Initialize the background-error covariance, also sets up control vector module (cvm)
258 call bmat1D_bsetup(vco_anl, hco_anl, obsSpaceData)
259 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
260
261 ! Initialize the gridded variable transform module
262 call gvt_setup(hco_anl, hco_core, vco_anl)
263
264 ! Set up the minimization module, now that the required parameters are known
265 ! NOTE: some global variables remain in minimization_mod that must be initialized before
266 ! inn_setupBackgroundColumns
267 call min_setup(cvm_nvadim, hco_anl, oneDVarMode_opt=.true.) ! IN
268 allocate(controlVectorIncr(cvm_nvadim),stat=ierr)
269 if (ierr /= 0) then
270 write(*,*) 'midas-var1D: Problem allocating memory for ''controlVectorIncr''',ierr
271 call utl_abort('aborting in VAR1D')
272 end if
273 call utl_reallocate(controlVectorIncrSum,cvm_nvadim)
274 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
275
276 ! Horizontally interpolate high-resolution stateVectorUpdate to trial columns
277 call inn_setupColumnsOnTrlLev(columnTrlOnTrlLev, obsSpaceData, hco_core, &
278 stateVectorTrialHighRes )
279
280 ! Interpolate trial columns to analysis levels and setup for linearized H
281 call inn_setupColumnsOnAnlIncLev( columnTrlOnTrlLev,columnTrlOnAnlIncLev )
282
283 ! Compute observation innovations and prepare obsSpaceData for minimization
284 call inn_computeInnovation(columnTrlOnTrlLev, obsSpaceData)
285
286 ! Do minimization of cost function
287 outerLoopIndex = 1
288 numIterMaxInnerLoop = 0
289 controlVectorIncr(:) = 0.0d0
290 call min_minimize( outerLoopIndex, columnTrlOnAnlIncLev, obsSpaceData, controlVectorIncrSum, &
291 controlVectorIncr, numIterMaxInnerLoop )
292 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
293
294 ! Compute satellite bias correction increment and write to file
295 ! Is is still necessary ? (will do nothing, but does it make sense in 1DVar mode ?)
296 call bcs_writebias(controlVectorIncr)
297
298 call col_setVco(columnAnlInc, col_getVco(columnTrlOnAnlIncLev))
299 call col_allocate(columnAnlInc, col_getNumCol(columnTrlOnAnlIncLev), mpiLocal_opt=.true.)
300 call col_zero(columnAnlInc)
301 ! get final increment
302 call bmat1D_get1DvarIncrement(controlVectorIncr, columnAnlInc, columnTrlOnAnlIncLev, obsSpaceData, cvm_nvadim)
303 call var1D_transferColumnToYGrid(stateVectorIncr, obsSpaceData, columnAnlInc, bmat1D_includeAnlVar)
304
305 ! output the analysis increment
306 call inc_writeIncrement(stateVectorIncr)
307 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
308
309 ! compute and write the analysis (as well as the increment on the trial grid)
310 call var1d_transferColumnToYGrid(stateVectorAnalysis, obsSpaceData, columnTrlOnAnlIncLev, bmat1D_includeAnlVar)
311
312 if (mmpi_myId == 0) call gsv_add(statevectorIncr, statevectorAnalysis)
313
314 call inc_writeAnalysis(stateVectorAnalysis)
315
316 ! Deallocate memory related to B matrices
317 call var1D_finalize()
318 if (mmpi_myId == 0) then
319 call gsv_deallocate(stateVectorIncr)
320 call gsv_deallocate(stateVectorAnalysis)
321 end if
322
323 ! write the Hessian
324 call min_writeHessian(controlVectorIncr)
325 deallocate(controlVectorIncrSum)
326 deallocate(controlVectorIncr)
327
328 ! Deallocate memory related to variational bias correction
329 call bcs_finalize()
330 ! Now write out the observation data files
331 if (min_niter > 0) then
332 if ( .not. obsf_filesSplit() ) then
333 write(*,*) 'We read/write global observation files'
334 call obs_expandToMpiGlobal(obsSpaceData)
335 if (mmpi_myid == 0) call obsf_writeFiles(obsSpaceData)
336 else
337 ! redistribute obs data to how it was just after reading the files
338 call obs_MpiRedistribute(obsSpaceData, OBS_IPF)
339 call obsf_writeFiles(obsSpaceData)
340 end if
341 end if
342
343 ! Deallocate copied obsSpaceData
344 call obs_finalize(obsSpaceData)
345
346 ! Job termination
347 istamp = exfin('VAR1D','FIN','NON')
348
349 call utl_tmg_stop(0)
350
351 call tmg_terminate(mmpi_myid, 'TMG_INFO')
352
353 call rpn_comm_finalize(ierr)
354
355end program midas_var1D
