1module oMinusF_mod
2 ! MODULE oMinusF_mod (prefix='omf' category='1. High-level functionality')
3 !
4 !:Purpose: Module for Observation minus Forecast (O-F) computation
5 !
6 use codePrecision_mod
7 use ramDisk_mod
8 use midasMpi_mod
9 use mathPhysConstants_mod
10 use horizontalCoord_mod
11 use verticalCoord_mod
12 use timeCoord_mod
13 use obsSpaceData_mod
14 use columnData_mod
15 use gridStateVector_mod
16 use gridStateVectorFileIO_mod
17 use obsFiles_mod
18 use innovation_mod
19 use obsErrors_mod
20 use biasCorrectionConv_mod
21 use obsSpaceErrorStdDev_mod
22 use ensembleObservations_mod
23 use ensembleStateVector_mod
24 use fileNames_mod
25 use statetocolumn_mod
26 implicit none
27 private
28
29 ! public subroutines and functions
30 public :: omf_oMinusF, omf_oMinusFens
31
32 contains
33
34 !--------------------------------------------------------------------------
35 ! omf_oMinusF
36 !--------------------------------------------------------------------------
37 subroutine omf_oMinusF(columnTrlOnAnlIncLev, columnTrlOnTrlLev, obsSpaceData, &
38 midasMode, addHBHT, addSigmaO)
39 !
40 !:Purpose: compute Observation-minus-Forecast (OmF)
41 !
42 implicit none
43
44 ! Arguments:
45 type(struct_columnData),target, intent(inout) :: columnTrlOnAnlIncLev
46 type(struct_columnData),target, intent(inout) :: columnTrlOnTrlLev
47 type(struct_obs), target, intent(inout) :: obsSpaceData
48 character(len=*), intent(in) :: midasMode
49 logical, intent(in) :: addHBHT
50 logical, intent(in) :: addSigmaO
51
52 ! Locals:
53 type(struct_gsv) :: stateVectorTrialHighRes
54 type(struct_vco), pointer :: vco_anl => null()
55 type(struct_hco), pointer :: hco_anl => null()
56 type(struct_hco), pointer :: hco_trl => null()
57 type(struct_vco), pointer :: vco_trl => null()
58 type(struct_hco), pointer :: hco_core => null()
59 logical :: allocHeightSfc
60 character(len=48) :: obsMpiStrategy
61 character(len=3) :: obsColumnMode
62 character(len=10) :: trialFileName
63 integer :: dateStampFromObs, get_max_rss
64
65 write(*,*) " ---------------------------------------"
66 write(*,*) " --- START OF SUBROUTINE oMinusF ---"
67 write(*,*) " --- Computation of the innovation ---"
68 write(*,*) " ---------------------------------------"
69
70 !
71 !- 1. Settings and module initializations
72 !
73 write(*,*)
74 write(*,*) '> omf_oMinusF: setup - START'
75
76 obsMpiStrategy = 'LIKESPLITFILES'
77 obsColumnMode = 'VAR'
78 trialFileName = './trlm_01'
79
80 !- 1.3 RAM disk usage
81 call ram_setup
82
83 !- 1.4 Temporal grid and dateStamp from trial file
84 call tim_setup(fileNameForDate_opt=trim(trialFileName))
85
86 !- 1.5 Observation file names and get datestamp, but do not use it
87 call obsf_setup(dateStampFromObs, trim(midasMode))
88
89 !- 1.6 Constants
90 if ( mmpi_myid == 0 ) then
91 call mpc_printConstants(6)
92 call pre_printPrecisions
93 end if
94
95 !- 1.7 Variables of the model states
96 call gsv_setup
97
98 !- 1.8 Set the horizontal domain
99 if ( addHBHT ) then
100 call hco_SetupFromFile(hco_anl, './analysisgrid', 'ANALYSIS') ! IN
101 if ( hco_anl % global ) then
102 hco_core => hco_anl
103 else
104 !- Iniatilized the core (Non-Exteded) analysis grid
105 call hco_SetupFromFile( hco_core, './analysisgrid', 'COREGRID') ! IN
106 end if
107 else
108 call hco_SetupFromFile(hco_anl, trim(trialFileName), ' ') ! IN
109 hco_core => hco_anl
110 end if
111
112 ! 1.9 Setup a column vector following the analysis vertical grid
113 if ( addHBHT ) then
114 call vco_SetupFromFile(vco_anl, & ! OUT
115 './analysisgrid') ! IN
116 call col_setVco(columnTrlOnAnlIncLev,vco_anl)
117 end if
118
119 !- 1.10 Setup and read observations
120 call inn_setupObs(obsSpaceData, hco_anl, obsColumnMode, obsMpiStrategy,trim(midasMode))!IN
121
122 ! Apply optional bias corrections when namelist logicals {fam}BiasActive are TRUE
123 ! (Only reverse existing corrections when namelist logicals {fam}RevOnly are TRUE)
124 if (obs_famExist(obsSpaceData,'AI')) call bcc_applyAIBcor(obsSpaceData)
125 if (obs_famExist(obsSpaceData,'GP')) call bcc_applyGPBcor(obsSpaceData)
126 if (obs_famExist(obsSpaceData,'UA')) call bcc_applyUABcor(obsSpaceData)
127
128 !- 1.11 Basic setup of columnData module
129 call col_setup
130
131 !- 1.12 Memory allocation for background column data
132 if ( addHBHT ) then
133 call col_allocate(columnTrlOnAnlIncLev, obs_numheader(obsSpaceData),mpiLocal_opt=.true.)
134 end if
135
136 if ( addSigmaO ) then
137 !- 1.13 Initialize the observation error covariances
138 write(*,*)
139 write(*,*) '> omf_oMinusF: Adding sigma_O'
140 call oer_setObsErrors(obsSpaceData, trim(midasMode))
141 end if
142
143 ! Reading trials
144 call inn_getHcoVcoFromTrlmFile( hco_trl, vco_trl )
145 allocHeightSfc = ( vco_trl%Vcode /= 0 )
146
147 call gsv_allocate( stateVectorTrialHighRes, tim_nstepobs, hco_trl, vco_trl, &
148 dateStamp_opt=tim_getDateStamp(), mpi_local_opt=.true., &
149 mpi_distribution_opt='Tiles', dataKind_opt=4, &
150 allocHeightSfc_opt=allocHeightSfc, hInterpolateDegree_opt='LINEAR', &
151 beSilent_opt=.false. )
152 call gsv_zero( stateVectorTrialHighRes )
153 call gio_readTrials( stateVectorTrialHighRes )
154 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
155
156 ! Horizontally interpolate trials to trial columns
157 call inn_setupColumnsOnTrlLev( columnTrlOnTrlLev, obsSpaceData, hco_core, &
158 stateVectorTrialHighRes )
159
160 write(*,*)
161 write(*,*) '> omf_oMinusF: setup - END'
162 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
163
164 !
165 !- 2. O-F computation
166 !
167
168 !- 2.1 Compute observation innovations
169 write(*,*)
170 write(*,*) '> omf_oMinusF: compute innovation'
171 call inn_computeInnovation(columnTrlOnTrlLev, obsSpaceData, analysisMode_opt=.false.)
172
173 if ( addHBHT ) then
174 write(*,*)
175 write(*,*) '> omf_oMinusF: Adding HBH^T'
176 !- 2.2 Interpolate background columns to analysis levels and setup for linearized H
177 call inn_setupColumnsOnAnlIncLev( columnTrlOnTrlLev, columnTrlOnAnlIncLev )
178 !- 2.3 Compute the background errors in observation space
179 call ose_computeStddev(columnTrlOnAnlIncLev,hco_anl,obsSpaceData)
180 end if
181
182 end subroutine omf_oMinusF
183
184 !--------------------------------------------------------------------------
185 ! omf_oMinusFens
186 !--------------------------------------------------------------------------
187 subroutine omf_oMinusFens(ensObs, obsSpaceData, nEns, ensPathName, &
188 midasMode, addHBHT, addSigmaO)
189 !
190 !:Purpose: compute Observation-minus-Forecast (OmF) for ensembles
191 !
192 implicit none
193
194 ! Arguments:
195 type(struct_eob), target, intent(inout) :: ensObs
196 type(struct_obs), target, intent(inout) :: obsSpaceData
197 integer, intent(in) :: nEns
198 character(len=*), intent(in) :: ensPathName
199 character(len=*), intent(in) :: midasMode
200 logical, intent(in) :: addHBHT
201 logical, intent(in) :: addSigmaO
202
203 ! Locals:
204 type(struct_columnData) :: columTrlOnTrlLev
205 type(struct_columnData) :: columnTrlOnAnlIncLev
206 type(struct_ens) :: ensembleTrl4D
207 type(struct_gsv) :: stateVector4D
208 type(struct_gsv) :: stateVectorWithZandP4D
209 type(struct_gsv) :: stateVectorHeightSfc
210 type(struct_vco), pointer :: vco_anl => null()
211 type(struct_hco), pointer :: hco_anl => null()
212 type(struct_hco), pointer :: hco_ens => null()
213 type(struct_vco), pointer :: vco_ens => null()
214 type(struct_hco), pointer :: hco_core => null()
215 character(len=256) :: ensFileName
216 character(len=48) :: obsMpiStrategy
217 character(len=3) :: obsColumnMode
218 integer, allocatable :: dateStampList(:)
219 integer :: datestamp, get_max_rss, memberIndex
220
221 write(*,*) " -----------------------------------------"
222 write(*,*) " --- START OF SUBROUTINE oMinusFens ---"
223 write(*,*) " --- Computation of the innovation ---"
224 write(*,*) " -----------------------------------------"
225
226 !
227 !- 1. Settings and module initializations
228 !
229 write(*,*)
230 write(*,*) '> omf_oMinusFens: setup - START'
231
232 obsMpiStrategy = 'LIKESPLITFILES'
233 obsColumnMode = 'VAR'
234
235 !- 1.3 RAM disk usage
236 call ram_setup
237
238 !- 1.4 Horizontal, vertical and temporal dimensions
239 call fln_ensFileName(ensFileName, ensPathName, memberIndex_opt=1)
240
241 call tim_setup( fileNameForDate_opt=trim(ensFileName) )
242 allocate(dateStampList(tim_nstepobs))
243 call tim_getstamplist(dateStampList,tim_nstepobs,tim_getDatestamp())
244
245 call hco_setupFromFile(hco_ens, ensFileName, ' ', 'ENSFILEGRID')
246 call vco_setupFromFile(vco_ens, ensFileName)
247
248 !- 1.5 Observation file names and get datestamp, but do not use it
249 call obsf_setup(dateStamp, trim(midasMode))
250
251 !- 1.6 Constants
252 if ( mmpi_myid == 0 ) then
253 call mpc_printConstants(6)
254 call pre_printPrecisions
255 end if
256
257 !- 1.7 Variables of the model states
258 call gsv_setup
259
260 !- 1.8 Set the horizontal domain
261 if ( addHBHT ) then
262 call hco_SetupFromFile(hco_anl, './analysisgrid', 'ANALYSIS') ! IN
263 if ( hco_anl % global ) then
264 hco_core => hco_anl
265 else
266 !- Iniatilized the core (Non-Exteded) analysis grid
267 call hco_SetupFromFile( hco_core, './analysisgrid', 'COREGRID') ! IN
268 end if
269 else
270 hco_core => hco_ens
271 end if
272
273 ! 1.9 Setup a column vector following the analysis vertical grid
274 if ( addHBHT ) then
275 call vco_SetupFromFile(vco_anl, & ! OUT
276 './analysisgrid') ! IN
277 call col_setVco(columnTrlOnAnlIncLev,vco_anl)
278 end if
279
280 !- 1.10 Setup and read observations
281 call inn_setupObs(obsSpaceData, hco_core, obsColumnMode, obsMpiStrategy, trim(midasMode)) !IN
282
283 ! Apply optional bias corrections when namelist logicals {fam}BiasActive are TRUE
284 ! (Only reverse existing corrections when namelist logicals {fam}RevOnly are TRUE)
285 if (obs_famExist(obsSpaceData,'AI')) call bcc_applyAIBcor(obsSpaceData)
286 if (obs_famExist(obsSpaceData,'GP')) call bcc_applyGPBcor(obsSpaceData)
287 if (obs_famExist(obsSpaceData,'UA')) call bcc_applyUABcor(obsSpaceData)
288
289 !- 1.11 Basic setup of columnData module
290 call col_setup
291 call col_setVco(columTrlOnTrlLev, vco_ens)
292 call col_allocate(columTrlOnTrlLev, obs_numheader(obsSpaceData))
293 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
294
295 !- 1.12 Memory allocation for background column data
296 if ( addHBHT ) then
297 call col_allocate(columnTrlOnAnlIncLev, obs_numheader(obsSpaceData),mpiLocal_opt=.true.)
298 end if
299
300 if ( addSigmaO ) then
301 !- 1.13 Initialize the observation error covariances
302 write(*,*)
303 write(*,*) '> omf_oMinusF: Adding sigma_O'
304 call oer_setObsErrors(obsSpaceData, trim(midasMode))
305 end if
306
307 !- 1.14 Allocate and initialize eob object for storing y-HX values
308 call eob_allocate(ensObs, nEns, obs_numBody(obsSpaceData), obsSpaceData)
309 call eob_zero(ensObs)
310 call eob_setLatLonObs(ensObs)
311
312 !- 1.15 Allocate statevector for storing state with heights and pressures allocated (for s2c_nl)
313 call gsv_allocate(stateVectorWithZandP4D, tim_nstepobs, hco_ens, vco_ens, &
314 dateStamp_opt=tim_getDateStamp(), &
315 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
316 dataKind_opt=4, allocHeightSfc_opt=.true.)
317 call gsv_zero(stateVectorWithZandP4D)
318
319 call gsv_allocate(stateVector4D, tim_nstepobs, hco_ens, vco_ens, &
320 dateStamp_opt=tim_getDateStamp(), &
321 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
322 dataKind_opt=4, allocHeightSfc_opt=.true., &
323 allocHeight_opt=.false., allocPressure_opt=.false.)
324 call gsv_zero(stateVector4D)
325
326 !- 1.16 Read the sfc height from ensemble member 1
327 call gsv_allocate(stateVectorHeightSfc, 1, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
328 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
329 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0','TT'/))
330 call gio_readFromFile(stateVectorHeightSfc, ensFileName, ' ', ' ', &
331 containsFullField_opt=.true., readHeightSfc_opt=.true.)
332
333 !- 1.17 Reading ensemble
334 call ens_allocate(ensembleTrl4D, nEns, tim_nstepobs, hco_ens, vco_ens, dateStampList)
335 call ens_readEnsemble(ensembleTrl4D, ensPathName, biPeriodic=.false.)
336
337 write(*,*)
338 write(*,*) '> omf_oMinusFens: setup - END'
339 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
340
341 !
342 !- 2. O-F computation
343 !
344 do memberIndex = 1, nEns
345 write(*,*) ''
346 write(*,*) 'oMinusFens: compute O-P for ensemble member ', memberIndex
347 write(*,*) 'Memory Used: ', get_max_rss()/1024, 'Mb'
348
349 !- 2.1 Copy selected member to a stateVector
350 call ens_copyMember(ensembleTrl4D, stateVector4D, memberIndex)
351 call gsv_copy(stateVector4D, stateVectorWithZandP4D, allowVarMismatch_opt=.true., &
352 beSilent_opt=.true.)
353 call gsv_copyHeightSfc(stateVectorHeightSfc, stateVectorWithZandP4D)
354
355 !- 2.2 Compute and set Yb in ensObs
356 call s2c_nl(stateVectorWithZandP4D, obsSpaceData, columTrlOnTrlLev, hco_ens, &
357 timeInterpType='LINEAR', dealloc_opt=.false., &
358 beSilent_opt=.true.)
359
360 !- 2.3 Compute innovation
361 call inn_computeInnovation(columTrlOnTrlLev, obsSpaceData, analysisMode_opt=.false., beSilent_opt=.true.)
362
363 !- 2.4 Copy to ensObs
364 call eob_setYb(ensObs, memberIndex)
365
366 !- 2.5 Add background errors in observation space
367 if (addHBHT .and. memberIndex == 1) then
368 write(*,*)
369 write(*,*) '> omf_oMinusFens: Adding HBH^T'
370 ! Interpolate background columns to analysis levels and setup for linearized H
371 call inn_setupColumnsOnAnlIncLev(columTrlOnTrlLev, columnTrlOnAnlIncLev)
372 ! Compute the background errors in observation space
373 call ose_computeStddev(columnTrlOnAnlIncLev,hco_anl,obsSpaceData)
374 end if
375
376 end do
377
378 !
379 !- 3. Ending
380 !
381 call gsv_deallocate(stateVectorWithZandP4D)
382 call gsv_deallocate(stateVector4D)
383 call gsv_deallocate(stateVectorHeightSfc)
384
385 end subroutine omf_oMinusFens
386
387end module oMinusF_mod
