1module analysisErrorOI_mod
2 ! MODULE analysisErrorOI_mod (prefix='aer' category='1. High-level functionality')
3 !
4 !:Purpose: Calculate the analysis-error standard deviation.
5 ! The method used is Optimal Interpolation,
6 ! where it is assumed that only a subset of the
7 ! total number of observations influence the analysis at a given grid point.
8 !
9 use columnData_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use kdTree2_mod
13 use midasMpi_mod
14 use obsSpaceData_mod
15 use physicsFunctions_mod
16 use stateToColumn_mod
17 use varNamelist_mod
18 use utilities_mod
19 use horizontalCoord_mod
20 use verticalCoord_mod
21 use obsOperators_mod
22 use message_mod
23 use oceanMask_mod
24 use timeCoord_mod
25
26 implicit none
27
28 private
29
30 ! public subroutines and functions
31 public :: aer_analysisError
32
33 type struct_neighborhood
34 integer :: numObs
35 integer, pointer :: headerIndex(:)
36 integer, pointer :: bodyIndex(:)
37 integer, pointer :: procIndex(:)
38 end type struct_neighborhood
39
40 integer, external :: get_max_rss
41 integer, parameter :: maxNumLocalGridptsSearch = 1500
42
43contains
44
45 !--------------------------------------------------------------------------
46 ! aer_analysisError
47 !--------------------------------------------------------------------------
48 subroutine aer_analysisError(obsSpaceData, hco_ptr, vco_ptr)
49 !
50 !:Purpose: Calculate analysis-error variance.
51 !
52 implicit none
53
54 ! Arguments:
55 type(struct_obs), intent(inout) :: obsSpaceData ! observation data structure
56 type(struct_hco), pointer, intent(in) :: hco_ptr ! horizontal grid definition
57 type(struct_vco), pointer, intent(in) :: vco_ptr ! vertical grid definition
58
59 ! Locals:
60 integer :: fnom, fclos, nulnam, ierr
61 type(struct_gsv) :: stateVectorAnlErrorStd ! state vector for analysis error std deviation
62 type(struct_gsv) :: stateVectorTrlErrorStd ! state vector for background error std deviation
63 real(8), pointer :: anlErrorStdDev_ptr(:,:,:,:) ! pointer for analysis error std deviation
64 real(8), pointer :: trlErrorStdDev_ptr(:,:,:,:) ! pointer for background error std deviation
65 integer, allocatable :: numObs(:,:)
66 integer :: lonIndex, latIndex
67 character(len=4), pointer :: analysisVariable(:)
68 type(struct_gsv) :: statevectorLcorr
69 real(4), pointer :: field3D_r4_ptr(:,:,:)
70 type(struct_columnData) :: column
71 type(struct_columnData) :: columng
72 type(struct_ocm) :: oceanMask
73 real(8) :: leadTimeInHours
74 character(len=20), parameter :: errStddevFileName_in = 'errorstdev_in' ! input filename for anl ot trl error standard deviation
75 character(len=20), parameter :: anlErrStddevFileName_out = 'anlerrorstdev_out' ! output filename for anl error std deviation
76 character(len=20), parameter :: trlErrStddevFileName_out = 'trlerrorstdev_out' ! output filename for trl (background) error std deviation
77 type(struct_neighborhood), pointer :: influentObs(:,:)
78 real(8), allocatable :: Lcorr(:,:)
79
80 ! namelist variables:
81 real(8) :: maxAnalysisErrorStdDev ! maximum limit imposed on analysis error stddev
82 logical :: propagateAnalysisError ! choose to propagate analysis error
83 logical :: propagateDSLO ! choose to propagate Days Since Last Obs field
84 real(4) :: errorGrowth ! seaice: fraction of ice per hour, SST: estimated growth
85 character(len=12) :: analysisEtiket ! analysis field etiket in a standard file
86 character(len=12) :: anlErrorStdEtiket ! analysis error standard deviation field etiket in the input/output standard files
87 character(len=12) :: trlErrorStdEtiket ! background error standard deviation field etiket in the input/output standard files
88 integer :: hoursSinceLastAnalysis ! number of hours since the last analysis
89 logical :: saveTrlStdField ! choose to save trial standard deviation field
90 character(len=2) :: inputTypeVar ! typvar of the analysis error field in the input file
91 character(len=2) :: outputTypeVar ! typvar of the analysis error field for the output file
92 real(4) :: multFactorLcorr ! multiplication scaling factor to increase the correlation length scale field
93 namelist /namaer/ maxAnalysisErrorStdDev, propagateAnalysisError, propagateDSLO, &
94 errorGrowth, analysisEtiket, anlErrorStdEtiket, trlErrorStdEtiket, &
95 hoursSinceLastAnalysis, saveTrlStdField, inputTypeVar, outputTypeVar, &
96 multFactorLcorr
97
98 write(*,*) '**********************************************************'
99 write(*,*) '** aer_analysisError: Calculate analysis-error variance **'
100 write(*,*) '**********************************************************'
101
102 ! default namelist variable values
103 maxAnalysisErrorStdDev = 1.0d0
104 propagateAnalysisError = .false.
105 propagateDSLO = .false.
106 errorGrowth = 1.0
107 analysisEtiket = ''
108 anlErrorStdEtiket = 'A-ER STD DEV'
109 trlErrorStdEtiket = 'B-ER STD DEV'
110 hoursSinceLastAnalysis = 6
111 saveTrlStdField = .false.
112 inputTypeVar = 'P@'
113 outputTypeVar = 'A@'
114 multFactorLcorr = 1.0
115
116 ! read the namelist
117 if (.not. utl_isNamelistPresent('namaer','./flnml')) then
118 if (mmpi_myid == 0) then
119 call msg('aer_analysisError:', ' namaer is missing in the namelist.')
120 call msg('aer_analysisError:', ' the default values will be taken.')
121 end if
122 else
123 ! reading namelist variables
124 nulnam = 0
125 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
126 read(nulnam, nml = namaer, iostat = ierr)
127 if (ierr /= 0) call utl_abort('aer_analysisError:: Error reading namelist')
128 ierr = fclos(nulnam)
129 end if
130 write(*, nml = namaer)
131
132 nullify(analysisVariable)
133 call gsv_varNamesList(analysisVariable)
134
135 if (size(analysisVariable) > 1) then
136 call utl_abort('aer_analysisError: Check namelist NAMSTATE. analysisVariable is greater than 1.')
137 end if
138
139 if (analysisVariable(1) == 'GL') then
140 call msg('aer_analysisError:', ' computing seaice analysis error...')
141 else if(analysisVariable(1) == 'TM') then
142 call msg('aer_analysisError:', ' computing SST analysis error...')
143 else
144 call utl_abort('aer_analysisError:: The current code does not work with '&
145 //trim(analysisVariable(1))//' analysis variable.')
146 end if
147
148 call gsv_allocate(stateVectorAnlErrorStd, 1, hco_ptr, vco_ptr, dateStamp_opt = -1, &
149 mpi_local_opt = .false., mpi_distribution_opt = 'None', &
150 varNames_opt = (/analysisVariable(1)/), dataKind_opt = 8)
151 call gsv_getField(stateVectorAnlErrorStd, anlErrorStdDev_ptr)
152 call gsv_allocate(stateVectorTrlErrorStd, 1, hco_ptr, vco_ptr, dateStamp_opt = -1, &
153 mpi_local_opt = .false., mpi_distribution_opt = 'None', &
154 varNames_opt = (/analysisVariable(1)/), dataKind_opt = 8)
155 call gsv_getField(stateVectorTrlErrorStd, trlErrorStdDev_ptr)
156
157 if (propagateAnalysisError) then
158 call msg('aer_analysisError:', &
159 ' analysis error std field is read from: '//trim(errStddevFileName_in))
160 call gio_readFromFile(stateVectorAnlErrorStd, errStddevFileName_in, &
161 etiket_in = anlErrorStdEtiket, typvar_in = inputTypeVar, &
162 containsFullField_opt = .false.)
163
164 ! initialize trl error std deviation field:
165 trlErrorStdDev_ptr(:,:,:,:) = anlErrorStdDev_ptr(:,:,:,:)
166
167 call ocm_readMaskFromFile (oceanMask, hco_ptr, vco_ptr, errStddevFileName_in)
168 call aer_propagateAnalysisError (stateVectorTrlErrorStd, oceanMask, &
169 analysisVariable(1), &
170 analysisEtiket, errorGrowth, &
171 hco_ptr, vco_ptr)
172
173 ! impose maximum value on trial error standard deviation field
174 trlErrorStdDev_ptr(:,:,:,:) = min(trlErrorStdDev_ptr(:,:,:,:), &
175 maxAnalysisErrorStdDev)
176
177 if (saveTrlStdField) then
178 ! zap analysis error etiket with background error etiket
179 stateVectorTrlErrorStd%etiket = trlErrorStdEtiket
180
181 ! copy mask from analysis error std deviation field to trl error std field
182 call gsv_copyMask(stateVectorAnlErrorStd, stateVectorTrlErrorStd)
183
184 ! update dateStamp from env variable
185 call gsv_modifyDate(stateVectorTrlErrorStd, tim_getDateStamp(), &
186 modifyDateOrigin_opt = .true.)
187
188 ! save background error (increased analysis error) standard deviation field
189 call gio_writeToFile(stateVectorTrlErrorStd, trlErrStddevFileName_out, &
190 trlErrorStdEtiket, typvar_opt = outputTypeVar, &
191 containsFullField_opt = .false.)
192 end if
193
194 else
195 call msg('aer_analysisError:', &
196 ' trial error std field is read from: '//trim(errStddevFileName_in))
197 call gio_readFromFile(stateVectorTrlErrorStd, errStddevFileName_in, &
198 etiket_in = trlErrorStdEtiket, typvar_in = inputTypeVar, &
199 containsFullField_opt = .false.)
200 call gsv_copyMask(stateVectorTrlErrorStd, stateVectorAnlErrorStd)
201 end if
202
203 leadTimeInHours = real(stateVectorTrlErrorStd%deet * stateVectorTrlErrorStd%npasList(1), 8) / 3600.0d0
204 call incdatr(stateVectorAnlErrorStd%dateOriginList(1), &
205 stateVectorTrlErrorStd%dateOriginList(1), leadTimeInHours)
206
207 stateVectorAnlErrorStd%etiket = anlErrorStdEtiket
208
209 call col_setVco(column, vco_ptr)
210 call col_allocate(column, obs_numHeader(obsSpaceData), varNames_opt=(/analysisVariable(1)/))
211 call col_setVco(columng, vco_ptr)
212 call col_allocate(columng, obs_numHeader(obsSpaceData), varNames_opt=(/analysisVariable(1)/))
213 call s2c_tl(stateVectorTrlErrorStd, column, columng, obsSpaceData)
214
215 allocate(Lcorr(hco_ptr%ni, hco_ptr%nj))
216
217 ! get correlation length scale field
218 write(*,*) 'aer_analysisError: get correlation length scale field...'
219 call gsv_allocate(statevectorLcorr, 1, hco_ptr, vco_ptr, dateStamp_opt = -1, &
220 mpi_local_opt = .false., mpi_distribution_opt = 'None', &
221 dataKind_opt = 4, hInterpolateDegree_opt = 'LINEAR', &
222 varNames_opt = (/analysisVariable(1)/))
223 call gsv_zero(statevectorLcorr)
224 call gio_readFromFile(statevectorLcorr, './bgstddev', 'CORRLEN', ' ', &
225 unitConversion_opt = .false.)
226 call gsv_getField(statevectorLcorr, field3D_r4_ptr, analysisVariable(1))
227
228 ! apply multiplication scaling factor
229 field3D_r4_ptr(:, :, 1) = field3D_r4_ptr(:, :, 1) * multFactorLcorr
230
231 ! Convert from km to meters
232 Lcorr(:,:) = 1000.0d0 * real(field3D_r4_ptr(:, :, 1), 8)
233
234 write(*,*) 'aer_analysisError: min/max correlation length scale 2D field: ', &
235 minval(Lcorr(:,:) ), maxval(Lcorr(:,:))
236 call gsv_deallocate(statevectorLcorr)
237
238 write(*,*) 'aer_analysisError: done creating kdtree for stateVectorTrlErrorStd'
239 write(*,*) 'Memory Used: ', get_max_rss() / 1024, 'Mb'
240
241 ! Go through all observations a first time to get
242 ! the number of influential observations
243 ! in order to allocate memory appropriately.
244
245 allocate(influentObs(hco_ptr%ni, hco_ptr%nj))
246 allocate(numObs(hco_ptr%ni, hco_ptr%nj))
247
248 call msg('aer_analysisError:', ' looking for observations...')
249 call findObs(obsSpaceData, stateVectorTrlErrorStd, numObs, &
250 analysisVariable(1), Lcorr, influentObs)
251 write(*,*) 'aer_analysisError: min/max of numObs = ', minval(numObs), maxval(numObs)
252
253 ! Memory allocation
254
255 do latIndex = 1, hco_ptr%nj
256 do lonIndex = 1, hco_ptr%ni
257 influentObs(lonIndex, latIndex)%numObs = numObs(lonIndex,latIndex)
258 allocate(influentObs(lonIndex, latIndex)%headerIndex(numObs(lonIndex, latIndex)))
259 allocate(influentObs(lonIndex, latIndex)%bodyIndex(numObs(lonIndex, latIndex)))
260 allocate(influentObs(lonIndex, latIndex)%procIndex(numObs(lonIndex, latIndex)))
261 end do
262 end do
263
264 ! Go through all observations a second time to get
265 ! the indexes of the influential observations.
266 ! This is not efficient to go through all observations twice
267 ! but it saves lot of memory space.
268
269 call msg('aer_analysisError:', ' go through all observations a second time...')
270 call findObs(obsSpaceData, stateVectorTrlErrorStd, numObs, &
271 analysisVariable(1), Lcorr, influentObs)
272
273 deallocate(numObs)
274
275 write(*,*) 'aer_analysisError:: obs found'
276 write(*,*) 'Memory Used: ',get_max_rss() / 1024, 'Mb'
277
278 ! compute analysis error standard deviation
279 call aer_computeAnlErrorStd(obsSpaceData, stateVectorAnlErrorStd, stateVectorTrlErrorStd, &
280 analysisVariable(1), maxAnalysisErrorStdDev, influentObs, Lcorr)
281
282 deallocate(influentObs)
283 deallocate(Lcorr)
284
285 ! update dateStamp from env variable
286 call gsv_modifyDate(stateVectorAnlErrorStd, tim_getDateStamp(), &
287 modifyDateOrigin_opt = .true.)
288
289 ! save analysis error
290 call msg('aer_analysisError:', ' writing analysis error std field to output file...')
291 if (mmpi_myid == 0) then
292 call gio_writeToFile(stateVectorAnlErrorStd, anlErrStddevFileName_out, &
293 anlErrorStdEtiket, typvar_opt = outputTypeVar, &
294 containsFullField_opt = .false.)
295 end if
296
297 call col_deallocate(columng)
298 call col_deallocate(column)
299 call gsv_deallocate(stateVectorAnlErrorStd)
300 call gsv_deallocate(stateVectorTrlErrorStd)
301
302 if (analysisVariable(1) == 'GL') then
303 ! Update the Days Since Last Obs
304 call aer_daysSinceLastObs(obsSpaceData, hco_ptr, vco_ptr, &
305 errStddevFileName_in, anlErrStddevFileName_out, &
306 analysisVariable(1), propagateDSLO, &
307 hoursSinceLastAnalysis)
308 end if
309
310 call msg('aer_analysisError:', ' finished.')
311
312 end subroutine aer_analysisError
313
314 !---------------------------------------------------------
315 ! findObs
316 !---------------------------------------------------------
317 subroutine findObs(obsSpaceData, stateVectorTrlErrorStd, numObs, variableName, Lcorr, influentObs)
318 !
319 !:Purpose: Find all observations used for the analysis.
320 !
321 implicit none
322
323 ! Arguments:
324 type(struct_obs), intent(in) :: obsSpaceData ! observation data structure
325 type(struct_gsv), intent(in) :: stateVectorTrlErrorStd ! state containing background error stddev
326 integer, intent(out) :: numObs(:,:) ! number of observations found
327 character(len=*), intent(in) :: variableName ! 'GL' for seaice or 'TM' for SST
328 real(8) , intent(in) :: Lcorr(:,:) ! horizontal background-error correlation length scale
329 type(struct_neighborhood), pointer, intent(inout) :: influentObs(:,:) ! details about observations to use in update
330
331 ! Locals:
332 integer :: headerIndex, bodyIndexBeg, bodyIndexEnd, bodyIndex
333 integer :: procIndex, kIndex, stepIndex
334 integer :: gridptCount, gridpt, numLocalGridptsFoundSearch
335 integer :: lonIndex, latIndex, resultsIndex, gridIndex, numStep
336 type(kdtree2_result) :: searchResults(maxNumLocalGridptsSearch)
337 real(kdkind) :: refPosition(3), maxRadiusSquared
338 real(4) :: footprintRadius_r4 ! (metres) used for seaice observations only
339 real(4) :: influenceRadius_r4 ! (metres)
340 real(8) :: obsLonInRad, obsLatInRad, maxLcorr
341 type(kdtree2), pointer :: tree
342 real(kdkind), allocatable :: positionArray(:,:)
343 real(8), allocatable :: latInRad(:,:), lonInRad(:,:)
344 real(8) :: interpWeight(maxNumLocalGridptsSearch)
345 integer :: obsLatIndex(maxNumLocalGridptsSearch), obsLonIndex(maxNumLocalGridptsSearch)
346 integer, allocatable :: obsAss(:), allObsAss(:,:)
347 integer, allocatable :: obsRln(:), allObsRln(:,:)
348 integer, allocatable :: obsNlv(:), allObsNlv(:,:)
349 integer :: ierr, numHeaderMax, allNumHeader(mmpi_nprocs), numBodyMax, allNumBody(mmpi_nprocs)
350 real(4), allocatable :: footprintRadiusVec_r4(:), allFootprintRadius_r4(:,:)
351 real(8), allocatable :: obsLon(:), allObsLon(:,:), obsLat(:), allObsLat(:,:)
352
353 numStep = stateVectorTrlErrorStd%numStep
354 if (numStep > 1) call utl_abort('aer findObs: Code is not adapted for numStep > 1')
355
356 ! Communicate some quantities to all MPI tasks
357
358 call rpn_comm_allgather(obs_numHeader(obsSpaceData), 1, 'mpi_integer', &
359 allNumHeader, 1, 'mpi_integer', 'GRID', ierr)
360 numHeaderMax = maxval(allNumHeader)
361 call rpn_comm_allgather(obs_numBody(obsSpaceData), 1, 'mpi_integer', &
362 allNumBody, 1, 'mpi_integer', 'GRID', ierr)
363 numBodyMax = maxval(allNumBody)
364
365 allocate(obsAss(numBodyMax))
366 obsAss(:) = 0
367 allocate(allObsAss(numBodyMax, mmpi_nprocs))
368 do bodyIndex = 1, obs_numBody(obsSpaceData)
369 obsAss(bodyIndex) = obs_notassimilated
370 if (obs_bodyElem_i(obsSpaceData, obs_ass, bodyIndex) == obs_assimilated) then
371 obsAss(bodyIndex) = obs_assimilated
372 end if
373 end do
374 call rpn_comm_gather(obsAss, numBodyMax, 'mpi_integer', &
375 allObsAss, numBodyMax, 'mpi_integer', 0, 'grid', ierr)
376
377 allocate(obsRln(numHeaderMax))
378 obsRln(:) = 0
379 allocate(obsNlv(numHeaderMax))
380 obsNlv(:) = 0
381 allocate(obsLon(numHeaderMax))
382 obsLon(:) = 0.0d0
383 allocate(obsLat(numHeaderMax))
384 obsLat(:) = 0.0d0
385 do headerIndex = 1, obs_numHeader(obsSpaceData)
386 obsRln(headerIndex) = obs_headElem_i(obsSpaceData, obs_rln, headerIndex)
387 obsNlv(headerIndex) = obs_headElem_i(obsSpaceData, obs_nlv, headerIndex)
388 obsLon(headerIndex) = obs_headElem_r(obsSpaceData, obs_lon, headerIndex)
389 obsLat(headerIndex) = obs_headElem_r(obsSpaceData, obs_lat, headerIndex)
390 end do
391 allocate(allObsRln(numHeaderMax,mmpi_nprocs))
392 allocate(allObsNlv(numHeaderMax,mmpi_nprocs))
393 allocate(allObsLon(numHeaderMax,mmpi_nprocs))
394 allocate(allObsLat(numHeaderMax,mmpi_nprocs))
395 call rpn_comm_gather(obsRln, numHeaderMax, 'mpi_integer', &
396 allObsRln, numHeaderMax, 'mpi_integer', 0, 'grid', ierr)
397 call rpn_comm_gather(obsNlv, numHeaderMax, 'mpi_integer', &
398 allObsNlv, numHeaderMax, 'mpi_integer', 0, 'grid', ierr)
399 call rpn_comm_gather(obsLon, numHeaderMax, 'mpi_real8', &
400 allObsLon, numHeaderMax, 'mpi_real8', 0, 'grid', ierr)
401 call rpn_comm_gather(obsLat, numHeaderMax, 'mpi_real8', &
402 allObsLat, numHeaderMax, 'mpi_real8', 0, 'grid', ierr)
403
404 allocate(footprintRadiusVec_r4(numHeaderMax))
405 do headerIndex = 1, obs_numHeader(obsSpaceData)
406 footprintRadiusVec_r4(headerIndex) = s2c_getFootprintRadius(obsSpaceData, stateVectorTrlErrorStd, headerIndex)
407 end do
408 allocate(allFootprintRadius_r4(numHeaderMax,mmpi_nprocs))
409 call rpn_comm_gather(footprintRadiusVec_r4, numHeaderMax, 'MPI_REAL4', &
410 allFootprintRadius_r4(:,:), numHeaderMax, 'MPI_REAL4', &
411 0, 'GRID', ierr)
412
413 ! create kdtree
414 write(*,*) 'findObs: start creating kdtree for stateVectorTrlErrorStd'
415 write(*,*) 'Memory Used: ', get_max_rss() / 1024, 'Mb'
416
417 allocate(positionArray(3, stateVectorTrlErrorStd%hco%ni * stateVectorTrlErrorStd%hco%nj))
418 allocate(lonInRad(stateVectorTrlErrorStd%hco%ni, stateVectorTrlErrorStd%hco%nj))
419 allocate(latInRad(stateVectorTrlErrorStd%hco%ni, stateVectorTrlErrorStd%hco%nj))
420
421 gridIndex = 0
422 do latIndex = 1, stateVectorTrlErrorStd%hco%nj
423 do lonIndex = 1, stateVectorTrlErrorStd%hco%ni
424 gridIndex = gridIndex + 1
425 latInRad(lonIndex, latIndex) = real(stateVectorTrlErrorStd%hco%lat2d_4(lonIndex, latIndex), 8)
426 lonInRad(lonIndex, latIndex) = real(stateVectorTrlErrorStd%hco%lon2d_4(lonIndex, latIndex), 8)
427 positionArray(:, gridIndex) = kdtree2_3dPosition(lonInRad(lonIndex, latIndex), &
428 latInRad(lonIndex, latIndex))
429 end do
430 end do
431
432 write(*,*) 'findObs: latInRad min/max: ', minval(latInRad(:,:)), maxval(latInRad(:,:))
433 write(*,*) 'findObs: lonInRad min/max: ', minval(lonInRad(:,:)), maxval(lonInRad(:,:))
434
435 nullify(tree)
436 tree => kdtree2_create(positionArray, sort = .false., rearrange = .true.)
437
438 deallocate(positionArray)
439 deallocate(lonInRad)
440 deallocate(latInRad)
441
442 maxLcorr = maxval(Lcorr(:,:))
443 numObs(:,:) = 0
444
445 if (mmpi_myid == 0) then
446 PROC_LOOP: do procIndex = 1, mmpi_nprocs
447 HEADER_LOOP: do headerIndex = 1, allNumHeader(procIndex)
448
449 bodyIndexBeg = allObsRln(headerIndex,procIndex)
450 bodyIndexEnd = allObsNlv(headerIndex,procIndex) + bodyIndexBeg - 1
451
452 BODY_LOOP: do bodyIndex = bodyIndexBeg, bodyIndexEnd
453
454 if (allObsAss(bodyIndex,procIndex) /= obs_assimilated) then
455 cycle BODY_LOOP
456 end if
457
458 if (trim(variableName) == 'GL') then
459 footprintRadius_r4 = allFootPrintRadius_r4(headerIndex, procIndex)
460 influenceRadius_r4 = max(0.0, footprintRadius_r4) + maxLcorr
461 else if (trim(variableName) == 'TM') then
462 influenceRadius_r4 = maxLcorr
463 else
464 call utl_abort('findObs: The current code does not work with '&
465 //trim(variableName)//' analysis variable.')
466 end if
467
468 if (maxLcorr == 0.0d0) then
469
470 do kIndex = stateVectorTrlErrorStd%mykBeg, stateVectorTrlErrorStd%mykEnd
471 do stepIndex = 1, stateVectorTrlErrorStd%numStep
472
473 call s2c_getWeightsAndGridPointIndexes(headerIndex, kIndex, stepIndex, &
474 procIndex, interpWeight, obsLatIndex, &
475 obsLonIndex, gridptCount)
476
477 GRIDPT_LOOP: do gridpt = 1, gridptCount
478
479 if (interpWeight(gridpt) == 0.0d0) cycle GRIDPT_LOOP
480
481 lonIndex = obsLonIndex(gridpt)
482 latIndex = obsLatIndex(gridpt)
483
484 numObs(lonIndex, latIndex) = numObs(lonIndex, latIndex) + 1
485 if(associated(influentObs(lonIndex, latIndex)%bodyIndex)) then
486 if(numObs(lonIndex, latIndex) > influentObs(lonIndex, latIndex)%numObs) then
487 call utl_abort('findObs: Array too small')
488 end if
489 influentObs(lonIndex, latIndex)%headerIndex(numObs(lonIndex, latIndex)) = headerIndex
490 influentObs(lonIndex, latIndex)%bodyIndex(numObs(lonIndex, latIndex)) = bodyIndex
491 influentObs(lonIndex, latIndex)%procIndex(numObs(lonIndex, latIndex)) = procIndex
492 end if
493
494 end do GRIDPT_LOOP
495
496 end do
497 end do
498
499 else
500
501 ! Determine the grid point nearest the observation.
502
503 obsLonInRad = allObsLon(headerIndex, procIndex)
504 obsLatInRad = allObsLat(headerIndex, procIndex)
505
506 ! do the search
507 maxRadiusSquared = real(influenceRadius_r4, 8) ** 2
508 refPosition(:) = kdtree2_3dPosition(obsLonInRad, obsLatInRad)
509 call kdtree2_r_nearest(tp = tree, qv = refPosition, r2 = maxRadiusSquared, &
510 nfound = numLocalGridptsFoundSearch, &
511 nalloc = maxNumLocalGridptsSearch, &
512 results = searchResults)
513
514 if (numLocalGridptsFoundSearch > maxNumLocalGridptsSearch) then
515 call utl_abort('findObs: parameter maxNumLocalGridptsSearch must be increased')
516 end if
517
518 do resultsIndex = 1, numLocalGridptsFoundSearch
519
520 gridIndex = searchResults(resultsIndex)%idx
521 if (gridIndex < 1 .or. &
522 gridIndex > stateVectorTrlErrorStd%hco%ni * stateVectorTrlErrorStd%hco%nj) then
523 write(*,*) 'analysisErrorStdDev: gridIndex = ', gridIndex
524 call utl_abort('findObs: gridIndex out of bound.')
525 end if
526
527 latIndex = (gridIndex - 1) / stateVectorTrlErrorStd%hco%ni + 1
528 lonIndex = gridIndex - (latIndex - 1) * stateVectorTrlErrorStd%hco%ni
529 if (lonIndex < 1 .or. lonIndex > stateVectorTrlErrorStd%hco%ni .or. &
530 latIndex < 1 .or. latIndex > stateVectorTrlErrorStd%hco%nj) then
531 write(*,*) 'analysisErrorStdDev: lonIndex = ', lonIndex, &
532 ', latIndex = ', latIndex
533 call utl_abort('findObs: lonIndex/latIndex out of bound.')
534 end if
535
536 numObs(lonIndex, latIndex) = numObs(lonIndex, latIndex) + 1
537 if (associated(influentObs(lonIndex, latIndex)%bodyIndex)) then
538 if (numObs(lonIndex, latIndex) > influentObs(lonIndex, latIndex)%numObs) then
539 call utl_abort('findObs: Array too small in subroutine findObs')
540 end if
541 influentObs(lonIndex, latIndex)%headerIndex(numObs(lonIndex, latIndex)) = headerIndex
542 influentObs(lonIndex, latIndex)%bodyIndex(numObs(lonIndex, latIndex)) = bodyIndex
543 influentObs(lonIndex, latIndex)%procIndex(numObs(lonIndex, latIndex)) = procIndex
544 end if
545
546 end do
547
548 end if
549
550 end do BODY_LOOP
551
552 end do HEADER_LOOP
553
554 end do PROC_LOOP
555 end if
556
557 ! Communicate values from proc 0 to others
558 if (mmpi_nprocs > 1) then
559 call rpn_comm_bcast(numObs, size(numObs), 'MPI_INTEGER', 0, 'GRID', ierr)
560 if (associated(influentObs(1,1)%bodyIndex)) then
561 write(*,*) 'findObs: communicate bodyIndex, headerIndex'
562 do latIndex = 1, stateVectorTrlErrorStd%hco%nj
563 do lonIndex = 1, stateVectorTrlErrorStd%hco%ni
564 call rpn_comm_bcast(influentObs(lonIndex, latIndex)%headerIndex, &
565 influentObs(lonIndex, latIndex)%numObs, 'MPI_INTEGER', 0, 'GRID', ierr)
566 call rpn_comm_bcast(influentObs(lonIndex, latIndex)%bodyIndex, &
567 influentObs(lonIndex, latIndex)%numObs, 'MPI_INTEGER', 0, 'GRID', ierr)
568 call rpn_comm_bcast(influentObs(lonIndex, latIndex)%procIndex, &
569 influentObs(lonIndex, latIndex)%numObs, 'MPI_INTEGER', 0, 'GRID', ierr)
570 end do
571 end do
572 end if
573 end if
574
575 deallocate(obsAss)
576 deallocate(allObsAss)
577 deallocate(obsRln)
578 deallocate(obsNlv)
579 deallocate(obsLon)
580 deallocate(obsLat)
581 deallocate(allObsRln)
582 deallocate(allObsNlv)
583 deallocate(allObsLon)
584 deallocate(allObsLat)
585 deallocate(footprintRadiusVec_r4)
586 deallocate(allFootprintRadius_r4)
587
588 end subroutine findObs
589
590 !---------------------------------------------------------
591 ! aer_daysSinceLastObs
592 !---------------------------------------------------------
593 subroutine aer_daysSinceLastObs(obsSpaceData, hco_ptr, vco_ptr, &
594 inputFileName, outputFileName, &
595 variableName, propagateDSLO, &
596 hoursSinceLastAnalysis)
597 !
598 !:Purpose: Update the field "days since last obs" with the newly assimilated obs.
599 !
600 implicit none
601
602 ! Arguments:
603 type(struct_obs), intent(inout) :: obsSpaceData ! observation data structure
604 type(struct_hco), pointer, intent(in) :: hco_ptr ! horizontal grid definition
605 type(struct_vco), pointer, intent(in) :: vco_ptr ! vertical grid definition
606 character(len=*), intent(in) :: inputFileName ! input file name
607 character(len=*), intent(in) :: outputFileName ! output file name
608 character(len=*), intent(in) :: variableName ! name of variable being treated
609 logical , intent(in) :: propagateDSLO ! propagate (increase) Days Since Last Obs in time
610 integer , intent(in) :: hoursSinceLastAnalysis ! number of hours between analysis times
611
612 ! Locals:
613 type(struct_gsv) :: stateVectorTrlDSLO, stateVectorAnlDSLO
614 real(8), pointer :: trlDSLO_ptr(:,:,:,:), anlDSLO_ptr(:,:,:,:)
615 integer :: stepIndex, levIndex, lonIndex, latIndex, headerIndex
616 integer :: bodyIndexBeg, bodyIndexEnd, bodyIndex, kIndex, procIndex
617 integer :: gridptCount, gridpt
618 type(struct_columnData) :: column, columng
619 real(8) :: leadTimeInHours, interpWeight(maxNumLocalGridptsSearch)
620 integer :: obsLatIndex(maxNumLocalGridptsSearch), obsLonIndex(maxNumLocalGridptsSearch)
621 integer :: ierr, numHeaderMax, allNumHeader(mmpi_nprocs)
622 integer, allocatable :: obsAss(:), obsAssMpiGlobal(:,:)
623
624 write(*,*) '**********************************************************'
625 write(*,*) '** aer_daysSinceLastObs: Update the days since last obs **'
626 write(*,*) '**********************************************************'
627
628 write(*,*) 'aer_daysSinceLastObs: input variable: ', trim(variableName)
629
630 call gsv_allocate(stateVectorTrlDSLO, 1, hco_ptr, vco_ptr, dateStamp_opt = -1, &
631 mpi_local_opt = .false., mpi_distribution_opt = 'None', &
632 varNames_opt=(/'DSLO'/), dataKind_opt = 8)
633 call gsv_allocate(stateVectorAnlDSLO, 1, hco_ptr, vco_ptr, dateStamp_opt = -1, &
634 mpi_local_opt = .false., mpi_distribution_opt = 'None', &
635 varNames_opt = (/'DSLO'/), dataKind_opt = 8)
636
637 call rpn_comm_allgather(obs_numHeader(obsSpaceData), 1, 'mpi_integer', &
638 allNumHeader, 1, 'mpi_integer', 'GRID', ierr)
639
640 if (propagateDSLO) then
641 call aer_propagateDSLO(stateVectorTrlDSLO, inputFileName, outputFileName, &
642 hoursSinceLastAnalysis, hco_ptr)
643 else
644 write(*,*) 'aer_daysSinceLastObs: DSLO trial field is read from: ', trim(inputFileName)
645 call gio_readFromFile(stateVectorTrlDSLO, inputFileName, &
646 etiket_in = ' ', typvar_in = 'P@')
647 end if
648
649 leadTimeInHours = real(stateVectorTrlDSLO%deet * stateVectorTrlDSLO%npasList(1), 8) / 3600.0d0
650 call incdatr(stateVectorAnlDSLO%dateOriginList(1), stateVectorTrlDSLO%dateOriginList(1), &
651 leadTimeInHours)
652
653 call gsv_copyMask(stateVectorTrlDSLO, stateVectorAnlDSLO)
654 stateVectorAnlDSLO%etiket = stateVectorTrlDSLO%etiket
655
656 call col_setVco(column, vco_ptr)
657 call col_allocate(column, obs_numHeader(obsSpaceData), varNames_opt=(/'DSLO'/))
658 call col_setVco(columng, vco_ptr)
659 call col_allocate(columng, obs_numHeader(obsSpaceData), varNames_opt=(/'DSLO'/))
660 call s2c_tl(stateVectorTrlDSLO, column, columng, obsSpaceData)
661
662 call gsv_getField(stateVectorTrlDSLO, trlDSLO_ptr, 'DSLO')
663 call gsv_getField(stateVectorAnlDSLO, anlDSLO_ptr, 'DSLO')
664
665 ! Initialisation
666 do stepIndex = 1, stateVectorTrlDSLO%numStep
667 do levIndex = 1, gsv_getNumLev(stateVectorTrlDSLO, vnl_varLevelFromVarname('DSLO'))
668 do latIndex = 1, hco_ptr%nj
669 do lonIndex = 1, hco_ptr%ni
670 anlDSLO_ptr(lonIndex, latIndex, levIndex, stepIndex) = &
671 trlDSLO_ptr(lonIndex, latIndex, levIndex, stepIndex)
672 end do
673 end do
674 end do
675 end do
676
677 numHeaderMax = maxval(allNumHeader)
678 allocate(obsAss(numHeaderMax))
679 allocate(obsAssMpiGlobal(numHeaderMax, mmpi_nprocs))
680 do headerIndex = 1, obs_numHeader(obsSpaceData)
681 obsAss(headerIndex) = obs_notassimilated
682 bodyIndexBeg = obs_headElem_i(obsSpaceData, obs_rln, headerIndex)
683 bodyIndexEnd = obs_headElem_i(obsSpaceData, obs_nlv, headerIndex) + bodyIndexBeg - 1
684
685 do bodyIndex = bodyIndexBeg, bodyIndexEnd
686 if (obs_bodyElem_i(obsSpaceData, obs_ass, bodyIndex) == obs_assimilated) then
687 obsAss(headerIndex) = obs_assimilated
688 end if
689 end do
690 end do
691 call rpn_comm_gather(obsAss, numHeaderMax, 'mpi_integer', &
692 obsAssMpiGlobal, numHeaderMax, 'mpi_integer', 0, 'grid', ierr)
693
694 if (mmpi_myid == 0) then
695 do procIndex = 1, mmpi_nprocs
696 HEADER_LOOP: do headerIndex = 1, allNumHeader(procIndex)
697
698 if (obsAssMpiGlobal(headerIndex,procIndex) /= obs_assimilated) then
699 cycle HEADER_LOOP
700 end if
701
702 do kIndex = stateVectorTrlDSLO%mykBeg, stateVectorTrlDSLO%mykEnd
703 do stepIndex = 1, stateVectorTrlDSLO%numStep
704 call s2c_getWeightsAndGridPointIndexes(headerIndex, kIndex, stepIndex, &
705 procIndex, interpWeight, obsLatIndex, &
706 obsLonIndex, gridptCount)
707
708 GRIDPT_LOOP: do gridpt = 1, gridptCount
709
710 if (interpWeight(gridpt) == 0.0d0) cycle GRIDPT_LOOP
711
712 lonIndex = obsLonIndex(gridpt)
713 latIndex = obsLatIndex(gridpt)
714
715 do levIndex = 1, gsv_getNumLev(stateVectorTrlDSLO, vnl_varLevelFromVarname('DSLO'))
716 anlDSLO_ptr(lonIndex, latIndex, levIndex, stepIndex) = 0.0
717 end do
718
719 end do GRIDPT_LOOP
720 end do
721 end do
722
723 end do HEADER_LOOP
724 end do
725
726 call gio_writeToFile(stateVectorAnlDSLO, outputFileName, '', typvar_opt = 'A@', &
727 containsFullField_opt = .false. )
728 end if
729
730 deallocate(obsAss)
731 deallocate(obsAssMpiGlobal)
732 call col_deallocate(columng)
733 call col_deallocate(column)
734 call gsv_deallocate(stateVectorTrlDSLO)
735 call gsv_deallocate(stateVectorAnlDSLO)
736
737 end subroutine aer_daysSinceLastObs
738
739 !---------------------------------------------------------
740 ! aer_propagateAnalysisError
741 !---------------------------------------------------------
742 subroutine aer_propagateAnalysisError(stateVectorErrorStd, oceanMask, variableName, &
743 analysisEtiket, errorGrowth, hco_ptr, vco_ptr)
744 !
745 !:Purpose: read analysis error standard deviation field and propagate it forward in time
746 !
747 implicit none
748
749 ! Arguments:
750 type(struct_gsv), intent(inout) :: stateVectorErrorStd ! Input: analysis std error; Output: background std error.
751 type(struct_ocm), intent(in) :: oceanMask ! ocean-land mask (1=water, 0=land)
752 character(len=*), intent(in) :: variableName ! variable name
753 character(len=*), intent(in) :: analysisEtiket ! analysis etiket in the input std file
754 real(4) , intent(in) :: errorGrowth ! seaice: fraction of ice per hour, SST: estimated growth
755 type(struct_hco), pointer, intent(in) :: hco_ptr ! horizontal coordinates structure, pointer
756 type(struct_vco), pointer, intent(in) :: vco_ptr ! vertical coordinates structure, pointer
757
758 ! Locals:
759 type(struct_gsv) :: stateVectorAnalysis
760 integer :: latIndex, lonIndex, localLatIndex, localLonIndex, pointCount
761 real(8), pointer :: stateVectorStdError_ptr(:,:,:), stateVectorAnalysis_ptr(:,:,:)
762 real(4) :: totalLocalVariance
763
764 write(*,*) ''
765 write(*,*) 'aer_propagateAnalysisError: propagate analysis error forward in time for: ', &
766 trim(variableName)
767
768 ! read analysis itself (seaice concentration or SST analysis)
769 call msg('aer_propagateAnalysisError:', ' reading analysis field...')
770 call gsv_allocate(stateVectorAnalysis, 1, hco_ptr, vco_ptr, dateStamp_opt = -1, &
771 mpi_local_opt = .false., mpi_distribution_opt = 'None', &
772 varNames_opt = (/trim(variableName)/), dataKind_opt = 8)
773 call gio_readFromFile(stateVectorAnalysis, './anlm_000m', analysisEtiket, &
774 'A@', containsFullField_opt = .true.)
775 call gsv_getField(stateVectorAnalysis, stateVectorAnalysis_ptr)
776
777 ! initialize pointer for the error standard deviation field
778 call gsv_getField(stateVectorErrorStd, stateVectorStdError_ptr)
779
780 ! calculation in variance unit
781 stateVectorStdError_ptr(:,:,1) = stateVectorStdError_ptr(:,:,1)**2
782
783 pointCount = 0
784 totalLocalVariance = 0.0d0
785
786 do latIndex = 1, hco_ptr%nj
787 do lonIndex = 1, hco_ptr%ni
788
789 OCEANPOINTS: if (oceanMask%mask(lonIndex, latIndex, 1)) then
790 do localLatIndex = latIndex - 1, latIndex + 1
791 if (localLatIndex >= 1 .and. localLatIndex <= hco_ptr%nj) then
792 do localLonIndex = lonIndex - 1, lonIndex + 1
793 if (localLonIndex >= 1 .and. localLonIndex <= hco_ptr%ni) then
794 if (oceanMask%mask(localLonIndex, localLatIndex, 1) .and. &
795 (localLonIndex /= lonIndex .or. localLatIndex /= latIndex)) then
796
797 pointCount = pointCount + 1
798 totalLocalVariance = totalLocalVariance + &
799 (stateVectorAnalysis_ptr(localLonIndex, localLatIndex, 1) - &
800 stateVectorAnalysis_ptr(lonIndex, latIndex, 1))**2
801
802 end if
803 end if
804 end do
805 end if
806 end do
807 end if OCEANPOINTS
808
809 if (pointCount > 0) then
810 totalLocalVariance = totalLocalVariance / real(pointCount)
811 if (stateVectorStdError_ptr(lonIndex, latIndex,1) < totalLocalVariance) then
812 stateVectorStdError_ptr(lonIndex, latIndex, 1) = totalLocalVariance
813 end if
814 end if
815
816 end do
817 end do
818
819 call gsv_deallocate(stateVectorAnalysis)
820
821 ! Back to standard deviation
822 stateVectorStdError_ptr(:, :, 1) = sqrt(stateVectorStdError_ptr(:, :, 1))
823
824 ! Add the standard deviation increment
825 do latIndex = 1, hco_ptr%nj
826 do lonIndex = 1, hco_ptr%ni
827 if (trim(variableName) == 'GL') then
828 stateVectorStdError_ptr(lonIndex, latIndex, 1) = &
829 min(stateVectorStdError_ptr(lonIndex, latIndex, 1) + &
830 errorGrowth * tim_dstepobs, 1.0)
831 else if (trim(variableName) == 'TM') then
832 stateVectorStdError_ptr(lonIndex, latIndex, 1) = &
833 stateVectorStdError_ptr(lonIndex, latIndex, 1) + &
834 errorGrowth * tim_dstepobs / 2.0
835 end if
836 end do
837 end do
838
839 end subroutine aer_propagateAnalysisError
840
841 !---------------------------------------------------------
842 ! aer_propagateDSLO
843 !---------------------------------------------------------
844 subroutine aer_propagateDSLO(stateVectorErrorStd, inputFileName, outputFileName, &
845 hoursSinceLastAnalysis, hco_ptr)
846 !
847 !:Purpose: propagate the field "days since last obs" in time.
848 !
849 implicit none
850
851 ! Arguments:
852 type(struct_gsv), intent(inout) :: stateVectorErrorStd ! read "analysis" into it and increase by hoursSinceLastAnalysis
853 character(len=*), intent(in) :: inputFileName ! input file name
854 character(len=*), intent(in) :: outputFileName ! output file name
855 integer , intent(in) :: hoursSinceLastAnalysis ! hours since last analysis (namelist variable)
856 type(struct_hco), pointer, intent(in) :: hco_ptr ! horizontal grid definition
857
858 ! Locals:
859 real(8), pointer :: analysisDaysSinceLastObs_ptr(:,:,:)
860 real(8) :: daysSinceLastAnalysis
861 integer :: latIndex, lonIndex
862
863 write(*,*) 'aer_propagateDSLO: propagating in time the Days Since Last Obs field...'
864 write(*,*) 'aer_propagateDSLO: hours since last analysis: ', hoursSinceLastAnalysis
865
866 daysSinceLastAnalysis = real(hoursSinceLastAnalysis, 8) / 24.0d0
867
868 ! get DSLO field
869 call gio_readFromFile(stateVectorErrorStd, inputFileName, '', 'A@')
870 call gsv_getField(stateVectorErrorStd, analysisDaysSinceLastObs_ptr)
871
872 ! Add number of days since last obs
873 do latIndex = 1, hco_ptr%nj
874 do lonIndex = 1, hco_ptr%ni
875 analysisDaysSinceLastObs_ptr(lonIndex, latIndex, 1) = analysisDaysSinceLastObs_ptr(lonIndex, latIndex, 1) + &
876 daysSinceLastAnalysis
877 end do
878 end do
879
880 ! update dateStamp from env variable
881 call gsv_modifyDate(stateVectorErrorStd, tim_getDateStamp(), &
882 modifyDateOrigin_opt = .true.)
883
884 ! save increased DSLO field into an fst-file
885 call gio_writeToFile(stateVectorErrorStd, outputFileName, &
886 '', typvar_opt = 'P@', &
887 containsFullField_opt = .false.)
888
889 end subroutine aer_propagateDSLO
890
891 !---------------------------------------------------------
892 ! aer_computeAnlErrorStd
893 !---------------------------------------------------------
894 subroutine aer_computeAnlErrorStd(obsSpaceData, stateVectorAnlErrorStd, stateVectorTrlErrorStd, &
895 analysisVariable, maxAnalysisErrorStdDev, influentObs, Lcorr)
896 !
897 !:Purpose: compute analysis error standard deviation for
898 ! one analysis variable (grid point) at a time
899 !
900 implicit none
901
902 ! Arguments:
903 type(struct_obs) , intent(in) :: obsSpaceData ! obsSpaceData structure
904 type(struct_gsv) , intent(inout) :: stateVectorAnlErrorStd ! state vector for analysis error std deviation
905 type(struct_gsv) , intent(in) :: stateVectorTrlErrorStd ! state vector for background error std deviation
906 character(len=*) , intent(in) :: analysisVariable ! variable name ('GL' or 'TM')
907 real(8) , intent(in) :: maxAnalysisErrorStdDev ! maximum limit imposed on analysis error stddev
908 type(struct_neighborhood), pointer, intent(in) :: influentObs(:,:) ! details about observations to use in update
909 real(8) , intent(in) :: Lcorr(:,:) ! horizontal background-error length scale
910
911 ! Locals:
912 integer :: latIndex, lonIndex, stepIndex, levIndex, kIndex, numStep, numLev
913 integer :: numInfluentObs, bodyIndex, numVariables
914 integer :: influentObsIndex2, influentObsIndex
915 integer :: xStateIndex, yStateIndex, procIndex
916 integer :: xIndex1, yIndex1, xIndex2, yIndex2
917 integer :: gridptCount, gridpt, headerIndex
918 integer :: varIndex1, varIndex2, currentAnalVarIndex, ni, nj
919 real(8), allocatable :: latInRad(:,:), lonInRad(:,:)
920 real(8), pointer :: anlErrorStdDev_ptr(:,:,:,:) ! pointer for analysis error std deviation
921 real(8), pointer :: trlErrorStdDev_ptr(:,:,:,:) ! pointer for background error std deviation
922 real(8), allocatable :: obsOperator(:,:), Bmatrix(:,:), PHiA(:)
923 real(8), allocatable :: innovCovariance(:,:), obsErrorVariance(:)
924 real(8), allocatable :: PH(:,:), KH(:), IKH(:), innovCovarianceInverse(:,:)
925 integer, parameter :: maxvar = 15000
926 integer :: statei(maxvar), statej(maxvar) ! Model grid coordinate i,j of neighbors
927 real(8) :: scaling, distance
928 logical :: found
929 real(8) :: interpWeight(maxNumLocalGridptsSearch)
930 integer :: obsLatIndex(maxNumLocalGridptsSearch), obsLonIndex(maxNumLocalGridptsSearch)
931 integer :: ierr
932 integer :: numBodyMax, allNumBody(mmpi_nprocs)
933 real(8), allocatable :: allObsOer(:,:), obsOer(:), iceScaling(:), allIceScaling(:,:), anlErrorStdDevMpiGlobal(:,:,:,:)
934
935 call msg('aer_computeAnlErrorStd:', ' computing analysis error std field initialization...')
936
937 call utl_tmg_start(124,'--AnalErrOI_computeAnlErrStd')
938
939 ni = stateVectorTrlErrorStd%hco%ni
940 nj = stateVectorTrlErrorStd%hco%nj
941
942 ! compute lon/lat in radians
943 allocate(lonInRad(ni, nj))
944 allocate(latInRad(ni, nj))
945
946 do latIndex = 1, nj
947 do lonIndex = 1, ni
948 latInRad(lonIndex, latIndex) = real(stateVectorTrlErrorStd%hco%lat2d_4(lonIndex, latIndex), 8)
949 lonInRad(lonIndex, latIndex) = real(stateVectorTrlErrorStd%hco%lon2d_4(lonIndex, latIndex), 8)
950 end do
951 end do
952
953 ! Initialisation of pointers
954 call gsv_getField(stateVectorAnlErrorStd, anlErrorStdDev_ptr)
955 call gsv_getField(stateVectorTrlErrorStd, trlErrorStdDev_ptr)
956 anlErrorStdDev_ptr(:,:,:,:) = 0.0d0
957
958 ! Communicate some values from proc 0 to all others
959 call rpn_comm_allgather(obs_numBody(obsSpaceData), 1, 'mpi_integer', &
960 allNumBody, 1, 'mpi_integer', 'GRID', ierr)
961 numBodyMax = maxval(allNumBody)
962 write(*,*) 'aer_computeAnlErrorStd: numBodyMax = ', numBodyMax
963
964 if (analysisVariable == 'GL') then
965 allocate(iceScaling(numBodyMax))
966 allocate(allIceScaling(numBodyMax,mmpi_nprocs))
967 do bodyIndex = 1, obs_numBody(obsSpaceData)
968 iceScaling(bodyIndex) = oop_iceScaling(obsSpaceData, bodyIndex)
969 end do
970 call rpn_comm_allGather(iceScaling, numBodyMax, 'MPI_REAL8', &
971 allIceScaling, numBodyMax, 'MPI_REAL8', 'GRID', ierr)
972 end if
973
974 allocate(obsOer(numBodyMax))
975 allocate(allObsOer(numBodyMax,mmpi_nprocs))
976 do bodyIndex = 1, obs_numBody(obsSpaceData)
977 obsOer(bodyIndex) = obs_bodyElem_r(obsSpaceData, OBS_OER, &
978 bodyIndex)
979 end do
980 call rpn_comm_allGather(obsOer, numBodyMax, 'MPI_REAL8', &
981 allObsOer, numBodyMax, 'MPI_REAL8', 'GRID', ierr)
982
983 numStep = stateVectorTrlErrorStd%numStep
984 numLev = gsv_getNumLev(stateVectorTrlErrorStd, vnl_varLevelFromVarname(analysisVariable))
985
986 ! Only variables assigned within or by the loop can be private.
987 STEP: do stepIndex = 1, numStep
988 LEVEL: do levIndex = 1, numLev
989
990 !$omp parallel do default(shared) schedule(dynamic) private(obsOperator, Bmatrix, &
991 !$omp PHiA, innovCovariance, innovCovarianceInverse, obsErrorVariance, &
992 !$omp PH, KH, IKH, &
993 !$omp statei, statej, &
994 !$omp numVariables, varIndex1, varIndex2, currentAnalVarIndex, found, &
995 !$omp influentObsIndex2, influentObsIndex, &
996 !$omp scaling, xStateIndex, yStateIndex, lonIndex, latIndex, headerIndex, &
997 !$omp bodyIndex, kIndex, procIndex, &
998 !$omp numInfluentObs, xIndex1, yIndex1, xIndex2, yIndex2, distance, &
999 !$omp interpWeight, obsLatIndex, obsLonIndex, gridptCount, gridpt)
1000 YINDEX: do latIndex = mmpi_myidy+1, nj, mmpi_npey
1001 XINDEX: do lonIndex = mmpi_myidx+1, ni, mmpi_npex
1002
1003 ! Initialize to trial value
1004 anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex) = &
1005 trlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex)
1006
1007 numInfluentObs = influentObs(lonIndex, latIndex)%numObs
1008
1009 if (numInfluentObs == 0 .or. &
1010 .not. stateVectorTrlErrorStd%oceanMask%mask(lonIndex, latIndex, levIndex)) cycle XINDEX
1011
1012 ! form the observation-error covariance (diagonal) matrix
1013
1014 allocate(obsErrorVariance(numInfluentObs))
1015
1016 do influentObsIndex = 1, numInfluentObs
1017 bodyIndex = influentObs(lonIndex, latIndex)%bodyIndex(influentObsIndex)
1018 procIndex = influentObs(lonIndex, latIndex)%procIndex(influentObsIndex)
1019 obsErrorVariance(influentObsIndex) = allObsOer(bodyIndex,procIndex)**2
1020 end do
1021
1022 ! find all model variables involved here
1023
1024 numVariables = 0
1025 statei(:) = 0
1026 statej(:) = 0
1027
1028 INFLUENTOBSCYCLE: do influentObsIndex = 1, numInfluentObs
1029
1030 headerIndex = influentObs(lonIndex, latIndex)%headerIndex(influentObsIndex)
1031 bodyIndex = influentObs(lonIndex, latIndex)%bodyIndex(influentObsIndex)
1032 procIndex = influentObs(lonIndex, latIndex)%procIndex(influentObsIndex)
1033
1034 if (analysisVariable == 'GL') then
1035 scaling = allIceScaling(bodyIndex,procIndex)
1036 else if (analysisVariable == 'TM') then
1037 scaling = 1.0d0
1038 end if
1039
1040 if (scaling == 0.0d0) cycle INFLUENTOBSCYCLE
1041 KINDEXCYCLE: do kIndex = stateVectorTrlErrorStd%mykBeg, stateVectorTrlErrorStd%mykEnd
1042
1043 call s2c_getWeightsAndGridPointIndexes(headerIndex, kIndex, stepIndex, procIndex, &
1044 interpWeight, obsLatIndex, obsLonIndex, &
1045 gridptCount)
1046
1047 GRIDPTCYCLE: do gridpt = 1, gridptCount
1048
1049 if (interpWeight(gridpt) == 0.0d0) cycle GRIDPTCYCLE
1050
1051 xStateIndex = obsLonIndex(gridpt)
1052 yStateIndex = obsLatIndex(gridpt)
1053
1054 found = .false.
1055 do varIndex2=1, numVariables
1056 if(xStateIndex == statei(varIndex2) .and. &
1057 yStateIndex == statej(varIndex2)) then
1058 found = .true.
1059 exit
1060 end if
1061 end do
1062
1063 if (found) cycle GRIDPTCYCLE
1064
1065 numVariables = numVariables + 1
1066 if(numVariables > maxvar) then
1067 call utl_abort('aer_computeAnlErrorStd: Structure state'// &
1068 ' too small in subroutine'// &
1069 ' analysis_error_mod. Increase maxvar')
1070 end if
1071 statei(numVariables) = xStateIndex
1072 statej(numVariables) = yStateIndex
1073
1074 end do GRIDPTCYCLE
1075 end do KINDEXCYCLE
1076 end do INFLUENTOBSCYCLE
1077
1078 ! make sure that current analysis variable is part of the state
1079 ! vector even if it does not participate in obs calculation
1080
1081 found = .false.
1082 do varIndex2 = 1, numVariables
1083 if(lonIndex == statei(varIndex2) .and. &
1084 latIndex == statej(varIndex2)) then
1085 currentAnalVarIndex = varIndex2
1086 found = .true.
1087 exit
1088 end if
1089 end do
1090
1091 if(.not. found) then
1092 numVariables = numVariables + 1
1093 if(numVariables > maxvar) then
1094 call utl_abort('aer_computeAnlErrorStd: Structure state too small in'// &
1095 ' subroutine analysis_error_mod'// &
1096 ' increase maxvar')
1097 end if
1098 currentAnalVarIndex = numVariables
1099 statei(numVariables) = lonIndex
1100 statej(numVariables) = latIndex
1101 end if
1102
1103 ! form the observation operator matrix (obsOperator)
1104
1105 allocate(obsOperator(numVariables, numInfluentObs))
1106
1107 obsOperator(:,:) = 0.0d0
1108
1109 varIndex2 = 0
1110
1111 INFLUENTOBSCYCLE2: do influentObsIndex = 1, numInfluentObs
1112
1113 headerIndex = influentObs(lonIndex, latIndex)%headerIndex(influentObsIndex)
1114 bodyIndex = influentObs(lonIndex, latIndex)%bodyIndex(influentObsIndex)
1115 procIndex = influentObs(lonIndex, latIndex)%procIndex(influentObsIndex)
1116
1117 if (analysisVariable == 'GL') then
1118 scaling = allIceScaling(bodyIndex,procIndex)
1119 else if (analysisVariable == 'TM') then
1120 scaling = 1.0d0
1121 end if
1122
1123 if (scaling == 0.0d0) cycle INFLUENTOBSCYCLE2
1124
1125 KINDEXCYCLE2: do kIndex = stateVectorTrlErrorStd%mykBeg, stateVectorTrlErrorStd%mykEnd
1126
1127 call s2c_getWeightsAndGridPointIndexes(headerIndex, kIndex, stepIndex, procIndex, &
1128 interpWeight, obsLatIndex, obsLonIndex, &
1129 gridptCount)
1130
1131 GRIDPTCYCLE2: do gridpt = 1, gridptCount
1132
1133 if (interpWeight(gridpt) == 0.0d0) cycle GRIDPTCYCLE2
1134
1135 xStateIndex = obsLonIndex(gridpt)
1136 yStateIndex= obsLatIndex(gridpt)
1137
1138 found = .false.
1139 do while (.not. found)
1140 if(varIndex2 < numVariables) then
1141 varIndex2 = varIndex2 + 1
1142 else
1143 varIndex2 = 1
1144 end if
1145 if(xStateIndex == statei(varIndex2) .and. &
1146 yStateIndex == statej(varIndex2)) then
1147 found = .true.
1148 obsOperator(varIndex2, influentObsIndex) = scaling * &
1149 interpWeight(gridpt)
1150 end if
1151 end do
1152
1153 if(found) cycle GRIDPTCYCLE2
1154
1155 write(*,*) 'xStateIndex = ', xStateIndex
1156 write(*,*) 'yStateIndex = ', yStateIndex
1157 write(*,*) 'lonIndex = ', lonIndex
1158 write(*,*) 'latIndex = ', latIndex
1159 write(*,*) 'gridptCount = ', gridptCount
1160 write(*,*) 'gridpt = ', gridpt
1161 write(*,*) 'numVariables = ', numVariables
1162 do varIndex2=1, numVariables
1163 write(*,*) 'varIndex2 statei statej = ',statei(varIndex2), &
1164 statej(varIndex2)
1165 end do
1166 call utl_abort('aer_computeAnlErrorStd: not found in state vect.')
1167
1168 end do GRIDPTCYCLE2
1169 end do KINDEXCYCLE2
1170 end do INFLUENTOBSCYCLE2
1171
1172 ! form the background-error covariance matrix
1173
1174 allocate(Bmatrix(numVariables, numVariables))
1175
1176 Bmatrix(:,:) = 0.0d0
1177
1178 do varIndex2 = 1, numVariables
1179
1180 xIndex2 = statei(varIndex2)
1181 yIndex2 = statej(varIndex2)
1182
1183 do varIndex1 = varIndex2, numVariables
1184
1185 xIndex1 = statei(varIndex1)
1186 yIndex1 = statej(varIndex1)
1187
1188 if(xIndex2 == xIndex1 .and. yIndex2 == yIndex1) then
1189 Bmatrix(varIndex1,varIndex2) = trlErrorStdDev_ptr(xIndex1, yIndex1, levIndex, stepIndex)**2
1190 else if(Lcorr(xIndex2,yIndex2) > 0.0d0) then
1191 distance = phf_calcDistance(latInRad(xIndex2, yIndex2), lonInRad(xIndex2, yIndex2), &
1192 latInRad(xIndex1, yIndex1), lonInRad(xIndex1, yIndex1))
1193 Bmatrix(varIndex1,varIndex2) = trlErrorStdDev_ptr(xIndex2, yIndex2, levIndex, stepIndex) * &
1194 trlErrorStdDev_ptr(xIndex1, yIndex1, levIndex, stepIndex) * &
1195 exp(-0.5 * (distance / Lcorr(xIndex2, yIndex2))**2)
1196 end if
1197
1198 ! symmetric matrix !
1199 Bmatrix(varIndex2, varIndex1) = Bmatrix(varIndex1, varIndex2)
1200
1201 end do
1202
1203 end do
1204
1205 ! form the observation background error covariance matrix (PHT)
1206
1207 allocate(PH(numInfluentObs, numVariables))
1208
1209 ! PH = matmul (Bmatrix, transpose(obsOperator))
1210 do varIndex1 = 1, numVariables
1211 do influentObsIndex = 1, numInfluentObs
1212 PH(influentObsIndex,varIndex1) = dot_product(Bmatrix(:, varIndex1), &
1213 obsOperator(:, influentObsIndex))
1214 end do
1215 end do
1216
1217 ! covariance matrix of the innovation (HPHT + R)
1218
1219 allocate(innovCovariance(numInfluentObs, numInfluentObs), &
1220 innovCovarianceInverse(numInfluentObs, numInfluentObs))
1221
1222 do influentObsIndex = 1, numInfluentObs
1223
1224 do influentObsIndex2 = 1, numInfluentObs
1225 if (influentObsIndex == influentObsIndex2) then
1226 innovCovariance(influentObsIndex2, influentObsIndex) = &
1227 dot_product(obsOperator(:, influentObsIndex), &
1228 PH(influentObsIndex2,:)) + &
1229 obsErrorVariance(influentObsIndex)
1230 else
1231 innovCovariance(influentObsIndex2, influentObsIndex) = &
1232 dot_product(obsOperator(:, influentObsIndex), &
1233 PH(influentObsIndex2,:))
1234 end if
1235 end do
1236
1237 end do
1238
1239 ! Inverse of the covariance matrix of the innovation
1240
1241 innovCovarianceInverse(:,:) = innovCovariance(:,:)
1242 call utl_matInverse(innovCovarianceInverse, numInfluentObs, &
1243 printInformation_opt = .false.)
1244
1245 ! Kalman gain; this is the row corresponding to the analysis variable
1246
1247 allocate(PHiA(numInfluentObs))
1248
1249 do influentObsIndex = 1, numInfluentObs
1250 PHiA(influentObsIndex) = dot_product(PH(:, currentAnalVarIndex), &
1251 innovCovarianceInverse(influentObsIndex,:))
1252 end do
1253
1254 ! compute the error variance of the analysis
1255
1256 allocate(KH(numVariables))
1257
1258 ! KH = matmul (PHiA, obsOperator)
1259 do varIndex1 = 1, numVariables
1260 KH(varIndex1) = dot_product(PHiA, obsOperator(varIndex1,:))
1261 end do
1262
1263 ! IKH = I - KH
1264
1265 allocate(IKH(numVariables))
1266
1267 IKH = -KH
1268
1269 IKH(currentAnalVarIndex) = 1.0d0 - KH(currentAnalVarIndex)
1270
1271 anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex) = &
1272 sqrt(dot_product (IKH, Bmatrix(:, currentAnalVarIndex)))
1273
1274 if(anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex) < 0.0) then
1275 write(*,*) 'aer_computeAnlErrorStd: negative analysis-error Std dev. = ', &
1276 anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex), &
1277 ' reset to zero at grid point (',lonIndex, latIndex,')'
1278 anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex) = &
1279 max(anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex), 0.0)
1280 end if
1281
1282 if(anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex) > &
1283 trlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex)) then
1284 write(*,*) 'aer_computeAnlErrorStd: analysis-error Std dev. = ', &
1285 anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex), &
1286 ' is larger than background-error ', &
1287 'Std dev. is kept at = ', &
1288 trlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex)
1289 anlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex) = &
1290 min(trlErrorStdDev_ptr(lonIndex, latIndex, levIndex, stepIndex), &
1291 maxAnalysisErrorStdDev)
1292 end if
1293
1294 deallocate(Bmatrix)
1295 deallocate(obsErrorVariance)
1296 deallocate(innovCovariance)
1297 deallocate(innovCovarianceInverse)
1298 deallocate(obsOperator)
1299 deallocate(PH)
1300 deallocate(PHiA)
1301 deallocate(KH)
1302 deallocate(IKH)
1303 deallocate(influentObs(lonIndex, latIndex)%headerIndex)
1304 deallocate(influentObs(lonIndex, latIndex)%bodyIndex)
1305 deallocate(influentObs(lonIndex, latIndex)%procIndex)
1306
1307 end do XINDEX
1308 end do YINDEX
1309 !$omp end parallel do
1310
1311 end do LEVEL
1312 end do STEP
1313
1314 deallocate(lonInRad)
1315 deallocate(latInRad)
1316
1317 ! Combine results from all MPI tasks onto task 0
1318 if (mmpi_nprocs > 1) then
1319 write(*,*) 'aer_computeAnlErrorStd: do mpi communication of anlErrorStdDev'
1320 allocate(anlErrorStdDevMpiGlobal(ni,nj,numStep,numLev))
1321 call rpn_comm_reduce(anlErrorStdDev_ptr, anlErrorStdDevMpiGlobal, &
1322 size(anlErrorStdDev_ptr), "mpi_real8", "MPI_SUM", 0, "GRID", ierr)
1323 anlErrorStdDev_ptr(:,:,:,:) = anlErrorStdDevMpiGlobal(:,:,:,:)
1324 deallocate(anlErrorStdDevMpiGlobal)
1325 end if
1326
1327 call utl_tmg_stop(124)
1328
1329 write(*,*) 'aer_computeAnlErrorStd: finished'
1330
1331 end subroutine aer_computeAnlErrorStd
1332
1333end module analysisErrorOI_mod
