1module enkf_mod
2 ! MODULE enkf_mod (prefix='enkf' category='1. High-level functionality')
3 !
4 !:Purpose: Various routines that are useful for implementing
5 ! an EnKF in MIDAS, including the LETKF.
6 !
7 use mpi, only : mpi_statuses_ignore ! this is the mpi library module
8 use midasMpi_mod
9 use utilities_mod
10 use mathPhysConstants_mod
11 use timeCoord_mod
12 use verticalCoord_mod
13 use horizontalCoord_mod
14 use ensembleStateVector_mod
15 use gridStateVector_mod
16 use obsSpaceData_mod
17 use tovsNL_mod
18 use ensembleObservations_mod
19 use gridVariableTransforms_mod
20 use localizationFunction_mod
21 use varNameList_mod
22 use codePrecision_mod
23 use codTyp_mod
24 use calcHeightAndPressure_mod
25 implicit none
26 save
27 private
28
29 ! public types
30 public :: struct_enkfInterpInfo
31
32 ! public procedures
33 public :: enkf_setupInterpInfo, enkf_LETKFanalyses, enkf_modifyAMSUBobsError
34 public :: enkf_rejectHighLatIR, enkf_getModulatedState, enkf_setupModulationFactor
35
36 ! for weight interpolation
37 type struct_enkfInterpInfo
38 integer :: latLonStep
39 integer, allocatable :: numIndexes(:,:)
40 integer, allocatable :: lonIndexes(:,:,:)
41 integer, allocatable :: latIndexes(:,:,:)
42 real(8), allocatable :: interpWeights(:,:,:)
43 integer :: myLonBegHalo
44 integer :: myLonEndHalo
45 integer :: myLatBegHalo
46 integer :: myLatEndHalo
47 integer :: myLonBeg
48 integer :: myLonEnd
49 integer :: myLatBeg
50 integer :: myLatEnd
51 end type struct_enkfInterpInfo
52
53 integer, external :: get_max_rss
54
55contains
56
57 !----------------------------------------------------------------------
58 ! enkf_LETKFanalyses
59 !----------------------------------------------------------------------
60 subroutine enkf_LETKFanalyses(algorithm, numSubEns, randomShuffleSubEns, &
61 ensembleAnl, ensembleTrl, &
62 ensObs_mpiglobal, ensObsGain_mpiglobal, &
63 stateVectorMeanAnl, &
64 wInterpInfo, maxNumLocalObs, &
65 hLocalize, hLocalizePressure, vLocalize, &
66 mpiDistribution, numRetainedEigen)
67 !
68 !:Purpose: Local subroutine containing the code for computing
69 ! the LETKF analyses for all ensemble members, ensemble
70 ! mean.
71 !
72 implicit none
73
74 ! Arguments:
75 character(len=*), intent(in) :: algorithm
76 integer , intent(in) :: numSubEns
77 logical , intent(in) :: randomShuffleSubEns
78 type(struct_ens), pointer, intent(inout) :: ensembleTrl
79 type(struct_ens), intent(inout) :: ensembleAnl
80 type(struct_eob), target, intent(in) :: ensObs_mpiglobal
81 type(struct_eob), intent(in) :: ensObsGain_mpiglobal
82 type(struct_gsv), intent(in) :: stateVectorMeanAnl
83 type(struct_enkfInterpInfo), intent(in) :: wInterpInfo
84 integer, intent(in) :: maxNumLocalObs
85 real(8), intent(in) :: hLocalize(:)
86 real(8), intent(in) :: hLocalizePressure(:)
87 real(8), intent(in) :: vLocalize
88 character(len=*), intent(in) :: mpiDistribution
89 integer, intent(in) :: numRetainedEigen
90
91 ! Locals:
92 integer :: nEns, nEnsPerSubEns, nEnsPerSubEns_mod, nEnsIndependentPerSubEns
93 integer :: nLev_M, nLev_depth, nLev_weights
94 integer :: memberIndex, memberIndex1, memberIndex2, ierr, matrixRank
95 integer :: memberIndexCV, memberIndexCV1, memberIndexCV2
96 integer :: procIndex, procIndexSend, hLocIndex
97 integer :: latIndex, lonIndex, stepIndex, varLevIndex, levIndex, levIndex2
98 integer :: bodyIndex, localObsIndex, numLocalObs, numLocalObsFound
99 integer :: countMaxExceeded, maxCountMaxExceeded, numGridPointWeights
100 integer :: myNumLatLonRecv, myNumLatLonSend
101 integer :: latLonIndex, subEnsIndex, subEnsIndex2
102 integer :: sendTag, recvTag, nsize, numRecv, numSend
103 integer :: myLonBeg, myLonEnd, myLatBeg, myLatEnd, numVarLev
104 integer :: myLonBegHalo, myLonEndHalo, myLatBegHalo, myLatEndHalo
105 integer :: imode, dateStamp, timePrint, datePrint, randomSeed, newDate
106 integer :: nEnsGain, eigenVectorColumnIndex
107 integer :: memberIndexInModEns
108 real(8) :: anlLat, anlLon, anlVertLocation
109 real(8) :: distance, tolerance, localization
110 real(4) :: modulationFactor_r4
111
112 integer, allocatable :: localBodyIndices(:), levFromK(:)
113 integer, allocatable :: myLatIndexesRecv(:), myLonIndexesRecv(:)
114 integer, allocatable :: myLatIndexesSend(:), myLonIndexesSend(:)
115 integer, allocatable :: myNumProcIndexesSend(:)
116 integer, allocatable :: myProcIndexesRecv(:), myProcIndexesSend(:,:)
117 integer, allocatable :: requestIdRecv(:), requestIdSend(:)
118 integer, allocatable :: memberIndexSubEns(:,:), memberIndexSubEns_mod(:,:)
119 integer, allocatable :: memberIndexSubEnsComp(:,:)
120 integer, allocatable :: randomMemberIndexArray(:), latLonTagMpiGlobal(:,:)
121
122 real(8), pointer :: PaInv_mean(:,:), Pa_mean(:,:)
123 real(8), pointer :: YbTinvR_mean(:,:), YbTinvRCopy_mean(:,:), YbTinvRYb_mean(:,:)
124 real(8), pointer :: eigenValues_mean(:), eigenVectors_mean(:,:)
125
126 real(8), allocatable, target :: PaInv_pert(:,:), Pa_pert(:,:)
127 real(8), allocatable, target :: YbTinvR_pert(:,:),YbTinvRCopy_pert(:,:), YbTinvRYb_pert(:,:)
128 real(8), allocatable, target :: eigenValues_pert(:), eigenVectors_pert(:,:)
129
130 real(8), allocatable :: distances(:), PaSqrt_pert(:,:)
131 real(8), allocatable :: YbTinvRYb_CV(:,:), YbTinvRYb_mod(:,:)
132 real(8), allocatable :: eigenValues_CV(:), eigenVectors_CV(:,:)
133 real(8), allocatable :: weightsTemp(:), weightsTemp2(:)
134 real(8), allocatable :: weightsMembers(:,:,:,:), weightsMembersLatLon(:,:,:)
135 real(8), allocatable :: weightsMean(:,:,:,:), weightsMeanLatLon(:,:,:)
136 real(8), allocatable :: memberAnlPert(:)
137 real(4), allocatable :: vertLocation_r4(:,:,:), YbCopy_r4(:,:), YbGainCopy_r4(:,:)
138
139 real(4), pointer :: meanTrl_ptr_r4(:,:,:,:), meanAnl_ptr_r4(:,:,:,:), meanInc_ptr_r4(:,:,:,:)
140 real(4), pointer :: memberTrl_ptr_r4(:,:,:,:), memberAnl_ptr_r4(:,:,:,:)
141 real(4) :: pert_r4
142
143 character(len=4) :: varName
144 character(len=2), allocatable :: varKindFromK(:)
145 character(len=4), allocatable :: varLevelFromK(:)
146
147 type(struct_hco), pointer :: hco_ens
148 type(struct_vco), pointer :: vco_ens
149 type(struct_gsv) :: stateVectorMeanInc
150 type(struct_gsv) :: stateVectorMeanTrl
151
152 logical :: hLocalizeIsConstant, useModulatedEns, firstTime = .true.
153
154 call utl_tmg_start(131,'--LETKFanalysis')
155
156 write(*,*) 'enkf_LETKFanalyses: starting'
157 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
158
159 !
160 ! Set things up for the redistribution of work across mpi tasks
161 !
162 call enkf_LETKFsetupMpiDistribution(myNumLatLonRecv, myNumLatLonSend, &
163 myLatIndexesRecv, myLonIndexesRecv, &
164 myLatIndexesSend, myLonIndexesSend, &
165 myProcIndexesRecv, myProcIndexesSend, &
166 myNumProcIndexesSend, mpiDistribution, wInterpInfo)
167 allocate(requestIdSend(3*myNumLatLonSend*maxval(myNumProcIndexesSend)))
168 allocate(requestIdRecv(3*myNumLatLonRecv))
169
170 nEns = ens_getNumMembers(ensembleAnl)
171 useModulatedEns = (numRetainedEigen > 0)
172 if ( useModulatedEns ) then
173 nEnsGain = nEns * numRetainedEigen
174 else
175 nEnsGain = nEns
176 end if
177 nLev_M = ens_getNumLev(ensembleAnl, 'MM')
178 nLev_depth = ens_getNumLev(ensembleAnl, 'DP')
179 nLev_weights = max(nLev_M,nLev_depth)
180 if ( useModulatedEns ) nLev_weights = 1
181 hco_ens => ens_getHco(ensembleAnl)
182 vco_ens => ens_getVco(ensembleAnl)
183 myLonBeg = stateVectorMeanAnl%myLonBeg
184 myLonEnd = stateVectorMeanAnl%myLonEnd
185 myLatBeg = stateVectorMeanAnl%myLatBeg
186 myLatEnd = stateVectorMeanAnl%myLatEnd
187 numVarLev = stateVectorMeanAnl%nk
188 myLonBegHalo = wInterpInfo%myLonBegHalo
189 myLonEndHalo = wInterpInfo%myLonEndHalo
190 myLatBegHalo = wInterpInfo%myLatBegHalo
191 myLatEndHalo = wInterpInfo%myLatEndHalo
192
193 !
194 ! Compute gridded 3D ensemble weights
195 !
196 allocate(localBodyIndices(maxNumLocalObs))
197 allocate(distances(maxNumLocalObs))
198 allocate(YbTinvR_pert(nEnsGain,maxNumLocalObs))
199 allocate(YbTinvRCopy_pert(maxNumLocalObs,nEnsGain))
200 allocate(YbGainCopy_r4(maxNumLocalObs,nEnsGain))
201 allocate(YbTinvRYb_pert(nEnsGain,nEnsGain))
202 if ( trim(algorithm) == 'CVLETKF-ME' .or. &
203 trim(algorithm) == 'LETKF-Gain-ME' ) then
204 allocate(YbTinvRYb_mod(nEnsGain,nEns))
205 allocate(YbCopy_r4(maxNumLocalObs,nEns))
206 end if
207 allocate(eigenValues_pert(nEnsGain))
208 allocate(eigenVectors_pert(nEnsGain,nEnsGain))
209 allocate(PaInv_pert(nEnsGain,nEnsGain))
210 allocate(PaSqrt_pert(nEnsGain,nEnsGain))
211 allocate(Pa_pert(nEnsGain,nEnsGain))
212
213 if (eob_simObsAssim) then
214 allocate(YbTinvR_mean(nEnsGain,maxNumLocalObs))
215 allocate(YbTinvRCopy_mean(maxNumLocalObs,nEnsGain))
216 allocate(YbTinvRYb_mean(nEnsGain,nEnsGain))
217 allocate(eigenValues_mean(nEnsGain))
218 allocate(eigenVectors_mean(nEnsGain,nEnsGain))
219 allocate(PaInv_mean(nEnsGain,nEnsGain))
220 allocate(Pa_mean(nEnsGain,nEnsGain))
221 else
222 YbTinvR_mean => YbTinvR_pert
223 YbTinvRCopy_mean => YbTinvRCopy_pert
224 YbTinvRYb_mean => YbTinvRYb_pert
225 eigenValues_mean => eigenValues_pert
226 eigenVectors_mean => eigenVectors_pert
227 PaInv_mean => PaInv_pert
228 Pa_mean => Pa_pert
229 end if
230
231 allocate(memberAnlPert(nEns))
232 allocate(weightsTemp(nEnsGain))
233 allocate(weightsTemp2(nEnsGain))
234 weightsTemp(:) = 0.0d0
235 weightsTemp2(:) = 0.0d0
236 ! Weights for mean analysis
237 allocate(weightsMean(nEnsGain,1,myLonBegHalo:myLonEndHalo,myLatBegHalo:myLatEndHalo))
238 weightsMean(:,:,:,:) = 0.0d0
239 allocate(weightsMeanLatLon(nEnsGain,1,myNumLatLonSend))
240 weightsMeanLatLon(:,:,:) = 0.0d0
241 ! Weights for member analyses
242 allocate(weightsMembers(nEnsGain,nEns,myLonBegHalo:myLonEndHalo,myLatBegHalo:myLatEndHalo))
243 weightsMembers(:,:,:,:) = 0.0d0
244 allocate(weightsMembersLatLon(nEnsGain,nEns,myNumLatLonSend))
245 weightsMembersLatLon(:,:,:) = 0.0d0
246
247 call gsv_allocate( stateVectorMeanTrl, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
248 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
249 dataKind_opt=4, allocHeightSfc_opt=.true., &
250 allocHeight_opt=.false., allocPressure_opt=.false. )
251 call gsv_zero(stateVectorMeanTrl)
252 call gsv_allocate( stateVectorMeanInc, tim_nstepobsinc, hco_ens, vco_ens, dateStamp_opt=tim_getDateStamp(), &
253 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
254 dataKind_opt=4, allocHeightSfc_opt=.true., &
255 allocHeight_opt=.false., allocPressure_opt=.false. )
256 call gsv_zero(stateVectorMeanInc)
257
258 call ens_computeMean(ensembleTrl)
259 call ens_copyEnsMean(ensembleTrl, stateVectorMeanTrl)
260
261 ! Quantities needed for CVLETKF and CVLETKF-PERTOBS and CVLETKF-ME
262 if ( trim(algorithm) == 'CVLETKF' .or. trim(algorithm) == 'CVLETKF-PERTOBS' .or. &
263 trim(algorithm) == 'CVLETKF-ME' ) then
264 nEnsPerSubEns = nEns / numSubEns
265 if ( (nEnsPerSubEns * numSubEns) /= nEns ) then
266 call utl_abort('enkf_LETKFanalyses: ensemble size not divisible by numSubEnsembles')
267 end if
268 if (numSubEns <= 1) then
269 call utl_abort('enkf_LETKFanalyses: for CVLETKF(-PERTOBS)(-ME) algorithm, numSubEns must be greater than 1')
270 end if
271 if ( .not. useModulatedEns ) then
272 nEnsIndependentPerSubEns = nEns - nEnsPerSubEns
273 else
274 nEnsPerSubEns_mod = nEnsPerSubEns * numRetainedEigen
275 nEnsIndependentPerSubEns = nEnsGain - nEnsPerSubEns_mod
276 end if
277 allocate(YbTinvRYb_CV(nEnsIndependentPerSubEns,nEnsIndependentPerSubEns))
278 allocate(eigenValues_CV(nEnsIndependentPerSubEns))
279 allocate(eigenVectors_CV(nEnsIndependentPerSubEns,nEnsIndependentPerSubEns))
280 allocate(memberIndexSubEns(nEnsPerSubEns,numSubEns))
281 allocate(memberIndexSubEnsComp(nEnsIndependentPerSubEns,numSubEns))
282 if ( useModulatedEns ) allocate(memberIndexSubEns_mod(nEnsPerSubEns_mod,numSubEns))
283 if (.not.randomShuffleSubEns) then
284 ! form subensembles with contiguous sequential groups of members
285 do subEnsIndex = 1, numSubEns
286 do memberIndex = 1, nEnsPerSubEns
287 memberIndexSubEns(memberIndex,subEnsIndex) = &
288 (subEnsIndex-1)*nEnsPerSubEns + memberIndex
289 end do
290 end do
291 if ( useModulatedEns ) then
292 do subEnsIndex = 1, numSubEns
293 memberIndex2 = 0
294 do memberIndex = 1, nEnsPerSubEns
295 do eigenVectorColumnIndex = 1, numRetainedEigen
296 memberIndex2 = memberIndex2 + 1
297 memberIndexInModEns = (eigenVectorColumnIndex - 1) * nEns + &
298 memberIndex
299 memberIndexSubEns_mod(memberIndex2,subEnsIndex) = &
300 (subEnsIndex-1)*nEnsPerSubEns + memberIndexInModEns
301 end do
302 end do
303 end do
304 end if
305 else
306 ! compute random seed from the date for randomly forming subensembles
307 imode = -3 ! stamp to printable date and time: YYYYMMDD, HHMMSShh
308 dateStamp = tim_getDateStamp()
309 ierr = newdate(dateStamp, datePrint, timePrint, imode)
310 timePrint = timePrint/1000000
311 datePrint = datePrint*100 + timePrint
312 ! Remove the century, keeping 2 digits of the year
313 randomSeed = datePrint - 100000000*(datePrint/100000000)
314 allocate(randomMemberIndexArray(nEns))
315 do memberIndex = 1, nEns
316 randomMemberIndexArray(memberIndex) = memberIndex
317 end do
318 call utl_randomOrderInt(randomMemberIndexArray,randomSeed)
319 write(*,*) 'enkf_LETKFanalyses: seed for random shuffle of sub ens = ', randomSeed
320 write(*,*) 'enkf_LETKFanalyses: randomOrder = ', randomMemberIndexArray(:)
321 do subEnsIndex = 1, numSubEns
322 do memberIndex = 1, nEnsPerSubEns
323 memberIndexSubEns(memberIndex,subEnsIndex) = &
324 randomMemberIndexArray((subEnsIndex-1)*nEnsPerSubEns + memberIndex)
325 end do
326 end do
327 if ( useModulatedEns ) then
328 do subEnsIndex = 1, numSubEns
329 memberIndex2 = 0
330 do memberIndex = 1, nEnsPerSubEns
331 do eigenVectorColumnIndex = 1, numRetainedEigen
332 memberIndex2 = memberIndex2 + 1
333 memberIndexSubEns_mod(memberIndex2,subEnsIndex) = &
334 randomMemberIndexArray((subEnsIndex-1)*nEnsPerSubEns + memberIndex) + &
335 (eigenVectorColumnIndex - 1) * nEns
336 end do
337 end do
338 end do
339 end if
340 end if
341
342 do subEnsIndex = 1, numSubEns
343 memberIndex = 1
344 do subEnsIndex2 = 1, numSubEns
345 if (subEnsIndex2 == subEnsIndex) cycle
346
347 if ( .not. useModulatedEns ) then
348 memberIndexSubEnsComp(memberIndex:memberIndex+nEnsPerSubEns-1,subEnsIndex) = &
349 memberIndexSubEns(:,subEnsIndex2)
350 memberIndex = memberIndex + nEnsPerSubEns
351 else
352 memberIndexSubEnsComp(memberIndex:memberIndex+nEnsPerSubEns_mod-1,subEnsIndex) = &
353 memberIndexSubEns_mod(:,subEnsIndex2)
354 memberIndex = memberIndex + nEnsPerSubEns_mod
355 end if
356 end do
357 end do
358
359 if ( mmpi_myid == 0 ) then
360 write(*,*) 'nEns, numSubEns, nEnsPerSubEns, nEnsIndependentPerSubEns = ', &
361 nEns, numSubEns, nEnsPerSubEns, nEnsIndependentPerSubEns
362 do subEnsIndex = 1, numSubEns
363 write(*,*) 'memberIndexSubEns = '
364 write(*,*) memberIndexSubEns(:,subEnsIndex)
365 if ( useModulatedEns ) then
366 write(*,*) 'memberIndexSubEns_mod = '
367 write(*,*) memberIndexSubEns_mod(:,subEnsIndex)
368 end if
369 write(*,*) 'memberIndexSubEnsComp = '
370 write(*,*) memberIndexSubEnsComp(:,subEnsIndex)
371 end do
372 end if
373
374 end if ! if CVLETKF(-PERTOBS)(-ME) algorithm
375
376 call lfn_Setup(LocFunctionWanted='FifthOrder')
377
378 ! compute 3D field of vertical location needed for localization
379 if (vLocalize > 0.0d0) then
380 call enkf_computeVertLocation(vertLocation_r4,stateVectorMeanTrl)
381 end if
382
383 allocate(varLevelFromK(numVarLev))
384 allocate(levFromK(numVarLev))
385 allocate(varKindFromK(numVarLev))
386 do varLevIndex = 1, numVarLev
387 varName = gsv_getVarNameFromK(stateVectorMeanInc,varLevIndex)
388 varLevelFromK(varLevIndex) = vnl_varLevelFromVarname(varName)
389 levFromK(varLevIndex) = gsv_getLevFromK(stateVectorMeanInc,varLevIndex)
390 varKindFromK(varLevIndex) = vnl_varKindFromVarname(varName)
391 end do
392
393 call utl_tmg_start(141,'----Barr')
394 call rpn_comm_barrier('GRID',ierr)
395 call utl_tmg_stop(141)
396
397 ! get mpi global list of tags used for mpi send/recv
398 call utl_tmg_start(142, '----GetGlobalTags')
399 allocate(latLonTagMpiGlobal(stateVectorMeanAnl%ni,stateVectorMeanAnl%nj))
400 call enkf_LETKFgetMpiGlobalTags(latLonTagMpiGlobal,myLatIndexesRecv,myLonIndexesRecv)
401 call utl_tmg_stop(142)
402
403 ! Compute the weights for ensemble mean and members
404 countMaxExceeded = 0
405 maxCountMaxExceeded = 0
406 numGridPointWeights = 0
407 LEV_LOOP: do levIndex = 1, nLev_weights
408 write(*,*) 'computing ensemble updates for vertical level = ', levIndex
409
410 !
411 ! First post all recv instructions for communication of weights
412 !
413 call utl_tmg_start(132,'----CommWeights')
414 numSend = 0
415 numRecv = 0
416 do latLonIndex = 1, myNumLatLonRecv
417 latIndex = myLatIndexesRecv(latLonIndex)
418 lonIndex = myLonIndexesRecv(latLonIndex)
419 procIndex = myProcIndexesRecv(latLonIndex)
420 recvTag = latLonTagMpiGlobal(lonIndex,latIndex)
421
422 nsize = nEnsGain
423 numRecv = numRecv + 1
424 call mpi_irecv( weightsMean(:,1,lonIndex,latIndex), &
425 nsize, mmpi_datyp_real8, procIndex-1, recvTag, &
426 mmpi_comm_grid, requestIdRecv(numRecv), ierr )
427 nsize = nEnsGain * nEns
428 numRecv = numRecv + 1
429 recvTag = recvTag + maxval(latLonTagMpiGlobal(:,:))
430 call mpi_irecv( weightsMembers(:,:,lonIndex,latIndex), &
431 nsize, mmpi_datyp_real8, procIndex-1, recvTag, &
432 mmpi_comm_grid, requestIdRecv(numRecv), ierr )
433 end do
434 call utl_tmg_stop(132)
435
436 LATLON_LOOP: do latLonIndex = 1, myNumLatLonSend
437 latIndex = myLatIndexesSend(latLonIndex)
438 lonIndex = myLonIndexesSend(latLonIndex)
439
440 numGridPointWeights = numGridPointWeights + 1
441
442 ! lat-lon of the grid point for which we are doing the analysis
443 anlLat = hco_ens%lat2d_4(lonIndex,latIndex)
444 anlLon = hco_ens%lon2d_4(lonIndex,latIndex)
445 hLocalizeIsConstant = all(hLocalize(:) == hLocalize(1))
446 if (vLocalize > 0.0d0 .or. .not.hLocalizeIsConstant) then
447 anlVertLocation = real(vertLocation_r4(lonIndex,latIndex,levIndex),8)
448 end if
449
450 ! Find which horizontal localization value to use for this analysis level
451 if (hLocalizeIsConstant) then
452 hLocIndex = 1
453 else
454 hLocIndex = 1 + count(anlVertLocation > hLocalizePressure(:))
455 end if
456
457 ! Get list of nearby observations and distances to gridpoint. With modulated-ensembles,
458 ! we get observations in entire column.
459 call utl_tmg_start(133,'----GetLocalBodyIndices')
460 if ( useModulatedEns ) anlVertLocation = MPC_missingValue_R8
461 numLocalObs = eob_getLocalBodyIndices(ensObs_mpiglobal, localBodyIndices, &
462 distances, anlLat, anlLon, anlVertLocation, &
463 hLocalize(hLocIndex), vLocalize, numLocalObsFound)
464 if (numLocalObsFound > maxNumLocalObs) then
465 countMaxExceeded = countMaxExceeded + 1
466 maxCountMaxExceeded = max(maxCountMaxExceeded, numLocalObsFound)
467 end if
468 call utl_tmg_stop(133)
469
470 call utl_tmg_start(134,'----CalculateWeights')
471
472 ! Extract initial quantities YbTinvR and first term of PaInv (YbTinvR*Yb)
473 do localObsIndex = 1, numLocalObs
474 bodyIndex = localBodyIndices(localObsIndex)
475
476 ! Compute value of localization function
477 ! Horizontal
478 localization = lfn_Response(distances(localObsIndex),hLocalize(hLocIndex))
479 ! Vertical when NOT using modulated ensembles - use pressures at the grid point (not obs) location
480 if (vLocalize > 0.0d0 .and. .not. useModulatedEns) then
481 distance = abs( anlVertLocation - ensObs_mpiglobal%vertLocation(bodyIndex) )
482 localization = localization * lfn_Response(distance,vLocalize)
483 end if
484 do memberIndex = 1, nEnsGain
485 ! YbTinvR for updating ensemble perturbations
486 YbTinvR_pert(memberIndex,localObsIndex) = &
487 ensObsGain_mpiglobal%Yb_r4(memberIndex, bodyIndex) * &
488 localization * ensObsGain_mpiglobal%obsErrInv(bodyIndex)
489 end do
490 if (eob_simObsAssim) then
491 do memberIndex = 1, nEnsGain
492 ! YbTinvR for the ensemble mean update for EDA observation simulation experiment
493 YbTinvR_mean(memberIndex,localObsIndex) = &
494 ensObsGain_mpiglobal%Yb_r4(memberIndex, bodyIndex) * &
495 localization * ensObsGain_mpiglobal%obsErrInv_sim(bodyIndex)
496 end do
497 end if
498 end do ! localObsIndex
499
500 call utl_tmg_start(136,'------CalcYbTinvRYb')
501
502 ! make copy of YbTinvR, and ensObsGain_mpiglobal%Yb_r4
503 call utl_tmg_start(137,'--------YbArraysCopy')
504 !$OMP PARALLEL DO PRIVATE (localObsIndex, bodyIndex, memberIndex2)
505 do localObsIndex = 1, numLocalObs
506 bodyIndex = localBodyIndices(localObsIndex)
507 do memberIndex2 = 1, nEnsGain
508 YbGainCopy_r4(localObsIndex,memberIndex2) = ensObsGain_mpiglobal%Yb_r4(memberIndex2,bodyIndex)
509 YbTinvRCopy_pert(localObsIndex,memberIndex2) = YbTinvR_pert(memberIndex2,localObsIndex)
510 end do
511 end do
512 !$OMP END PARALLEL DO
513 if (eob_simObsAssim) then
514 !$OMP PARALLEL DO PRIVATE (localObsIndex, bodyIndex, memberIndex2)
515 do localObsIndex = 1, numLocalObs
516 bodyIndex = localBodyIndices(localObsIndex)
517 do memberIndex2 = 1, nEnsGain
518 YbTinvRCopy_mean(localObsIndex,memberIndex2) = YbTinvR_mean(memberIndex2,localObsIndex)
519 end do
520 end do
521 !$OMP END PARALLEL DO
522 end if
523 call utl_tmg_stop(137)
524
525 call utl_tmg_start(138,'--------YbTinvRYb1')
526
527 YbTinvRYb_pert(:,:) = 0.0D0
528 !$OMP PARALLEL DO PRIVATE (memberIndex1, memberIndex2)
529 do memberIndex2 = 1, nEnsGain
530 do memberIndex1 = 1, memberIndex2 ! compute only upper triangle
531 YbTinvRYb_pert(memberIndex1,memberIndex2) = &
532 YbTinvRYb_pert(memberIndex1,memberIndex2) + &
533 sum(YbTinvRCopy_pert(1:numLocalObs,memberIndex1) * YbGainCopy_r4(1:numLocalObs,memberIndex2))
534 end do
535 end do
536 !$OMP END PARALLEL DO
537 ! copy upper triangle to lower triangle (symmetric matrix)
538 do memberIndex2 = 1, nEnsGain
539 do memberIndex1 = memberIndex2+1, nEnsGain
540 YbTinvRYb_pert(memberIndex1,memberIndex2) = &
541 YbTinvRYb_pert(memberIndex2,memberIndex1)
542 end do
543 end do
544 if (eob_simObsAssim) then
545 YbTinvRYb_mean(:,:) = 0.0D0
546 !$OMP PARALLEL DO PRIVATE (memberIndex1, memberIndex2)
547 do memberIndex2 = 1, nEnsGain
548 do memberIndex1 = 1, memberIndex2 ! compute only upper triangle
549 YbTinvRYb_mean(memberIndex1,memberIndex2) = &
550 YbTinvRYb_mean(memberIndex1,memberIndex2) + &
551 sum(YbTinvRCopy_mean(1:numLocalObs,memberIndex1) * YbGainCopy_r4(1:numLocalObs,memberIndex2))
552 end do
553 end do
554 !$OMP END PARALLEL DO
555 ! copy upper triangle to lower triangle (symmetric matrix)
556 do memberIndex2 = 1, nEnsGain
557 do memberIndex1 = memberIndex2+1, nEnsGain
558 YbTinvRYb_mean(memberIndex1,memberIndex2) = &
559 YbTinvRYb_mean(memberIndex2,memberIndex1)
560 end do
561 end do
562 end if
563 call utl_tmg_stop(138)
564
565 ! computing YbTinvRYb that uses modulated and original ensembles for perturbation update
566 if ( trim(algorithm) == 'CVLETKF-ME' .or. &
567 trim(algorithm) == 'LETKF-Gain-ME' ) then
568 ! make copy of ensObs_mpiglobal%Yb_r4
569 call utl_tmg_start(137,'--------YbArraysCopy')
570 YbCopy_r4(:,:) = 0.0
571 do localObsIndex = 1, numLocalObs
572 bodyIndex = localBodyIndices(localObsIndex)
573 do memberIndex2 = 1, nEns
574 YbCopy_r4(localObsIndex,memberIndex2) = ensObs_mpiglobal%Yb_r4(memberIndex2,bodyIndex)
575 end do
576 end do
577 call utl_tmg_stop(137)
578
579 YbTinvRYb_mod(:,:) = 0.0D0
580 call utl_tmg_start(139,'--------YbTinvRYb2')
581 !$OMP PARALLEL DO PRIVATE (memberIndex1, memberIndex2)
582 do memberIndex2 = 1, nEns
583 do memberIndex1 = 1, nEnsGain
584 YbTinvRYb_mod(memberIndex1,memberIndex2) = &
585 YbTinvRYb_mod(memberIndex1,memberIndex2) + &
586 sum(YbTinvRCopy_pert(1:numLocalObs,memberIndex1) * YbCopy_r4(1:numLocalObs,memberIndex2))
587 end do
588 end do
589 !$OMP END PARALLEL DO
590 call utl_tmg_stop(139)
591 end if !CVLETKF-ME or LETKF-GAIN-ME
592
593 call utl_tmg_stop(136)
594
595 ! Rest of the computation of local weights for this grid point
596 if (numLocalObs > 0) then
597
598 if (trim(algorithm) == 'LETKF') then
599 !
600 ! Weight calculation for standard LETKF algorithm
601 !
602
603 ! Add second term of PaInv
604 PaInv_pert(:,:) = YbTinvRYb_pert(:,:)
605 do memberIndex = 1, nEns
606 PaInv_pert(memberIndex,memberIndex) = PaInv_pert(memberIndex,memberIndex) + real(nEns - 1,8)
607 end do
608 if (eob_simObsAssim) then
609 PaInv_mean(:,:) = YbTinvRYb_mean(:,:)
610 do memberIndex = 1, nEns
611 PaInv_mean(memberIndex,memberIndex) = PaInv_mean(memberIndex,memberIndex) + real(nEns - 1,8)
612 end do
613 end if
614
615 ! Compute Pa and sqrt(Pa) matrices from PaInv
616 Pa_pert(:,:) = PaInv_pert(:,:)
617 call utl_tmg_start(135,'------EigenDecomp')
618 call utl_matInverse(Pa_pert, nEns, inverseSqrt_opt=PaSqrt_pert)
619 if (eob_simObsAssim) then
620 Pa_mean(:,:) = PaInv_mean(:,:)
621 call utl_matInverse(Pa_mean, nEns)
622 end if
623 call utl_tmg_stop(135)
624
625 ! Compute ensemble mean local weights as Pa * YbTinvR * (obs - meanYb)
626 weightsTemp(:) = 0.0d0
627 do localObsIndex = 1, numLocalObs
628 bodyIndex = localBodyIndices(localObsIndex)
629 do memberIndex = 1, nEns
630 weightsTemp(memberIndex) = weightsTemp(memberIndex) + &
631 YbTinvR_mean(memberIndex,localObsIndex) * &
632 ( ensObs_mpiglobal%obsValue(bodyIndex) - &
633 ensObs_mpiglobal%meanYb(bodyIndex) )
634 end do
635 end do
636
637 weightsMeanLatLon(:,1,latLonIndex) = 0.0d0
638 do memberIndex2 = 1, nEns
639 do memberIndex1 = 1, nEns
640 weightsMeanLatLon(memberIndex1,1,latLonIndex) = &
641 weightsMeanLatLon(memberIndex1,1,latLonIndex) + &
642 Pa_mean(memberIndex1,memberIndex2)*weightsTemp(memberIndex2)
643 end do
644 end do
645
646 ! Compute ensemble perturbation weights: [(Nens-1)^1/2*PaSqrt]
647 weightsMembersLatLon(:,:,latLonIndex) = sqrt(real(nEns - 1,8)) * PaSqrt_pert(:,:)
648
649 else if (trim(algorithm) == 'LETKF-Gain') then
650 !
651 ! Weight calculation for standard LETKF algorithm
652 !
653
654 ! Compute eigenValues/Vectors of Yb^T R^-1 Yb = E * Lambda * E^T
655 call utl_tmg_start(135,'------EigenDecomp')
656 tolerance = 1.0D-50
657 call utl_eigenDecomp(YbTinvRYb_pert, eigenValues_pert, eigenVectors_pert, tolerance, matrixRank)
658 if (eob_simObsAssim) then
659 call utl_eigenDecomp(YbTinvRYb_mean, eigenValues_mean, eigenVectors_mean, tolerance, matrixRank)
660 end if
661 call utl_tmg_stop(135)
662
663 ! Compute ensemble mean local weights as E * (Lambda + (Nens-1)*I)^-1 * E^T * YbTinvR * (obs - meanYb)
664 weightsTemp(:) = 0.0d0
665 do localObsIndex = 1, numLocalObs
666 bodyIndex = localBodyIndices(localObsIndex)
667 do memberIndex = 1, nEns
668 weightsTemp(memberIndex) = weightsTemp(memberIndex) + &
669 YbTinvR_mean(memberIndex,localObsIndex) * &
670 ( ensObs_mpiglobal%obsValue(bodyIndex) - &
671 ensObs_mpiglobal%meanYb(bodyIndex) )
672 end do
673 end do
674 weightsTemp2(:) = 0.0d0
675 do memberIndex2 = 1, matrixRank
676 do memberIndex1 = 1, nEns
677 weightsTemp2(memberIndex2) = weightsTemp2(memberIndex2) + &
678 eigenVectors_mean(memberIndex1,memberIndex2) * &
679 weightsTemp(memberIndex1)
680 end do
681 end do
682 do memberIndex = 1, matrixRank
683 weightsTemp2(memberIndex) = weightsTemp2(memberIndex) * &
684 1.0D0/(eigenValues_mean(memberIndex) + real(nEns - 1,8))
685 end do
686 weightsMeanLatLon(:,1,latLonIndex) = 0.0d0
687 do memberIndex2 = 1, matrixRank
688 do memberIndex1 = 1, nEns
689 weightsMeanLatLon(memberIndex1,1,latLonIndex) = &
690 weightsMeanLatLon(memberIndex1,1,latLonIndex) + &
691 eigenVectors_mean(memberIndex1,memberIndex2) * &
692 weightsTemp2(memberIndex2)
693 end do
694 end do
695
696 ! Compute ensemble perturbation weights:
697 ! Wa = [ - (Nens-1)^1/2 * E *
698 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} * Lambda^-1 *
699 ! E^T * YbTinvRYb ]
700 ! Loop over members within the current sub-ensemble being updated
701 do memberIndex = 1, nEns
702
703 ! E^T * YbTinvRYb
704 weightsTemp(:) = 0.0d0
705 do memberIndex2 = 1, matrixRank
706 do memberIndex1 = 1, nEns
707 weightsTemp(memberIndex2) = weightsTemp(memberIndex2) + &
708 eigenVectors_pert(memberIndex1,memberIndex2) * &
709 YbTinvRYb_pert(memberIndex1,memberIndex)
710 end do
711 end do
712
713 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} Lambda^-1 * previous_result
714
715 do memberIndex1 = 1, matrixRank
716 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) * &
717 ( 1.0D0/sqrt(real(nEns - 1,8)) - &
718 1.0D0/sqrt(eigenValues_pert(memberIndex1) + &
719 real(nEns - 1,8)) )
720 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) / &
721 eigenValues_pert(memberIndex1)
722 end do
723
724 ! E * previous_result
725 weightsMembersLatLon(:,memberIndex,latLonIndex) = 0.0d0
726 do memberIndex2 = 1, matrixRank
727 do memberIndex1 = 1, nEns
728 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) = &
729 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) + &
730 eigenVectors_pert(memberIndex1,memberIndex2) * &
731 weightsTemp(memberIndex2)
732 end do
733 end do
734
735 ! -1 * (Nens-1)^1/2 * previous_result
736 weightsMembersLatLon(:,memberIndex,latLonIndex) = &
737 -1.0D0 * sqrt(real(nEns - 1,8)) * &
738 weightsMembersLatLon(:,memberIndex,latLonIndex)
739
740 ! I + previous_result
741 weightsMembersLatLon(memberIndex,memberIndex,latLonIndex) = &
742 1.0D0 + weightsMembersLatLon(memberIndex,memberIndex,latLonIndex)
743
744 end do
745
746 ! Remove the weights mean computed over the columns
747 do memberIndex = 1, nEns
748 weightsMembersLatLon(memberIndex,:,latLonIndex) = &
749 weightsMembersLatLon(memberIndex,:,latLonIndex) - &
750 sum(weightsMembersLatLon(memberIndex,:,latLonIndex))/real(nEns,8)
751 end do
752
753 else if (trim(algorithm) == 'LETKF-Gain-ME') then
754 !
755 ! Weight calculation for standard LETKF algorithm with modulated ensemble
756 !
757
758 ! Compute eigenValues/Vectors of Yb^T R^-1 Yb = E * Lambda * E^T
759 call utl_tmg_start(135,'------EigenDecomp')
760 tolerance = 1.0D-50
761 call utl_eigenDecomp(YbTinvRYb_pert, eigenValues_pert, eigenVectors_pert, tolerance, matrixRank)
762 if (eob_simObsAssim) then
763 call utl_eigenDecomp(YbTinvRYb_mean, eigenValues_mean, eigenVectors_mean, tolerance, matrixRank)
764 end if
765 call utl_tmg_stop(135)
766
767 ! Compute ensemble mean local weights as E * (Lambda + (Nens-1)*I)^-1 * E^T * YbTinvR * (obs - meanYb)
768 weightsTemp(:) = 0.0d0
769 do localObsIndex = 1, numLocalObs
770 bodyIndex = localBodyIndices(localObsIndex)
771 do memberIndex = 1, nEnsGain
772 weightsTemp(memberIndex) = weightsTemp(memberIndex) + &
773 YbTinvR_mean(memberIndex,localObsIndex) * &
774 (ensObs_mpiglobal%obsValue(bodyIndex) - &
775 ensObs_mpiglobal%meanYb(bodyIndex))
776 end do
777 end do
778 weightsTemp2(:) = 0.0d0
779 do memberIndex2 = 1, matrixRank
780 do memberIndex1 = 1, nEnsGain
781 weightsTemp2(memberIndex2) = weightsTemp2(memberIndex2) + &
782 eigenVectors_mean(memberIndex1,memberIndex2) * &
783 weightsTemp(memberIndex1)
784 end do
785 end do
786 do memberIndex = 1, matrixRank
787 weightsTemp2(memberIndex) = weightsTemp2(memberIndex) * &
788 1.0D0/(eigenValues_mean(memberIndex) + real(nEnsGain - 1,8))
789 end do
790 weightsMeanLatLon(:,1,latLonIndex) = 0.0d0
791 do memberIndex2 = 1, matrixRank
792 do memberIndex1 = 1, nEnsGain
793 weightsMeanLatLon(memberIndex1,1,latLonIndex) = &
794 weightsMeanLatLon(memberIndex1,1,latLonIndex) + &
795 eigenVectors_mean(memberIndex1,memberIndex2) * &
796 weightsTemp2(memberIndex2)
797 end do
798 end do
799
800 ! Compute ensemble perturbation weights:
801 ! Wa = [ - (Nens-1)^1/2 * E *
802 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} * Lambda^-1 *
803 ! E^T * YbTinvRYb_mod ]
804 ! Loop over members within the current sub-ensemble being updated
805 do memberIndex = 1, nEns
806
807 ! E^T * YbTinvRYb_mod
808 weightsTemp(:) = 0.0d0
809 do memberIndex2 = 1, matrixRank
810 do memberIndex1 = 1, nEnsGain
811 weightsTemp(memberIndex2) = weightsTemp(memberIndex2) + &
812 eigenVectors_pert(memberIndex1,memberIndex2) * &
813 YbTinvRYb_mod(memberIndex1,memberIndex)
814 end do
815 end do
816
817 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} Lambda^-1 * previous_result
818
819 do memberIndex1 = 1, matrixRank
820 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) * &
821 ( 1.0D0/sqrt(real(nEnsGain - 1,8)) - &
822 1.0D0/sqrt(eigenValues_pert(memberIndex1) + &
823 real(nEnsGain - 1,8)) )
824 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) / &
825 eigenValues_pert(memberIndex1)
826 end do
827
828 ! E * previous_result
829 weightsMembersLatLon(:,memberIndex,latLonIndex) = 0.0d0
830 do memberIndex2 = 1, matrixRank
831 do memberIndex1 = 1, nEnsGain
832 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) = &
833 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) + &
834 eigenVectors_pert(memberIndex1,memberIndex2) * &
835 weightsTemp(memberIndex2)
836 end do
837 end do
838
839 ! -1 * (Nens-1)^1/2 * previous_result
840 weightsMembersLatLon(:,memberIndex,latLonIndex) = &
841 -1.0D0 * sqrt(real(nEnsGain - 1,8)) * &
842 weightsMembersLatLon(:,memberIndex,latLonIndex)
843
844 end do
845
846 ! Remove the weights mean computed over the columns
847 do memberIndex = 1, nEnsGain
848 weightsMembersLatLon(memberIndex,:,latLonIndex) = &
849 weightsMembersLatLon(memberIndex,:,latLonIndex) - &
850 sum(weightsMembersLatLon(memberIndex,:,latLonIndex))/real(nEns,8)
851 end do
852
853 else if (trim(algorithm) == 'CVLETKF') then
854 !
855 ! Weight calculation for cross-validation LETKF algorithm
856 !
857
858 ! Compute eigenValues/Vectors of Yb^T R^-1 Yb = E * Lambda * E^T
859 call utl_tmg_start(135,'------EigenDecomp')
860 tolerance = 1.0D-50
861 call utl_eigenDecomp(YbTinvRYb_mean, eigenValues_mean, eigenVectors_mean, tolerance, matrixRank)
862 call utl_tmg_stop(135)
863
864 ! Compute ensemble mean local weights as E * (Lambda + (Nens-1)*I)^-1 * E^T * YbTinvR * (obs - meanYb)
865 weightsTemp(:) = 0.0d0
866 do localObsIndex = 1, numLocalObs
867 bodyIndex = localBodyIndices(localObsIndex)
868 do memberIndex = 1, nEns
869 weightsTemp(memberIndex) = weightsTemp(memberIndex) + &
870 YbTinvR_mean(memberIndex,localObsIndex) * &
871 ( ensObs_mpiglobal%obsValue(bodyIndex) - &
872 ensObs_mpiglobal%meanYb(bodyIndex) )
873 end do
874 end do
875 weightsTemp2(:) = 0.0d0
876 do memberIndex2 = 1, matrixRank
877 do memberIndex1 = 1, nEns
878 weightsTemp2(memberIndex2) = weightsTemp2(memberIndex2) + &
879 eigenVectors_mean(memberIndex1,memberIndex2) * &
880 weightsTemp(memberIndex1)
881 end do
882 end do
883 do memberIndex = 1, matrixRank
884 weightsTemp2(memberIndex) = weightsTemp2(memberIndex) * &
885 1.0D0/(eigenValues_mean(memberIndex) + real(nEns - 1,8))
886 end do
887 weightsMeanLatLon(:,1,latLonIndex) = 0.0d0
888 do memberIndex2 = 1, matrixRank
889 do memberIndex1 = 1, nEns
890 weightsMeanLatLon(memberIndex1,1,latLonIndex) = &
891 weightsMeanLatLon(memberIndex1,1,latLonIndex) + &
892 eigenVectors_mean(memberIndex1,memberIndex2) * &
893 weightsTemp2(memberIndex2)
894 end do
895 end do
896
897 ! Compute ensemble perturbation weights:
898 ! Wa = [ I - (Nens-1)^1/2 * E *
899 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} * Lambda^-1 *
900 ! E^T * YbTinvRYb ]
901 ! Loop over sub-ensembles
902
903 !$OMP PARALLEL DO PRIVATE(subEnsIndex, memberIndexCV, memberIndexCV1, memberIndexCV2, &
904 !$OMP memberIndex, memberIndex1, memberIndex2, weightsTemp, tolerance, &
905 !$OMP YbTinvRYb_CV, eigenValues_CV, eigenVectors_CV, matrixRank)
906 do subEnsIndex = 1, numSubEns
907
908 ! Use complement (independent) ens to get eigenValues/Vectors of Yb^T R^-1 Yb = E*Lambda*E^T
909 call utl_tmg_start(135,'------EigenDecomp')
910 do memberIndexCV2 = 1, nEnsIndependentPerSubEns
911 memberIndex2 = memberIndexSubEnsComp(memberIndexCV2, subEnsIndex)
912 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
913 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
914 YbTinvRYb_CV(memberIndexCV1,memberIndexCV2) = YbTinvRYb_pert(memberIndex1,memberIndex2)
915 end do
916 end do
917 tolerance = 1.0D-50
918 call utl_eigenDecomp(YbTinvRYb_CV, eigenValues_CV, eigenVectors_CV, tolerance, matrixRank)
919 call utl_tmg_stop(135)
920
921 ! Loop over members within the current sub-ensemble being updated
922 do memberIndexCV = 1, nEnsPerSubEns
923
924 ! This is index of member being updated
925 memberIndex = memberIndexSubEns(memberIndexCV, subEnsIndex)
926
927 ! E^T * YbTinvRYb
928 weightsTemp(:) = 0.0d0
929 do memberIndex2 = 1, matrixRank
930 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
931 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
932 weightsTemp(memberIndex2) = weightsTemp(memberIndex2) + &
933 eigenVectors_CV(memberIndexCV1,memberIndex2) * &
934 YbTinvRYb_pert(memberIndex1,memberIndex)
935 end do
936 end do
937
938 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} Lambda^-1 * previous_result
939
940 do memberIndex1 = 1, matrixRank
941 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) * &
942 ( 1.0D0/sqrt(real(nEnsIndependentPerSubEns - 1,8)) - &
943 1.0D0/sqrt(eigenValues_CV(memberIndex1) + &
944 real(nEnsIndependentPerSubEns - 1,8)) )
945 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) / &
946 eigenValues_CV(memberIndex1)
947 end do
948
949 ! E * previous_result
950 weightsMembersLatLon(:,memberIndex,latLonIndex) = 0.0d0
951 do memberIndex2 = 1, matrixRank
952 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
953 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
954 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) = &
955 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) + &
956 eigenVectors_CV(memberIndexCV1,memberIndex2) * &
957 weightsTemp(memberIndex2)
958 end do
959 end do
960
961 ! -1 * (Nens-1)^1/2 * previous_result
962 weightsMembersLatLon(:,memberIndex,latLonIndex) = &
963 -1.0D0 * sqrt(real(nEnsIndependentPerSubEns - 1,8)) * &
964 weightsMembersLatLon(:,memberIndex,latLonIndex)
965
966 ! I + previous_result
967 weightsMembersLatLon(memberIndex,memberIndex,latLonIndex) = &
968 1.0D0 + weightsMembersLatLon(memberIndex,memberIndex,latLonIndex)
969
970 end do ! memberIndexCV
971 end do ! subEnsIndex
972 !$OMP END PARALLEL DO
973
974 ! Remove the weights mean computed over the columns
975 do memberIndex = 1, nEns
976 weightsMembersLatLon(memberIndex,:,latLonIndex) = &
977 weightsMembersLatLon(memberIndex,:,latLonIndex) - &
978 sum(weightsMembersLatLon(memberIndex,:,latLonIndex))/real(nEns,8)
979 end do
980
981 else if (trim(algorithm) == 'CVLETKF-ME') then
982 !
983 ! Weight calculation for cross-validation LETKF algorithm
984 !
985
986 ! Compute eigenValues/Vectors of Yb^T R^-1 Yb = E * Lambda * E^T
987 call utl_tmg_start(135,'------EigenDecomp')
988 tolerance = 1.0D-50
989 call utl_eigenDecomp(YbTinvRYb_mean, eigenValues_mean, eigenVectors_mean, tolerance, matrixRank)
990 call utl_tmg_stop(135)
991 !if (matrixRank < (nEns-1)) then
992 ! write(*,*) 'YbTinvRYb is rank deficient =', matrixRank, nEns, numLocalObs
993 !end if
994
995 ! Compute ensemble mean local weights as E * (Lambda + (Nens-1)*I)^-1 * E^T * YbTinvR * (obs - meanYb)
996 weightsTemp(:) = 0.0d0
997 do localObsIndex = 1, numLocalObs
998 bodyIndex = localBodyIndices(localObsIndex)
999 do memberIndex = 1, nEnsGain
1000 weightsTemp(memberIndex) = weightsTemp(memberIndex) + &
1001 YbTinvR_mean(memberIndex,localObsIndex) * &
1002 ( ensObs_mpiglobal%obsValue(bodyIndex) - &
1003 ensObs_mpiglobal%meanYb(bodyIndex) )
1004 end do
1005 end do
1006 weightsTemp2(:) = 0.0d0
1007 do memberIndex2 = 1, matrixRank
1008 do memberIndex1 = 1, nEnsGain
1009 weightsTemp2(memberIndex2) = weightsTemp2(memberIndex2) + &
1010 eigenVectors_mean(memberIndex1,memberIndex2) * &
1011 weightsTemp(memberIndex1)
1012 end do
1013 end do
1014 do memberIndex = 1, matrixRank
1015 weightsTemp2(memberIndex) = weightsTemp2(memberIndex) * &
1016 1.0D0/(eigenValues_mean(memberIndex) + real(nEnsGain - 1,8))
1017 end do
1018 weightsMeanLatLon(:,1,latLonIndex) = 0.0d0
1019 do memberIndex2 = 1, matrixRank
1020 do memberIndex1 = 1, nEnsGain
1021 weightsMeanLatLon(memberIndex1,1,latLonIndex) = &
1022 weightsMeanLatLon(memberIndex1,1,latLonIndex) + &
1023 eigenVectors_mean(memberIndex1,memberIndex2) * &
1024 weightsTemp2(memberIndex2)
1025 end do
1026 end do
1027
1028 ! Compute ensemble perturbation weights:
1029 ! Wa = [ - (Nens-1)^1/2 * E *
1030 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} * Lambda^-1 *
1031 ! E^T * YbTinvRYb_mod ]
1032 ! Loop over sub-ensembles
1033 !$OMP PARALLEL DO PRIVATE(subEnsIndex, memberIndexCV, memberIndexCV1, memberIndexCV2, &
1034 !$OMP memberIndex, memberIndex1, memberIndex2, weightsTemp, tolerance, &
1035 !$OMP YbTinvRYb_CV, eigenValues_CV, eigenVectors_CV, matrixRank)
1036 do subEnsIndex = 1, numSubEns
1037
1038 ! Use complement (independent) ens to get eigenValues/Vectors of Yb^T R^-1 Yb = E*Lambda*E^T
1039 call utl_tmg_start(135,'------EigenDecomp')
1040 do memberIndexCV2 = 1, nEnsIndependentPerSubEns
1041 memberIndex2 = memberIndexSubEnsComp(memberIndexCV2, subEnsIndex)
1042 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
1043 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
1044 YbTinvRYb_CV(memberIndexCV1,memberIndexCV2) = YbTinvRYb_pert(memberIndex1,memberIndex2)
1045 end do
1046 end do
1047 tolerance = 1.0D-50
1048 call utl_eigenDecomp(YbTinvRYb_CV, eigenValues_CV, eigenVectors_CV, tolerance, matrixRank)
1049 call utl_tmg_stop(135)
1050
1051 ! Loop over members within the current sub-ensemble being updated
1052 do memberIndexCV = 1, nEnsPerSubEns
1053
1054 ! This is index of member being updated
1055 memberIndex = memberIndexSubEns(memberIndexCV, subEnsIndex)
1056
1057 ! E^T * YbTinvRYb
1058 weightsTemp(:) = 0.0d0
1059 do memberIndex2 = 1, matrixRank
1060 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
1061 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
1062 weightsTemp(memberIndex2) = weightsTemp(memberIndex2) + &
1063 eigenVectors_CV(memberIndexCV1,memberIndex2) * &
1064 YbTinvRYb_mod(memberIndex1,memberIndex)
1065 end do
1066 end do
1067
1068 ! {(Nens-1)^-1/2*I - (Lambda + (Nens-1)*I)^-1/2} Lambda^-1 * previous_result
1069
1070 do memberIndex1 = 1, matrixRank
1071 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) * &
1072 ( 1.0D0/sqrt(real(nEnsIndependentPerSubEns - 1,8)) - &
1073 1.0D0/sqrt(eigenValues_CV(memberIndex1) + &
1074 real(nEnsIndependentPerSubEns - 1,8)) )
1075 weightsTemp(memberIndex1) = weightsTemp(memberIndex1) / &
1076 eigenValues_CV(memberIndex1)
1077 end do
1078
1079 ! E * previous_result
1080 weightsMembersLatLon(:,memberIndex,latLonIndex) = 0.0d0
1081 do memberIndex2 = 1, matrixRank
1082 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
1083 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
1084 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) = &
1085 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) + &
1086 eigenVectors_CV(memberIndexCV1,memberIndex2) * &
1087 weightsTemp(memberIndex2)
1088 end do
1089 end do
1090
1091 ! -1 * (Nens-1)^1/2 * previous_result
1092 weightsMembersLatLon(:,memberIndex,latLonIndex) = &
1093 -1.0D0 * sqrt(real(nEnsIndependentPerSubEns - 1,8)) * &
1094 weightsMembersLatLon(:,memberIndex,latLonIndex)
1095
1096 end do ! memberIndexCV
1097 end do ! subEnsIndex
1098 !$OMP END PARALLEL DO
1099
1100 ! Remove the weights mean computed over the columns
1101 do memberIndex = 1, nEnsGain
1102 weightsMembersLatLon(memberIndex,:,latLonIndex) = &
1103 weightsMembersLatLon(memberIndex,:,latLonIndex) - &
1104 sum(weightsMembersLatLon(memberIndex,:,latLonIndex))/real(nEns,8)
1105 end do
1106
1107 else if (trim(algorithm) == 'CVLETKF-PERTOBS') then
1108 !
1109 ! Weight calculation for perturbed-obs cross-validation LETKF algorithm
1110 !
1111
1112 ! Compute eigenValues/Vectors of Yb^T R^-1 Yb = E * Lambda * E^T
1113 call utl_tmg_start(135,'------EigenDecomp')
1114 tolerance = 1.0D-50
1115 call utl_eigenDecomp(YbTinvRYb_mean, eigenValues_mean, eigenVectors_mean, tolerance, matrixRank)
1116 call utl_tmg_stop(135)
1117
1118 ! Compute ensemble mean local weights as E * (Lambda + (Nens-1)*I)^-1 * E^T * YbTinvR * (obs - meanYb)
1119 weightsTemp(:) = 0.0d0
1120 do localObsIndex = 1, numLocalObs
1121 bodyIndex = localBodyIndices(localObsIndex)
1122 do memberIndex = 1, nEns
1123 weightsTemp(memberIndex) = weightsTemp(memberIndex) + &
1124 YbTinvR_mean(memberIndex,localObsIndex) * &
1125 ( ensObs_mpiglobal%obsValue(bodyIndex) - &
1126 ensObs_mpiglobal%meanYb(bodyIndex) )
1127 end do
1128 end do
1129 weightsTemp2(:) = 0.0d0
1130 do memberIndex2 = 1, matrixRank
1131 do memberIndex1 = 1, nEns
1132 weightsTemp2(memberIndex2) = weightsTemp2(memberIndex2) + &
1133 eigenVectors_mean(memberIndex1,memberIndex2) * &
1134 weightsTemp(memberIndex1)
1135 end do
1136 end do
1137 do memberIndex = 1, matrixRank
1138 weightsTemp2(memberIndex) = weightsTemp2(memberIndex) * &
1139 1.0D0/(eigenValues_mean(memberIndex) + real(nEns - 1,8))
1140 end do
1141 weightsMeanLatLon(:,1,latLonIndex) = 0.0d0
1142 do memberIndex2 = 1, matrixRank
1143 do memberIndex1 = 1, nEns
1144 weightsMeanLatLon(memberIndex1,1,latLonIndex) = &
1145 weightsMeanLatLon(memberIndex1,1,latLonIndex) + &
1146 eigenVectors_mean(memberIndex1,memberIndex2) * &
1147 weightsTemp2(memberIndex2)
1148 end do
1149 end do
1150
1151 ! Compute ensemble perturbation weights using mean increment weights
1152 ! formula, but with subset of members:
1153 ! wa_i = I_i + E * (Lambda + (Nens-1)*I)^-1 * E^T * YbTinvR * (obs + randpert_i - Yb_i)
1154 ! Wa = wa_i - mean_over_i(wa_i)
1155 !
1156 ! Loop over sub-ensembles
1157 !$OMP PARALLEL DO PRIVATE(subEnsIndex, memberIndexCV, memberIndexCV1, memberIndexCV2, &
1158 !$OMP memberIndex, memberIndex1, memberIndex2, weightsTemp, tolerance, &
1159 !$OMP YbTinvRYb_CV, eigenValues_CV, eigenVectors_CV, matrixRank, &
1160 !$OMP weightsTemp2, localObsIndex, bodyIndex)
1161 do subEnsIndex = 1, numSubEns
1162
1163 ! Use complement (independent) ens to get eigenValues/Vectors of Yb^T R^-1 Yb = E*Lambda*E^T
1164 call utl_tmg_start(135,'------EigenDecomp')
1165 do memberIndexCV2 = 1, nEnsIndependentPerSubEns
1166 memberIndex2 = memberIndexSubEnsComp(memberIndexCV2, subEnsIndex)
1167 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
1168 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
1169 YbTinvRYb_CV(memberIndexCV1,memberIndexCV2) = YbTinvRYb_pert(memberIndex1,memberIndex2)
1170 end do
1171 end do
1172 tolerance = 1.0D-50
1173 call utl_eigenDecomp(YbTinvRYb_CV, eigenValues_CV, eigenVectors_CV, tolerance, matrixRank)
1174 call utl_tmg_stop(135)
1175
1176 ! Loop over members within the current sub-ensemble being updated
1177 do memberIndexCV = 1, nEnsPerSubEns
1178
1179 ! This is index of member being updated (i'th member)
1180 memberIndex = memberIndexSubEns(memberIndexCV, subEnsIndex)
1181
1182 ! YbTinvRYb * (obsValue + randPert_i - Yb_i)
1183 weightsTemp(:) = 0.0d0
1184 do localObsIndex = 1, numLocalObs
1185 bodyIndex = localBodyIndices(localObsIndex)
1186 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
1187 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
1188 weightsTemp(memberIndexCV1) = &
1189 weightsTemp(memberIndexCV1) + &
1190 YbTinvR_pert(memberIndex1,localObsIndex) * &
1191 ( ensObs_mpiglobal%obsValue(bodyIndex) + &
1192 ensObs_mpiglobal%randPert_r4(memberIndex,bodyIndex) - &
1193 ( ensObs_mpiglobal%meanYb(bodyIndex) + &
1194 ensObs_mpiglobal%Yb_r4(memberIndex,bodyIndex) ) )
1195 end do
1196 end do
1197
1198 ! E^T * previous_result
1199 weightsTemp2(:) = 0.0d0
1200 do memberIndex2 = 1, matrixRank
1201 do memberIndex1 = 1, nEnsIndependentPerSubEns
1202 weightsTemp2(memberIndex2) = weightsTemp2(memberIndex2) + &
1203 eigenVectors_CV(memberIndex1,memberIndex2) * &
1204 weightsTemp(memberIndex1)
1205 end do
1206 end do
1207
1208 ! [lambda + (N_indep-1)*I]^-1 * previous_result
1209 do memberIndex1 = 1, matrixRank
1210 weightsTemp2(memberIndex1) = &
1211 weightsTemp2(memberIndex1) * &
1212 1.0D0/(eigenValues_CV(memberIndex1) + real(nEnsIndependentPerSubEns - 1,8))
1213 end do
1214
1215 ! E * previous_result
1216 weightsMembersLatLon(:,memberIndex,latLonIndex) = 0.0d0
1217 do memberIndex2 = 1, matrixRank
1218 do memberIndexCV1 = 1, nEnsIndependentPerSubEns
1219 memberIndex1 = memberIndexSubEnsComp(memberIndexCV1, subEnsIndex)
1220 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) = &
1221 weightsMembersLatLon(memberIndex1,memberIndex,latLonIndex) + &
1222 eigenVectors_CV(memberIndexCV1,memberIndex2) * &
1223 weightsTemp2(memberIndex2)
1224 end do
1225 end do
1226
1227 ! I + previous_result
1228 weightsMembersLatLon(memberIndex,memberIndex,latLonIndex) = &
1229 1.0D0 + weightsMembersLatLon(memberIndex,memberIndex,latLonIndex)
1230
1231 end do ! memberIndexCV
1232 end do ! subEnsIndex
1233 !$OMP END PARALLEL DO
1234
1235 ! Remove the weights mean computed over the columns
1236 do memberIndex = 1, nEns
1237 weightsMembersLatLon(memberIndex,:,latLonIndex) = &
1238 weightsMembersLatLon(memberIndex,:,latLonIndex) - &
1239 sum(weightsMembersLatLon(memberIndex,:,latLonIndex))/real(nEns,8)
1240 end do
1241
1242 else
1243
1244 call utl_abort('UNKNOWN LETKF ALGORITHM')
1245
1246 end if
1247
1248 else
1249
1250 ! no obs near this grid point, mean weights zero, member weights identity
1251 weightsMeanLatLon(:,1,latLonIndex) = 0.0d0
1252 weightsMembersLatLon(:,:,latLonIndex) = 0.0d0
1253 do memberIndex = 1, nEns
1254 if ( useModulatedEns ) then
1255 do eigenVectorColumnIndex = 1, numRetainedEigen
1256 memberIndexInModEns = (eigenVectorColumnIndex - 1) * nEns + memberIndex
1257 weightsMembersLatLon(memberIndexInModEns,memberIndex,latLonIndex) = 1.0d0
1258 end do
1259 else
1260 weightsMembersLatLon(memberIndex,memberIndex,latLonIndex) = 1.0d0
1261 end if
1262 end do
1263
1264 end if ! numLocalObs > 0
1265
1266 call utl_tmg_stop(134)
1267
1268 !
1269 ! Now post all send instructions (each lat-lon may be sent to multiple tasks)
1270 !
1271 call utl_tmg_start(132,'----CommWeights')
1272 latIndex = myLatIndexesSend(latLonIndex)
1273 lonIndex = myLonIndexesSend(latLonIndex)
1274 do procIndex = 1, myNumProcIndexesSend(latLonIndex)
1275 sendTag = latLonTagMpiGlobal(lonIndex,latIndex)
1276 procIndexSend = myProcIndexesSend(latLonIndex, procIndex)
1277
1278 nsize = nEnsGain
1279 numSend = numSend + 1
1280 call mpi_isend( weightsMeanLatLon(:,1,latLonIndex), &
1281 nsize, mmpi_datyp_real8, procIndexSend-1, sendTag, &
1282 mmpi_comm_grid, requestIdSend(numSend), ierr )
1283 nsize = nEnsGain * nEns
1284 numSend = numSend + 1
1285 sendTag = sendTag + maxval(latLonTagMpiGlobal(:,:))
1286 call mpi_isend( weightsMembersLatLon(:,:,latLonIndex), &
1287 nsize, mmpi_datyp_real8, procIndexSend-1, sendTag, &
1288 mmpi_comm_grid, requestIdSend(numSend), ierr )
1289 end do
1290 call utl_tmg_stop(132)
1291
1292 end do LATLON_LOOP
1293
1294 !
1295 ! Wait for communiations to finish before continuing
1296 !
1297 call utl_tmg_start(132,'----CommWeights')
1298 if (firstTime) write(*,*) 'numSend/Recv = ', numSend, numRecv
1299 firstTime = .false.
1300
1301 if ( numRecv > 0 ) then
1302 call mpi_waitAll(numRecv, requestIdRecv(1:numRecv), MPI_STATUSES_IGNORE, ierr)
1303 end if
1304
1305 if ( numSend > 0 ) then
1306 call mpi_waitAll(numSend, requestIdSend(1:numSend), MPI_STATUSES_IGNORE, ierr)
1307 end if
1308
1309 call utl_tmg_stop(132)
1310
1311 !
1312 ! Interpolate weights from coarse to full resolution
1313 !
1314 call utl_tmg_start(140,'----InterpolateWeights')
1315 if (wInterpInfo%latLonStep > 1) then
1316 call enkf_interpWeights(wInterpInfo, weightsMean)
1317 call enkf_interpWeights(wInterpInfo, weightsMembers)
1318 end if
1319 call utl_tmg_stop(140)
1320
1321 call utl_tmg_start(143,'----ApplyWeights')
1322
1323 !
1324 ! Apply the weights to compute the ensemble mean and members
1325 !
1326 call gsv_getField(stateVectorMeanInc,meanInc_ptr_r4)
1327 call gsv_getField(stateVectorMeanTrl,meanTrl_ptr_r4)
1328 call gsv_getField(stateVectorMeanAnl,meanAnl_ptr_r4)
1329
1330 !$OMP PARALLEL DO PRIVATE(latIndex, lonIndex, varLevIndex, levIndex2, memberTrl_ptr_r4, memberAnl_ptr_r4), &
1331 !$OMP PRIVATE(memberAnlPert, stepIndex, memberIndex, memberIndex2, memberIndex1, eigenVectorColumnIndex, pert_r4), &
1332 !$OMP PRIVATE(memberIndexInModEns, modulationFactor_r4)
1333 do latIndex = myLatBeg, myLatEnd
1334 LON_LOOP5: do lonIndex = myLonBeg, myLonEnd
1335
1336 ! skip this grid point if all weights zero (no nearby obs)
1337 if (all(weightsMean(:,1,lonIndex,latIndex) == 0.0d0)) cycle LON_LOOP5
1338
1339 ! Compute the ensemble mean increment and analysis
1340 do varLevIndex = 1, numVarLev
1341 ! Only treat varLevIndex values that correspond with current levIndex
1342 if (varLevelFromK(varLevIndex) == 'SF' .or. varLevelFromK(varLevIndex) == 'SFMM' .or. &
1343 varLevelFromK(varLevIndex) == 'SFTH' .or. varLevelFromK(varLevIndex) == 'SS') then
1344 if (varKindFromK(varLevIndex) == 'OC') then
1345 levIndex2 = 1
1346 else
1347 levIndex2 = max(nLev_M,nLev_depth)
1348 end if
1349 else if (varLevelFromK(varLevIndex) == 'MM' .or. varLevelFromK(varLevIndex) == 'TH' .or. varLevelFromK(varLevIndex) == 'DP') then
1350 levIndex2 = levFromK(varLevIndex)
1351 else if (varLevelFromK(varLevIndex) == 'OT') then
1352 ! Most (all?) variables using the 'other' coordinate are surface
1353 levIndex2 = max(nLev_M,nLev_depth)
1354 else
1355 write(*,*) 'varLevel = ', varLevelFromK(varLevIndex)
1356 call utl_abort('enkf_LETKFanalyses: unknown varLevel')
1357 end if
1358 if (levIndex2 /= levIndex .and. .not. useModulatedEns) cycle
1359 memberTrl_ptr_r4 => ens_getOneLev_r4(ensembleTrl,varLevIndex)
1360 do stepIndex = 1, tim_nstepobsinc
1361 ! mean increment
1362 if ( useModulatedEns ) then
1363 do eigenVectorColumnIndex = 1, numRetainedEigen
1364 call getModulationFactor( stateVectorMeanInc%vco, levIndex2, &
1365 eigenVectorColumnIndex, numRetainedEigen, &
1366 nEns, vLocalize, &
1367 modulationFactor_r4 )
1368
1369 do memberIndex = 1, nEns
1370 pert_r4 = modulationFactor_r4 * ( memberTrl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) - &
1371 meanTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) )
1372
1373 ! Index of the modulated ensemble member corresponding to original
1374 ! ensemble member index (memberIndex1) and eigenVectorColumnIndex.
1375 memberIndexInModEns = (eigenVectorColumnIndex - 1) * nEns + memberIndex
1376
1377 meanInc_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) = &
1378 meanInc_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) + &
1379 weightsMean(memberIndexInModEns,1,lonIndex,latIndex) * pert_r4
1380 end do
1381 end do
1382 else
1383 do memberIndex = 1, nEns
1384 meanInc_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) = &
1385 meanInc_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) + &
1386 weightsMean(memberIndex,1,lonIndex,latIndex) * &
1387 (memberTrl_ptr_r4(memberIndex,stepIndex,lonIndex,latIndex) - &
1388 meanTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex))
1389 end do
1390 end if
1391
1392 ! mean analysis
1393 meanAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) = &
1394 meanTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) + &
1395 meanInc_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex)
1396 end do ! stepIndex
1397 end do ! varLevIndex
1398
1399 ! Compute the ensemble member analyses
1400 do varLevIndex = 1, numVarLev
1401 ! Only treat varLevIndex values that correspond with current levIndex
1402 if (varLevelFromK(varLevIndex) == 'SF' .or. varLevelFromK(varLevIndex) == 'SFMM' .or. &
1403 varLevelFromK(varLevIndex) == 'SFTH' .or. varLevelFromK(varLevIndex) == 'SS') then
1404 if (varKindFromK(varLevIndex) == 'OC') then
1405 levIndex2 = 1
1406 else
1407 levIndex2 = max(nLev_M,nLev_depth)
1408 end if
1409 else if (varLevelFromK(varLevIndex) == 'MM' .or. varLevelFromK(varLevIndex) == 'TH' .or. varLevelFromK(varLevIndex) == 'DP') then
1410 levIndex2 = levFromK(varLevIndex)
1411 else if (varLevelFromK(varLevIndex) == 'OT') then
1412 ! Most (all?) variables using the 'other' coordinate are surface
1413 levIndex2 = max(nLev_M,nLev_depth)
1414 else
1415 write(*,*) 'varLevel = ', varLevelFromK(varLevIndex)
1416 call utl_abort('enkf_LETKFanalyses: unknown varLevel')
1417 end if
1418 if (levIndex2 /= levIndex .and. .not. useModulatedEns) cycle
1419 memberTrl_ptr_r4 => ens_getOneLev_r4(ensembleTrl,varLevIndex)
1420 memberAnl_ptr_r4 => ens_getOneLev_r4(ensembleAnl,varLevIndex)
1421 do stepIndex = 1, tim_nstepobsinc
1422
1423 call utl_tmg_start(144,'------ApplyWeightsMember')
1424
1425 ! Compute analysis member perturbation
1426 memberAnlPert(:) = 0.0d0
1427
1428 if ( useModulatedEns ) then
1429 do memberIndex2 = 1, nEns
1430 do eigenVectorColumnIndex = 1, numRetainedEigen
1431 call getModulationFactor( stateVectorMeanInc%vco, levIndex2, &
1432 eigenVectorColumnIndex, numRetainedEigen, &
1433 nEns, vLocalize, &
1434 modulationFactor_r4 )
1435
1436 do memberIndex1 = 1, nEns
1437 ! Compute background ensemble perturbations for the modulated ensemble (Xb_Mod)
1438 pert_r4 = modulationFactor_r4 * ( memberTrl_ptr_r4(memberIndex1,stepIndex,lonIndex,latIndex) - &
1439 meanTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) )
1440
1441 ! Index of the modulated ensemble member corresponding to original
1442 ! ensemble member index (memberIndex1) and eigenVectorColumnIndex.
1443 memberIndexInModEns = (eigenVectorColumnIndex - 1) * nEns + memberIndex1
1444
1445 ! sum Xb_Mod * Wa over all modulated ensembles to get member perturbations for
1446 ! original ensemble (memberIndex2)
1447 memberAnlPert(memberIndex2) = memberAnlPert(memberIndex2) + &
1448 weightsMembers(memberIndexInModEns,memberIndex2,lonIndex,latIndex) * pert_r4
1449 end do
1450 end do
1451
1452 ! Compute final member perturbations by removing background original ensemble perturbations
1453 memberAnlPert(memberIndex2) = (memberTrl_ptr_r4(memberIndex2,stepIndex,lonIndex,latIndex) - &
1454 meanTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex)) + &
1455 memberAnlPert(memberIndex2)
1456
1457 end do ! memberIndex2
1458 else
1459 do memberIndex2 = 1, nEns
1460 do memberIndex1 = 1, nEns
1461 memberAnlPert(memberIndex2) = memberAnlPert(memberIndex2) + &
1462 weightsMembers(memberIndex1,memberIndex2,lonIndex,latIndex) * &
1463 (memberTrl_ptr_r4(memberIndex1,stepIndex,lonIndex,latIndex) - &
1464 meanTrl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex))
1465 end do ! memberIndex1
1466 end do ! memberIndex2
1467 end if
1468
1469 ! Add analysis member perturbation to mean analysis
1470 memberAnl_ptr_r4(:,stepIndex,lonIndex,latIndex) = &
1471 meanAnl_ptr_r4(lonIndex,latIndex,varLevIndex,stepIndex) + memberAnlPert(:)
1472
1473 call utl_tmg_stop(144)
1474
1475 end do ! stepIndex
1476 end do ! varLevIndex
1477
1478 end do LON_LOOP5
1479 end do
1480 !$OMP END PARALLEL DO
1481
1482 call utl_tmg_stop(143)
1483
1484 end do LEV_LOOP
1485
1486 if (countMaxExceeded > 0) then
1487 write(*,*) 'enkf_LETKFanalyses: WARNING: Found more local obs than specified max number at ', &
1488 real(100*countMaxExceeded)/real(numGridPointWeights), '% of grid points.'
1489 write(*,*) ' Maximum number found was ', maxCountMaxExceeded, &
1490 ' which is greater than specified number ', maxNumLocalObs
1491 write(*,*) ' Therefore will keep closest obs only.'
1492 end if
1493
1494 call utl_tmg_start(141,'----Barr')
1495 call rpn_comm_barrier('GRID',ierr)
1496 call utl_tmg_stop(141)
1497
1498 call gsv_deallocate(stateVectorMeanInc)
1499 call gsv_deallocate(stateVectorMeanTrl)
1500
1501 write(*,*) 'enkf_LETKFanalyses: done'
1502 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1503
1504 call utl_tmg_stop(131)
1505
1506 end subroutine enkf_LETKFanalyses
1507
1508 !----------------------------------------------------------------------
1509 ! enkf_computeVertLocation (private subroutine)
1510 !----------------------------------------------------------------------
1511 subroutine enkf_computeVertLocation(vertLocation_r4,stateVectorMeanTrl)
1512 !
1513 !:Purpose: Compute extract global 3D vertical location field from supplied
1514 ! stateVector. Can be either logPressure or depth levels.
1515 !
1516 implicit none
1517
1518 ! Arguments:
1519 real(4), allocatable, intent(inout) :: vertLocation_r4(:,:,:)
1520 type(struct_gsv), intent(inout) :: stateVectorMeanTrl
1521
1522 ! Locals:
1523 integer :: nLev_M, nLev_depth, nLev_vertLocation, levIndex, nsize, ierr
1524 real(4), pointer :: vertLocation_ptr_r4(:,:,:)
1525 type(struct_gsv) :: stateVectorMeanTrlPressure
1526 type(struct_gsv) :: stateVectorMeanTrlPressure_1step
1527
1528 write(*,*) 'enkf_computeVertLocation: starting'
1529
1530 nLev_M = gsv_getNumLev(stateVectorMeanTrl, 'MM')
1531 nLev_depth = gsv_getNumLev(stateVectorMeanTrl, 'DP')
1532 if ( nLev_M > 0 .and. nLev_depth > 0 ) then
1533 call utl_abort('enkf_computeVertLocation: both momentum and depth levels exist.')
1534 else if ( nLev_M == 0 .and. nLev_depth == 0 ) then
1535 call utl_abort('enkf_computeVertLocation: neither momentum nor depth levels exist.')
1536 end if
1537 nLev_vertLocation = max(nLev_M, nLev_depth)
1538
1539 allocate(vertLocation_r4(stateVectorMeanTrl%hco%ni, &
1540 stateVectorMeanTrl%hco%nj, &
1541 nLev_vertLocation))
1542
1543 if ( nLev_M > 0 ) then ! log pressure for NWP fields
1544
1545 ! Compute background ens mean 3D log pressure and make mpiglobal for vertical localization
1546 call gsv_allocate( stateVectorMeanTrlPressure, tim_nstepobsinc, &
1547 stateVectorMeanTrl%hco, stateVectorMeanTrl%vco, dateStamp_opt=tim_getDateStamp(), &
1548 mpi_local_opt=.true., mpi_distribution_opt='Tiles', &
1549 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0','P_M','P_T'/) )
1550 call gsv_zero(stateVectorMeanTrlPressure)
1551 call gsv_copy(stateVectorMeanTrl, stateVectorMeanTrlPressure, allowVarMismatch_opt=.true.)
1552 call gvt_transform(stateVectorMeanTrlPressure,'ZandP_nl')
1553 if (mmpi_myid == 0) then
1554 call gsv_allocate( stateVectorMeanTrlPressure_1step, 1, &
1555 stateVectorMeanTrl%hco, stateVectorMeanTrl%vco, dateStamp_opt=tim_getDateStamp(), &
1556 mpi_local_opt=.false., &
1557 dataKind_opt=4, allocHeightSfc_opt=.true., varNames_opt=(/'P0','P_M','P_T'/) )
1558 end if
1559 call gsv_transposeTilesToStep(stateVectorMeanTrlPressure_1step, stateVectorMeanTrlPressure, (tim_nstepobsinc+1)/2)
1560 call gsv_deallocate(stateVectorMeanTrlPressure)
1561 if (mmpi_myid == 0) then
1562 call gsv_getField(stateVectorMeanTrlPressure_1step,vertLocation_ptr_r4,'P_M')
1563 vertLocation_r4(:,:,:) = log(vertLocation_ptr_r4(:,:,:))
1564 end if
1565 nsize = stateVectorMeanTrlPressure%ni * stateVectorMeanTrlPressure%nj * nLev_M
1566 call rpn_comm_bcast(vertLocation_r4, nsize, 'mpi_real4', 0, 'GRID', ierr)
1567
1568 else if ( nLev_depth > 0 ) then ! depth for ocean fields
1569
1570 ! fill in all horizontal grid points with the same profile of depth values
1571 do levIndex = 1, nLev_depth
1572 write(*,*) 'setting vertLocation for levIndex =', levIndex, &
1573 ', depth = ', stateVectorMeanTrl%vco%depths(levIndex)
1574 vertLocation_r4(:,:,levIndex) = stateVectorMeanTrl%vco%depths(levIndex)
1575 end do
1576
1577 end if
1578
1579 write(*,*) 'enkf_computeVertLocation: finished'
1580
1581 end subroutine enkf_computeVertLocation
1582
1583 !----------------------------------------------------------------------
1584 ! enkf_setupMpiDistribution (private subroutine)
1585 !----------------------------------------------------------------------
1586 subroutine enkf_LETKFsetupMpiDistribution(myNumLatLonRecv, myNumLatLonSend, &
1587 myLatIndexesRecv, myLonIndexesRecv, &
1588 myLatIndexesSend, myLonIndexesSend, &
1589 myProcIndexesRecv, myProcIndexesSend, &
1590 myNumProcIndexesSend, mpiDistribution, wInterpInfo)
1591 !
1592 !:Purpose: Setup for distribution of grid points over mpi tasks.
1593 !
1594 implicit none
1595
1596 ! Arguments:
1597 integer, intent(out) :: myNumLatLonRecv
1598 integer, intent(out) :: myNumLatLonSend
1599 integer, allocatable, intent(out) :: myLatIndexesRecv(:)
1600 integer, allocatable, intent(out) :: myLonIndexesRecv(:)
1601 integer, allocatable, intent(out) :: myLatIndexesSend(:)
1602 integer, allocatable, intent(out) :: myLonIndexesSend(:)
1603 integer, allocatable, intent(out) :: myProcIndexesRecv(:)
1604 integer, allocatable, intent(out) :: myProcIndexesSend(:,:)
1605 integer, allocatable, intent(out) :: myNumProcIndexesSend(:)
1606 character(len=*), intent(in) :: mpiDistribution
1607 type(struct_enkfInterpInfo), intent(in) :: wInterpInfo
1608
1609 ! Locals:
1610 integer :: latIndex, lonIndex, procIndex, procIndexSend, latLonIndex
1611 integer :: myLonBegHalo, myLonEndHalo, myLatBegHalo, myLatEndHalo
1612 integer :: numLatLonRecvMax, numLatLonTotalUnique, ierr
1613 integer, allocatable :: allLatIndexesRecv(:,:), allLonIndexesRecv(:,:)
1614 integer, allocatable :: allLatIndexesSend(:,:), allLonIndexesSend(:,:)
1615 integer, allocatable :: allNumLatLonRecv(:), allNumLatLonSend(:)
1616
1617 myLonBegHalo = wInterpInfo%myLonBegHalo
1618 myLonEndHalo = wInterpInfo%myLonEndHalo
1619 myLatBegHalo = wInterpInfo%myLatBegHalo
1620 myLatEndHalo = wInterpInfo%myLatEndHalo
1621
1622 write(*,*) 'enkf_LETKFsetupMpiDistribution: starting'
1623 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1624
1625 if (trim(mpiDistribution) == 'TILES') then
1626
1627 ! First, determine number of grid points needed locally (for recv-ing)
1628 myNumLatLonRecv = 0
1629 do latIndex = myLatBegHalo, myLatEndHalo
1630 LON_LOOP0: do lonIndex = myLonBegHalo, myLonEndHalo
1631 ! If this lat-lon is to be interpolated, then skip calculation
1632 if (wInterpInfo%numIndexes(lonIndex,latIndex) > 0) cycle LON_LOOP0
1633 myNumLatLonRecv = myNumLatLonRecv + 1
1634 end do LON_LOOP0
1635 end do
1636
1637 write(*,*) 'enkf_LETKFsetupMpiDistribution: myNumLatLonRecv =', myNumLatLonRecv
1638
1639 ! Determine list of grid point indexes where weights needed locally (for recv-ing)
1640 allocate(myLatIndexesRecv(myNumLatLonRecv))
1641 allocate(myLonIndexesRecv(myNumLatLonRecv))
1642 allocate(myProcIndexesRecv(myNumLatLonRecv))
1643 myNumLatLonRecv = 0
1644 do latIndex = myLatBegHalo, myLatEndHalo
1645 LON_LOOP1: do lonIndex = myLonBegHalo, myLonEndHalo
1646 ! If this lat-lon is to be interpolated, then skip calculation
1647 if (wInterpInfo%numIndexes(lonIndex,latIndex) > 0) cycle LON_LOOP1
1648 myNumLatLonRecv = myNumLatLonRecv + 1
1649
1650 myLatIndexesRecv(myNumLatLonRecv) = latIndex
1651 myLonIndexesRecv(myNumLatLonRecv) = lonIndex
1652 myProcIndexesRecv(myNumLatLonRecv) = mmpi_myid+1
1653 end do LON_LOOP1
1654 end do
1655
1656 ! No communication, so send info equals recv info
1657 myNumLatLonSend = myNumLatLonRecv
1658 allocate(myLatIndexesSend(myNumLatLonSend))
1659 allocate(myLonIndexesSend(myNumLatLonSend))
1660 allocate(myProcIndexesSend(myNumLatLonSend,1))
1661 allocate(myNumProcIndexesSend(myNumLatLonSend))
1662
1663 myLatIndexesSend(:) = myLatIndexesRecv(:)
1664 myLonIndexesSend(:) = myLonIndexesRecv(:)
1665 myProcIndexesSend(:,1) = myProcIndexesRecv(:)
1666 myNumProcIndexesSend(:) = 1
1667
1668 else if (trim(mpiDistribution) == 'ROUNDROBIN') then
1669
1670 ! First, determine number of grid points needed locally (for recv-ing)
1671 myNumLatLonRecv = 0
1672 do latIndex = myLatBegHalo, myLatEndHalo
1673 LON_LOOP2: do lonIndex = myLonBegHalo, myLonEndHalo
1674 ! If this lat-lon is to be interpolated, then skip calculation
1675 if (wInterpInfo%numIndexes(lonIndex,latIndex) > 0) cycle LON_LOOP2
1676 myNumLatLonRecv = myNumLatLonRecv + 1
1677 end do LON_LOOP2
1678 end do
1679
1680 ! Communicate to all mpi tasks
1681 allocate(allNumLatLonRecv(mmpi_nprocs))
1682 call rpn_comm_allgather(myNumLatLonRecv, 1, "mpi_integer", &
1683 allNumLatLonRecv, 1,"mpi_integer", "GRID", ierr)
1684 numLatLonRecvMax = maxval(allNumLatLonRecv)
1685 write(*,*) 'enkf_LETKFsetupMpiDistribution: allNumLatLonRecv =', allNumLatLonRecv(:)
1686 write(*,*) 'enkf_LETKFsetupMpiDistribution: numLatLonRecvSum =', sum(allNumLatLonRecv)
1687 write(*,*) 'enkf_LETKFsetupMpiDistribution: numLatLonRecvMax =', numLatLonRecvMax
1688
1689 ! Determine list of grid point indexes where weights needed locally (for recv-ing)
1690 allocate(myLatIndexesRecv(numLatLonRecvMax))
1691 allocate(myLonIndexesRecv(numLatLonRecvMax))
1692 allocate(myProcIndexesRecv(numLatLonRecvMax))
1693 myLatIndexesRecv(:) = -1
1694 myLonIndexesRecv(:) = -1
1695 myProcIndexesRecv(:) = -1
1696 myNumLatLonRecv = 0
1697 do latIndex = myLatBegHalo, myLatEndHalo
1698 LON_LOOP3: do lonIndex = myLonBegHalo, myLonEndHalo
1699 ! If this lat-lon is to be interpolated, then skip calculation
1700 if (wInterpInfo%numIndexes(lonIndex,latIndex) > 0) cycle LON_LOOP3
1701 myNumLatLonRecv = myNumLatLonRecv + 1
1702
1703 myLatIndexesRecv(myNumLatLonRecv) = latIndex
1704 myLonIndexesRecv(myNumLatLonRecv) = lonIndex
1705 end do LON_LOOP3
1706 end do
1707
1708 ! Communicate to all mpi tasks this list of grid point lat-lon indexes
1709 allocate(allLatIndexesRecv(numLatLonRecvMax, mmpi_nprocs))
1710 allocate(allLonIndexesRecv(numLatLonRecvMax, mmpi_nprocs))
1711 call rpn_comm_allgather(myLatIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1712 allLatIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1713 "GRID", ierr)
1714 call rpn_comm_allgather(myLonIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1715 allLonIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1716 "GRID", ierr)
1717
1718 ! From these lat-lon lists, create unique master list of all grid points where weights computed
1719 ! and assign to mpi tasks for doing the calculation and for send-ing
1720 allocate(myLatIndexesSend(numLatLonRecvMax))
1721 allocate(myLonIndexesSend(numLatLonRecvMax))
1722 myLatIndexesSend(:) = -1
1723 myLonIndexesSend(:) = -1
1724 numLatLonTotalUnique = 0
1725 myNumLatLonSend = 0
1726 do procIndex = 1, mmpi_nprocs
1727 WEIGHTS1LEV_LOOP: do latLonIndex = 1, allNumLatLonRecv(procIndex)
1728 if (enkf_latLonAlreadyFound(allLatIndexesRecv, allLonIndexesRecv, latLonIndex, procIndex)) &
1729 cycle WEIGHTS1LEV_LOOP
1730 ! Count the total number of weights
1731 numLatLonTotalUnique = numLatLonTotalUnique + 1
1732
1733 ! Round-robin distribution of master list across mpi tasks
1734 procIndexSend = 1 + mod(numLatLonTotalUnique-1, mmpi_nprocs)
1735
1736 ! Store the lat-lon indexes of the weights I am responsible for
1737 if (procIndexSend == (mmpi_myid+1)) then
1738 myNumLatLonSend = myNumLatLonSend + 1
1739 myLatIndexesSend(myNumLatLonSend) = &
1740 allLatIndexesRecv(latLonIndex, procIndex)
1741 myLonIndexesSend(myNumLatLonSend) = &
1742 allLonIndexesRecv(latLonIndex, procIndex)
1743 end if
1744 end do WEIGHTS1LEV_LOOP
1745 end do
1746 write(*,*) 'enkf_LETKFsetupMpiDistribution: number of lat/lon points where weights to be computed =', &
1747 numLatLonTotalUnique
1748
1749 ! Communicate to all mpi tasks this list of grid point lat-lon indexes
1750 allocate(allNumLatLonSend(mmpi_nprocs))
1751 call rpn_comm_allgather(myNumLatLonSend, 1, "mpi_integer", &
1752 allNumLatLonSend, 1,"mpi_integer", "GRID", ierr)
1753 allocate(allLatIndexesSend(numLatLonRecvMax, mmpi_nprocs))
1754 allocate(allLonIndexesSend(numLatLonRecvMax, mmpi_nprocs))
1755 call rpn_comm_allgather(myLatIndexesSend, numLatLonRecvMax, "mpi_integer", &
1756 allLatIndexesSend, numLatLonRecvMax, "mpi_integer", &
1757 "GRID", ierr)
1758 call rpn_comm_allgather(myLonIndexesSend, numLatLonRecvMax, "mpi_integer", &
1759 allLonIndexesSend, numLatLonRecvMax, "mpi_integer", &
1760 "GRID", ierr)
1761
1762 ! Figure out which mpi tasks I will need to send my results to
1763 allocate(myProcIndexesSend(myNumLatLonSend,mmpi_nprocs))
1764 allocate(myNumProcIndexesSend(myNumLatLonSend))
1765 myProcIndexesSend(:,:) = -1
1766 myNumProcIndexesSend(:) = 0
1767 do latLonIndex = 1, myNumLatLonSend
1768 do procIndex = 1, mmpi_nprocs
1769 if ( any( (myLatIndexesSend(latLonIndex) == allLatIndexesRecv(1:allNumLatLonRecv(procIndex), procIndex)) .and. &
1770 (myLonIndexesSend(latLonIndex) == allLonIndexesRecv(1:allNumLatLonRecv(procIndex), procIndex)) ) ) then
1771 myNumProcIndexesSend(latLonIndex) = myNumProcIndexesSend(latLonIndex) + 1
1772 myProcIndexesSend(latLonIndex,myNumProcIndexesSend(latLonIndex)) = procIndex
1773 end if
1774 end do
1775 end do
1776
1777 ! Figure out which mpi tasks I will receive the results from
1778 do latLonIndex = 1, myNumLatLonRecv
1779 do procIndex = 1, mmpi_nprocs
1780 if ( any( (myLatIndexesRecv(latLonIndex) == allLatIndexesSend(1:allNumLatLonSend(procIndex), procIndex)) .and. &
1781 (myLonIndexesRecv(latLonIndex) == allLonIndexesSend(1:allNumLatLonSend(procIndex), procIndex)) ) ) then
1782 myProcIndexesRecv(latLonIndex) = procIndex
1783 end if
1784 end do
1785 end do
1786
1787 else
1788 call utl_abort('enkf_LETKFsetupMpiDistribution: unknown MPI distribution selected')
1789 end if
1790
1791 write(*,*) 'enkf_LETKFsetupMpiDistribution: lat/lon/proc indexes I need to receive:'
1792 do latLonIndex = 1, myNumLatLonRecv
1793 write(*,*) myLatIndexesRecv(latLonIndex), myLonIndexesRecv(latLonIndex), &
1794 myProcIndexesRecv(latLonIndex)
1795 end do
1796
1797 write(*,*) 'enkf_LETKFsetupMpiDistribution: number of lat/lon indexes I am responsible for =', myNumLatLonSend
1798 write(*,*) 'enkf_LETKFsetupMpiDistribution: the lat/lon/proc indexes I am responsible for:'
1799 do latLonIndex = 1, myNumLatLonSend
1800 write(*,*) myLatIndexesSend(latLonIndex), myLonIndexesSend(latLonIndex), &
1801 myProcIndexesSend(latLonIndex,1:myNumProcIndexesSend(latLonIndex))
1802 end do
1803
1804 write(*,*) 'enkf_LETKFsetupMpiDistribution: done'
1805 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1806
1807 end subroutine enkf_LETKFsetupMpiDistribution
1808
1809 !----------------------------------------------------------------------
1810 ! enkf_LETKFgetMpiGlobalTags (private subroutine)
1811 !----------------------------------------------------------------------
1812 subroutine enkf_LETKFgetMpiGlobalTags(latLonTagMpiGlobal,myLatIndexesRecv,myLonIndexesRecv)
1813 implicit none
1814
1815 ! Arguments:
1816 integer, intent(out) :: latLonTagMpiGlobal(:,:)
1817 integer, intent(in) :: myLatIndexesRecv(:)
1818 integer, intent(in) :: myLonIndexesRecv(:)
1819
1820 ! Locals:
1821 integer :: ierr, ni, nj, lonIndex, latIndex
1822 integer :: countTags, myNumLatLonRecv, numLatLonRecvMax
1823 integer, allocatable :: allNumLatLonRecv(:)
1824 integer, allocatable :: allLatIndexesRecv(:,:), allLonIndexesRecv(:,:)
1825
1826 write(*,*) 'enkf_LETKFgetMpiGlobalTags: Starting'
1827
1828 ni = size(latLonTagMpiGlobal,1)
1829 nj = size(latLonTagMpiGlobal,2)
1830
1831 myNumLatLonRecv = size(myLatIndexesRecv)
1832 allocate(allNumLatLonRecv(mmpi_nprocs))
1833 call rpn_comm_allgather(myNumLatLonRecv, 1, "mpi_integer", &
1834 allNumLatLonRecv, 1,"mpi_integer", "GRID", ierr)
1835 numLatLonRecvMax = maxval(allNumLatLonRecv)
1836
1837 ! Communicate to all mpi tasks this list of grid point lat-lon indexes
1838 allocate(allLatIndexesRecv(numLatLonRecvMax, mmpi_nprocs))
1839 allocate(allLonIndexesRecv(numLatLonRecvMax, mmpi_nprocs))
1840 call rpn_comm_allgather(myLatIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1841 allLatIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1842 "GRID", ierr)
1843 call rpn_comm_allgather(myLonIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1844 allLonIndexesRecv, numLatLonRecvMax, "mpi_integer", &
1845 "GRID", ierr)
1846
1847 latLonTagMpiGlobal(:,:) = 0
1848 !$OMP PARALLEL DO PRIVATE(latIndex, lonIndex)
1849 do lonIndex = 1, ni
1850 do latIndex = 1, nj
1851 if (any(lonIndex == allLonIndexesRecv(:,:) .and. latIndex == allLatIndexesRecv(:,:))) then
1852 latLonTagMpiGlobal(lonIndex,latIndex) = 1
1853 end if
1854 end do
1855 end do
1856 !$OMP END PARALLEL DO
1857
1858 countTags = 0
1859 do lonIndex = 1, ni
1860 do latIndex = 1, nj
1861 if (latLonTagMpiGlobal(lonIndex,latIndex) == 1) then
1862 countTags = countTags + 1
1863 latLonTagMpiGlobal(lonIndex,latIndex) = countTags
1864 end if
1865 end do
1866 end do
1867 write(*,*) 'number of Recv grid points found = ', maxval(latLonTagMpiGlobal(:,:))
1868
1869 write(*,*) 'enkf_LETKFgetMpiGlobalTags: Finished'
1870
1871 end subroutine enkf_LETKFgetMpiGlobalTags
1872
1873 !----------------------------------------------------------------------
1874 ! enkf_latLonAlreadyFound (private function)
1875 !----------------------------------------------------------------------
1876 function enkf_latLonAlreadyFound(allLatIndexesRecv, allLonIndexesRecv, latLonIndex, procIndex) result(found)
1877 implicit none
1878
1879 ! Arguments:
1880 integer, intent(in) :: allLatIndexesRecv(:,:)
1881 integer, intent(in) :: allLonIndexesRecv(:,:)
1882 integer, intent(in) :: latLonIndex
1883 integer, intent(in) :: procIndex
1884 ! Result:
1885 logical :: found
1886
1887 ! Locals:
1888 integer :: latLonIndex2, procIndex2, numLatLonRecvMax
1889
1890 numLatLonRecvMax = size(allLatIndexesRecv, 1)
1891
1892 ! check on all previous mpi tasks if this lat/lon has already been encountered
1893 found = .false.
1894 do procIndex2 = 1, procIndex-1
1895 WEIGHTS1LEV_LOOP2: do latLonIndex2 = 1, numLatLonRecvMax
1896 if (allLatIndexesRecv(latLonIndex2, procIndex2) < 0) cycle WEIGHTS1LEV_LOOP2
1897 if ( (allLatIndexesRecv(latLonIndex, procIndex) == &
1898 allLatIndexesRecv(latLonIndex2, procIndex2)) .and. &
1899 (allLonIndexesRecv(latLonIndex, procIndex) == &
1900 allLonIndexesRecv(latLonIndex2, procIndex2)) ) then
1901 found = .true.
1902 exit WEIGHTS1LEV_LOOP2
1903 end if
1904 end do WEIGHTS1LEV_LOOP2
1905 end do
1906
1907 end function enkf_latLonAlreadyFound
1908
1909 !--------------------------------------------------------------------------
1910 ! enkf_setupInterpInfo
1911 !--------------------------------------------------------------------------
1912 subroutine enkf_setupInterpInfo(wInterpInfo, hco, weightLatLonStep, &
1913 myLonBeg,myLonEnd,myLatBeg,myLatEnd)
1914 !
1915 !:Purpose: Setup the weights and lat/lon indices needed to bilinearly
1916 ! interpolate the LETKF weights from a coarse grid to the full
1917 ! resolution grid. The coarseness of the grid is specified by
1918 ! the weightLatLonStep argument.
1919 !
1920 implicit none
1921
1922 ! Arguments:
1923 type(struct_enkfInterpInfo), intent(out) :: wInterpInfo
1924 type(struct_hco), intent(in) :: hco
1925 integer, intent(in) :: weightLatLonStep
1926 integer, intent(in) :: myLonBeg
1927 integer, intent(in) :: myLonEnd
1928 integer, intent(in) :: myLatBeg
1929 integer, intent(in) :: myLatEnd
1930
1931 ! Locals:
1932 integer :: lonIndex, latIndex, ni, nj
1933 integer :: myLonBegHalo, myLonEndHalo, myLatBegHalo, myLatEndHalo
1934 real(8) :: interpWeightLon, interpWeightLat
1935 logical :: includesYinYangBndry
1936
1937 ni = hco%ni
1938 nj = hco%nj
1939
1940 myLonBegHalo = 1 + weightLatLonStep * floor(real(myLonBeg - 1)/real(weightLatLonStep))
1941 myLonEndHalo = min(ni, 1 + weightLatLonStep * ceiling(real(myLonEnd - 1)/real(weightLatLonStep)))
1942 myLatBegHalo = 1 + weightLatLonStep * floor(real(myLatBeg - 1)/real(weightLatLonStep))
1943 myLatEndHalo = min(nj, 1 + weightLatLonStep * ceiling(real(myLatEnd - 1)/real(weightLatLonStep)))
1944 write(*,*) 'enkf_setupInterpInfo: myLonBeg/End, myLatBeg/End (original) = ', &
1945 myLonBeg, myLonEnd, myLatBeg, myLatEnd
1946 write(*,*) 'enkf_setupInterpInfo: myLonBeg/End, myLatBeg/End (with Halo) = ', &
1947 myLonBegHalo, myLonEndHalo, myLatBegHalo, myLatEndHalo
1948 write(*,*) 'enkf_setupInterpInfo: myLonCount, myLatCount (with Halo) = ', &
1949 myLonEndHalo-myLonBegHalo+1, myLatEndHalo-myLatBegHalo+1
1950 write(*,*) 'enkf_setupInterpInfo: number of local gridpts where weights computed = ', &
1951 ( 1 + ceiling(real(myLonEndHalo - myLonBegHalo) / real(weightLatLonStep)) ) * &
1952 ( 1 + ceiling(real(myLatEndHalo - myLatBegHalo) / real(weightLatLonStep)) )
1953 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1954
1955 wInterpInfo%latLonStep = weightLatLonStep
1956 wInterpInfo%myLonBegHalo = myLonBegHalo
1957 wInterpInfo%myLonEndHalo = myLonEndHalo
1958 wInterpInfo%myLatBegHalo = myLatBegHalo
1959 wInterpInfo%myLatEndHalo = myLatEndHalo
1960 wInterpInfo%myLonBeg = myLonBeg
1961 wInterpInfo%myLonEnd = myLonEnd
1962 wInterpInfo%myLatBeg = myLatBeg
1963 wInterpInfo%myLatEnd = myLatEnd
1964
1965 allocate(wInterpInfo%numIndexes(myLonBegHalo:myLonEndHalo,myLatBegHalo:myLatEndHalo))
1966 if (weightLatLonStep > 1) then
1967 ! Figure out if this tile straddles Yin-Yang boundary
1968 if (hco%grtyp == 'U' .and. myLatBegHalo <= nj/2 .and. myLatEndHalo >= ((nj/2)+1)) then
1969 includesYinYangBndry = .true.
1970 else
1971 includesYinYangBndry = .false.
1972 end if
1973 allocate(wInterpInfo%lonIndexes(4,myLonBegHalo:myLonEndHalo,myLatBegHalo:myLatEndHalo))
1974 allocate(wInterpInfo%latIndexes(4,myLonBegHalo:myLonEndHalo,myLatBegHalo:myLatEndHalo))
1975 allocate(wInterpInfo%interpWeights(4,myLonBegHalo:myLonEndHalo,myLatBegHalo:myLatEndHalo))
1976 wInterpInfo%lonIndexes(:,:,:) = 0
1977 wInterpInfo%latIndexes(:,:,:) = 0
1978 wInterpInfo%interpWeights(:,:,:) = 0.0D0
1979 ! Determine which lat-lon are interpolated (wInterpInfo%numIndexes>0)
1980 wInterpInfo%numIndexes(:,:) = 4
1981 do latIndex = myLatBegHalo, myLatEndHalo, weightLatLonStep
1982 do lonIndex = myLonBegHalo, myLonEndHalo, weightLatLonStep
1983 wInterpInfo%numIndexes(lonIndex,latIndex) = 0
1984 end do
1985 end do
1986 ! Ensure weights are computed along edge of domain
1987 if (myLonEndHalo == ni .and. &
1988 myLatEndHalo == nj) then
1989 wInterpInfo%numIndexes(myLonEndHalo,myLatEndHalo) = 0
1990 end if
1991 if (myLonEndHalo == ni) then
1992 do latIndex = myLatBegHalo, myLatEndHalo, weightLatLonStep
1993 wInterpInfo%numIndexes(myLonEndHalo,latIndex) = 0
1994 end do
1995 ! Ensure weights are computed along both sides of Yin-Yang boundary
1996 if (includesYinYangBndry) then
1997 wInterpInfo%numIndexes(ni,nj/2) = 0
1998 wInterpInfo%numIndexes(ni,(nj/2)+1) = 0
1999 write(*,*) 'enkf_setupInterpInfo: Yin-Yang boundary (lon,lat1,lat2) =', &
2000 ni, nj/2, (nj/2)+1
2001 end if
2002 end if
2003 if (myLatEndHalo == nj) then
2004 do lonIndex = myLonBegHalo, myLonEndHalo, weightLatLonStep
2005 wInterpInfo%numIndexes(lonIndex,myLatEndHalo) = 0
2006 end do
2007 end if
2008 ! Ensure weights are computed along both sides of Yin-Yang boundary
2009 if (includesYinYangBndry) then
2010 do lonIndex = myLonBegHalo, myLonEndHalo, weightLatLonStep
2011 wInterpInfo%numIndexes(lonIndex,nj/2) = 0
2012 wInterpInfo%numIndexes(lonIndex,(nj/2)+1) = 0
2013 write(*,*) 'enkf_setupInterpInfo: Yin-Yang boundary (lon,lat1,lat2) =', &
2014 lonIndex, nj/2, (nj/2)+1
2015 end do
2016 end if
2017 ! For lon-only interpolation
2018 do latIndex = myLatBegHalo, myLatEndHalo
2019 if (wInterpInfo%numIndexes(myLonBegHalo,latIndex) > 0) cycle
2020 do lonIndex = myLonBegHalo, myLonEndHalo
2021 if (wInterpInfo%numIndexes(lonIndex,latIndex) == 0) cycle
2022 ! Find nearest grid point with a value towards left
2023 wInterpInfo%numIndexes(lonIndex,latIndex) = 2
2024 wInterpInfo%lonIndexes(1,lonIndex,latIndex) = myLonBegHalo + &
2025 weightLatLonStep * floor(real(lonIndex - myLonBegHalo)/real(weightLatLonStep))
2026 wInterpInfo%lonIndexes(2,lonIndex,latIndex) = min(ni, &
2027 wInterpInfo%lonIndexes(1,lonIndex,latIndex) + weightLatLonStep)
2028 wInterpInfo%latIndexes(1,lonIndex,latIndex) = latIndex
2029 wInterpInfo%latIndexes(2,lonIndex,latIndex) = latIndex
2030 wInterpInfo%interpWeights(1,lonIndex,latIndex) = &
2031 real(wInterpInfo%lonIndexes(2,lonIndex,latIndex) - lonIndex, 8)/real(weightLatLonStep, 8)
2032 wInterpInfo%interpWeights(2,lonIndex,latIndex) = 1.0D0 - &
2033 wInterpInfo%interpWeights(1,lonIndex,latIndex)
2034 end do
2035 end do
2036 ! For lat-only interpolation
2037 do latIndex = myLatBegHalo, myLatEndHalo
2038 do lonIndex = myLonBegHalo, myLonEndHalo, weightLatLonStep
2039 if (wInterpInfo%numIndexes(lonIndex,myLatBegHalo) > 0) cycle
2040 if (wInterpInfo%numIndexes(lonIndex,latIndex) == 0) cycle
2041 ! Find nearest grid point with a value towards bottom
2042 wInterpInfo%numIndexes(lonIndex,latIndex) = 2
2043 wInterpInfo%lonIndexes(1,lonIndex,latIndex) = lonIndex
2044 wInterpInfo%lonIndexes(2,lonIndex,latIndex) = lonIndex
2045 wInterpInfo%latIndexes(1,lonIndex,latIndex) = myLatBegHalo + &
2046 weightLatLonStep * floor(real(latIndex - myLatBegHalo)/real(weightLatLonStep))
2047 wInterpInfo%latIndexes(2,lonIndex,latIndex) = min(nj, &
2048 wInterpInfo%latIndexes(1,lonIndex,latIndex) + weightLatLonStep)
2049 ! Ensure we do not interpolate values across Yin-Yang boundary
2050 if (includesYinYangBndry) then
2051 if (latIndex <= nj/2) then
2052 wInterpInfo%latIndexes(2,lonIndex,latIndex) = min(nj/2, wInterpInfo%latIndexes(2,lonIndex,latIndex))
2053 else if(latIndex >= (nj/2)+1) then
2054 wInterpInfo%latIndexes(1,lonIndex,latIndex) = max((nj/2)+1, wInterpInfo%latIndexes(1,lonIndex,latIndex))
2055 end if
2056 end if
2057 wInterpInfo%interpWeights(1,lonIndex,latIndex) = &
2058 real(wInterpInfo%latIndexes(2,lonIndex,latIndex) - latIndex, 8)/real(weightLatLonStep, 8)
2059 wInterpInfo%interpWeights(2,lonIndex,latIndex) = 1.0D0 - &
2060 wInterpInfo%interpWeights(1,lonIndex,latIndex)
2061 end do
2062 end do
2063 ! For interior points needing 2D interpolation
2064 do latIndex = myLatBegHalo, myLatEndHalo
2065 do lonIndex = myLonBegHalo, myLonEndHalo
2066 if (wInterpInfo%numIndexes(lonIndex,latIndex) == 0) cycle ! no interpolation
2067 if (wInterpInfo%lonIndexes(1,lonIndex,latIndex) /= 0) cycle ! already set up
2068 wInterpInfo%numIndexes(lonIndex,latIndex) = 4
2069 ! 1. bottom-left indexes
2070 wInterpInfo%lonIndexes(1,lonIndex,latIndex) = myLonBegHalo + &
2071 weightLatLonStep * floor(real(lonIndex - myLonBegHalo)/real(weightLatLonStep))
2072 wInterpInfo%latIndexes(1,lonIndex,latIndex) = myLatBegHalo + &
2073 weightLatLonStep * floor(real(latIndex - myLatBegHalo)/real(weightLatLonStep))
2074 ! 2. bottom-right indexes
2075 wInterpInfo%lonIndexes(2,lonIndex,latIndex) = min(ni, &
2076 wInterpInfo%lonIndexes(1,lonIndex,latIndex) + weightLatLonStep)
2077 wInterpInfo%latIndexes(2,lonIndex,latIndex) = wInterpInfo%latIndexes(1,lonIndex,latIndex)
2078 ! 3. top-left indexes
2079 wInterpInfo%lonIndexes(3,lonIndex,latIndex) = wInterpInfo%lonIndexes(1,lonIndex,latIndex)
2080 wInterpInfo%latIndexes(3,lonIndex,latIndex) = min(nj, &
2081 wInterpInfo%latIndexes(1,lonIndex,latIndex) + weightLatLonStep)
2082 ! 4. top-right indexes
2083 wInterpInfo%lonIndexes(4,lonIndex,latIndex) = wInterpInfo%lonIndexes(2,lonIndex,latIndex)
2084 wInterpInfo%latIndexes(4,lonIndex,latIndex) = wInterpInfo%latIndexes(3,lonIndex,latIndex)
2085 ! Ensure we do not interpolate values across Yin-Yang boundary
2086 if (includesYinYangBndry) then
2087 if (latIndex <= nj/2) then
2088 wInterpInfo%latIndexes(3,lonIndex,latIndex) = min(nj/2, wInterpInfo%latIndexes(3,lonIndex,latIndex))
2089 wInterpInfo%latIndexes(4,lonIndex,latIndex) = min(nj/2, wInterpInfo%latIndexes(4,lonIndex,latIndex))
2090 else if(latIndex >= (nj/2)+1) then
2091 wInterpInfo%latIndexes(1,lonIndex,latIndex) = max((nj/2)+1, wInterpInfo%latIndexes(1,lonIndex,latIndex))
2092 wInterpInfo%latIndexes(2,lonIndex,latIndex) = max((nj/2)+1, wInterpInfo%latIndexes(2,lonIndex,latIndex))
2093 end if
2094 end if
2095 ! one-dimensional weights in lon and lat directions
2096 interpWeightLon = real(wInterpInfo%lonIndexes(4,lonIndex,latIndex) - lonIndex, 8) / &
2097 real(weightLatLonStep, 8)
2098 interpWeightLat = real(wInterpInfo%latIndexes(4,lonIndex,latIndex) - latIndex, 8) / &
2099 real(weightLatLonStep, 8)
2100 ! four interpolation weights
2101 wInterpInfo%interpWeights(1,lonIndex,latIndex) = interpWeightLon * interpWeightLat
2102 wInterpInfo%interpWeights(2,lonIndex,latIndex) = (1.0D0 - interpWeightLon) * interpWeightLat
2103 wInterpInfo%interpWeights(3,lonIndex,latIndex) = interpWeightLon * (1.0D0 - interpWeightLat)
2104 wInterpInfo%interpWeights(4,lonIndex,latIndex) = (1.0D0 - interpWeightLon) * (1.0D0 - interpWeightLat)
2105 end do
2106 end do
2107 else
2108 ! no interpolation, all weights are computed
2109 wInterpInfo%numIndexes(:,:) = 0
2110 end if
2111 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2112
2113 end subroutine enkf_setupInterpInfo
2114
2115 !--------------------------------------------------------------------------
2116 ! enkf_interpWeights
2117 !--------------------------------------------------------------------------
2118 subroutine enkf_interpWeights(wInterpInfo, weights)
2119 !
2120 !:Purpose: Perform the bilinear interpolation of the weights
2121 ! using the precalculated interpolation info.
2122 !
2123 implicit none
2124
2125 ! Arguments:
2126 type(struct_enkfInterpInfo), intent(in) :: wInterpInfo
2127 real(8), intent(out) :: weights(1:,1:,wInterpInfo%myLonBegHalo:,wInterpInfo%myLatBegHalo:)
2128
2129 ! Locals:
2130 integer :: myLonBegHalo, myLonEndHalo, myLatBegHalo, myLatEndHalo
2131 integer :: lonIndex, latIndex, memberIndex1, memberIndex2, interpIndex
2132 integer :: interpLonIndex, interpLatIndex, numMembers1, numMembers2
2133 integer :: totalCount(mmpi_numthread)
2134 integer, external :: omp_get_thread_num
2135 logical, save :: firstCall = .true.
2136
2137 myLonBegHalo = wInterpInfo%myLonBegHalo
2138 myLonEndHalo = wInterpInfo%myLonEndHalo
2139 myLatBegHalo = wInterpInfo%myLatBegHalo
2140 myLatEndHalo = wInterpInfo%myLatEndHalo
2141 numMembers1 = size(weights,1)
2142 numMembers2 = size(weights,2)
2143 totalCount(:) = 0
2144
2145 !$OMP PARALLEL DO PRIVATE(latIndex, lonIndex, interpIndex, interpLatIndex, interpLonIndex, memberIndex1, memberIndex2)
2146 do latIndex = wInterpInfo%myLatBeg, wInterpInfo%myLatEnd
2147 do lonIndex = wInterpInfo%myLonBeg, wInterpInfo%myLonEnd
2148 if (wInterpInfo%numIndexes(lonIndex,latIndex) <= 0) cycle
2149 weights(:,:,lonIndex,latIndex) = 0.0D0
2150 if (wInterpInfo%lonIndexes(1,lonIndex,latIndex) == 0) cycle
2151
2152 totalCount(omp_get_thread_num()+1) = totalCount(omp_get_thread_num()+1) + wInterpInfo%numIndexes(lonIndex,latIndex)
2153
2154 do interpIndex = 1, wInterpInfo%numIndexes(lonIndex,latIndex)
2155 interpLonIndex = wInterpInfo%lonIndexes(interpIndex,lonIndex,latIndex)
2156 interpLatIndex = wInterpInfo%latIndexes(interpIndex,lonIndex,latIndex)
2157
2158 do memberIndex2 = 1, numMembers2
2159 do memberIndex1 = 1, numMembers1
2160 weights(memberIndex1,memberIndex2,lonIndex,latIndex) = &
2161 weights(memberIndex1,memberIndex2,lonIndex,latIndex) + &
2162 wInterpInfo%interpWeights(interpIndex,lonIndex,latIndex) * &
2163 weights(memberIndex1,memberIndex2,interpLonIndex,interpLatIndex)
2164 end do
2165 end do
2166
2167 end do ! interpIndex
2168
2169 end do ! lonIndex
2170 end do ! latIndex
2171 !$OMP END PARALLEL DO
2172
2173 if (firstCall) write(*,*) 'enkf_interpWeights: totalCount = ', totalCount(:)
2174 firstCall = .false.
2175
2176 end subroutine enkf_interpWeights
2177
2178 !--------------------------------------------------------------------------
2179 ! enkf_modifyAMSUBobsError
2180 !--------------------------------------------------------------------------
2181 subroutine enkf_modifyAMSUBobsError(obsSpaceData)
2182 implicit none
2183
2184 ! Arguments:
2185 type(struct_obs), target, intent(inout) :: obsSpaceData
2186
2187 ! Locals:
2188 real(pre_obsReal), parameter :: AMSUB_trop_oer = 1.0 ! assumed value for AMSU-B obs error in tropics
2189 integer :: headerIndex, bodyIndex, bodyIndexBeg, bodyIndexEnd, codeType
2190 real(pre_obsReal) :: lat_obs
2191
2192 ! for AMSUB observations set the observation error std dev equal to 1.0
2193 ! in the larger tropical area where the spread-skill correlation suggests
2194 ! that the data are accurate (.i.e |lat|<40. ). Otherwise don't reduce the
2195 ! observational error.
2196 do headerIndex = 1, obs_numheader(obsSpaceData)
2197 lat_obs = obs_headElem_r(obsSpaceData, obs_lat, headerIndex)
2198 codeType = obs_headElem_i(obsSpaceData, obs_ity, headerIndex)
2199 lat_obs = lat_obs * MPC_DEGREES_PER_RADIAN_R8
2200 if ( abs(lat_obs) < 40. .and. (codeType == codtyp_get_codtyp('amsub') .or. &
2201 codeType == codtyp_get_codtyp('mhs') .or. &
2202 codeType == codtyp_get_codtyp('mwhs2')) ) then
2203 bodyIndexBeg = obs_headElem_i(obsSpaceData, obs_rln, headerIndex)
2204 bodyIndexEnd = obs_headElem_i(obsSpaceData, obs_nlv, headerIndex) + bodyIndexBeg - 1
2205 do bodyIndex = bodyIndexBeg, bodyIndexEnd
2206 call obs_bodySet_r(obsSpaceData, obs_oer, bodyIndex, AMSUB_trop_oer)
2207 end do
2208 end if
2209 end do
2210
2211 end subroutine enkf_modifyAMSUBobsError
2212
2213 !--------------------------------------------------------------------------
2214 ! enkf_rejectHighLatIR
2215 !--------------------------------------------------------------------------
2216 subroutine enkf_rejectHighLatIR(obsSpaceData)
2217 implicit none
2218
2219 ! Arguments:
2220 type(struct_obs), target, intent(inout) :: obsSpaceData
2221
2222 ! Locals:
2223 integer :: headerIndex, bodyIndex, bodyIndexBeg, bodyIndexEnd, codeType
2224 real(pre_obsReal) :: lat_obs
2225
2226 ! reject all HIR radiance observation in arctic and antarctic (.i.e |lat|>60. )
2227 do headerIndex = 1, obs_numheader(obsSpaceData)
2228 lat_obs = obs_headElem_r(obsSpaceData, obs_lat, headerIndex)
2229 codeType = obs_headElem_i(obsSpaceData, obs_ity, headerIndex)
2230 lat_obs = lat_obs * MPC_DEGREES_PER_RADIAN_R8
2231 if ( abs(lat_obs) > 60. .and. tvs_isIdBurpHyperSpectral(codeType) ) then
2232 write(*,*) 'enkf_rejectHighLatIR: !!!!!!!!--------WARNING--------!!!!!!!!'
2233 write(*,*) 'enkf_rejectHighLatIR: This HIR radiance profile was rejected because |lat|>60.'
2234 write(*,*) 'enkf_rejectHighLatIR: latidude= ', lat_obs, 'codtyp= ', codeType
2235 bodyIndexBeg = obs_headElem_i(obsSpaceData, obs_rln, headerIndex)
2236 bodyIndexEnd = obs_headElem_i(obsSpaceData, obs_nlv, headerIndex) + bodyIndexBeg - 1
2237 do bodyIndex = bodyIndexBeg, bodyIndexEnd
2238 call obs_bodySet_i(obsSpaceData, obs_ass, bodyIndex, obs_notAssimilated)
2239 ! also set the 'rejected by selection process' flag (bit 11)
2240 call obs_bodySet_i( obsSpaceData, obs_flg, bodyIndex, &
2241 ibset( obs_bodyElem_i( obsSpaceData, obs_flg, bodyIndex ), 11) )
2242 end do
2243 end if
2244 end do
2245
2246 end subroutine enkf_rejectHighLatIR
2247
2248 !--------------------------------------------------------------------------
2249 ! enkf_getModulatedState
2250 !--------------------------------------------------------------------------
2251 subroutine enkf_getModulatedState( stateVector_in, stateVectorMeanTrl, &
2252 vLocalizeLengthScale, numRetainedEigen, nEns, &
2253 eigenVectorColumnIndex, stateVector_out, &
2254 beSilent )
2255 !
2256 !:Purpose: Compute vertical localization matrix, and the corresponding
2257 ! eigenvectors/eigenvalues, to obtain modulated stateVector.
2258 !
2259 implicit none
2260
2261 ! Arguments:
2262 type(struct_gsv), intent(in) :: stateVector_in
2263 type(struct_gsv), intent(in) :: stateVectorMeanTrl
2264 real(8), intent(in) :: vLocalizeLengthScale
2265 integer, intent(in) :: numRetainedEigen
2266 integer, intent(in) :: nEns
2267 integer, intent(in) :: eigenVectorColumnIndex
2268 type(struct_gsv), intent(inout) :: stateVector_out
2269 logical, intent(in) :: beSilent
2270
2271 ! Locals:
2272 real(4) :: modulationFactor_r4
2273 real(4), pointer :: field_out_r4(:,:,:,:)
2274 integer :: nLev, nlev_out, levIndex, latIndex, lonIndex
2275 integer :: lon1, lon2, lat1, lat2
2276 integer :: varIndex, stepIndex, eigenVectorLevelIndex
2277 character(len=4) :: varName
2278
2279 call utl_tmg_start(130,'--getModulatedState')
2280
2281 if ( .not. beSilent ) write(*,*) 'enkf_getModulatedState: START'
2282
2283 if ( stateVector_in%dataKind /= 4 ) then
2284 call utl_abort('enkf_getModulatedState: only dataKind=4 is implemented')
2285 end if
2286
2287 nLev = stateVector_in%vco%nLev_M
2288 if ( vLocalizeLengthScale <= 0.0d0 .or. nLev <= 1 ) then
2289 call utl_abort('enkf_getModulatedState: no vertical localization')
2290 end if
2291
2292 ! Compute perturbation by subtracting ensMean
2293 call gsv_copy(stateVector_in, stateVector_out, beSilent_opt=beSilent)
2294 call gsv_add(stateVectorMeanTrl, stateVector_out, scaleFactor_opt=-1.0d0)
2295
2296 lon1 = stateVector_out%myLonBeg
2297 lon2 = stateVector_out%myLonEnd
2298 lat1 = stateVector_out%myLatBeg
2299 lat2 = stateVector_out%myLatEnd
2300
2301 ! Compute modulated member perturbation from original member perturbation:
2302 ! v'_k = (Nens*nLamda/(Nens - 1))^1/2 * Lambda^1/2 * E * x'_k
2303 step_loop: do stepIndex = 1, stateVector_out%numStep
2304 var_loop: do varIndex = 1, vnl_numvarmax
2305 varName = vnl_varNameList(varIndex)
2306 if ( .not. gsv_varExist(stateVector_out,varName) ) cycle var_loop
2307
2308 nlev_out = stateVector_out%varNumLev(varIndex)
2309
2310 call gsv_getField(statevector_out,field_out_r4,varName)
2311
2312 do latIndex = lat1, lat2
2313 do lonIndex = lon1, lon2
2314 do levIndex = 1, nlev_out
2315 if ( nlev_out == 1 ) then
2316 eigenVectorLevelIndex = nLev
2317 else
2318 eigenVectorLevelIndex = levIndex
2319 end if
2320
2321 call getModulationFactor( stateVector_in%vco, eigenVectorLevelIndex, &
2322 eigenVectorColumnIndex, numRetainedEigen, &
2323 nEns, vLocalizeLengthScale, &
2324 modulationFactor_r4, beSilent_opt=beSilent )
2325
2326 field_out_r4(lonIndex,latIndex,levIndex,stepIndex) = &
2327 field_out_r4(lonIndex,latIndex,levIndex,stepIndex) * &
2328 modulationFactor_r4
2329 end do
2330 end do
2331 end do
2332
2333 end do var_loop
2334 end do step_loop
2335
2336 ! Now add to ensMean to get modulated member
2337 ! v_k = v'_k + v_mean
2338 call gsv_add(stateVectorMeanTrl, stateVector_out)
2339
2340 if ( .not. beSilent ) write(*,*) 'enkf_getModulatedState: END'
2341
2342 call utl_tmg_stop(130)
2343
2344 end subroutine enkf_getModulatedState
2345
2346 !--------------------------------------------------------------------------
2347 ! enkf_setupModulationFactor
2348 !--------------------------------------------------------------------------
2349 subroutine enkf_setupModulationFactor(vco, numRetainedEigen, nEns, vLocalizeLengthScale, &
2350 beSilent)
2351 !
2352 !:Purpose: setup modulationFactorArray by calling getModulationFactor for first time.
2353 !
2354 implicit none
2355
2356 ! Arguments:
2357 type(struct_vco), pointer, intent(in) :: vco
2358 integer, intent(in) :: numRetainedEigen
2359 integer, intent(in) :: nEns
2360 real(8), intent(in) :: vLocalizeLengthScale
2361 logical, intent(in) :: beSilent
2362
2363 ! Locals:
2364 integer :: eigenVectorColumnIndex
2365 integer :: eigenVectorLevelIndex
2366 real(4) :: modulationFactor_r4
2367
2368 eigenVectorColumnIndex = 1
2369 eigenVectorLevelIndex = 1
2370 call getModulationFactor(vco, eigenVectorLevelIndex, &
2371 eigenVectorColumnIndex, numRetainedEigen, &
2372 nEns, vLocalizeLengthScale, &
2373 modulationFactor_r4, beSilent_opt=beSilent)
2374
2375 end subroutine enkf_setupModulationFactor
2376
2377 !--------------------------------------------------------------------------
2378 ! getModulationFactor
2379 !--------------------------------------------------------------------------
2380 subroutine getModulationFactor( vco, eigenVectorLevelIndex, &
2381 eigenVectorColumnIndex, numRetainedEigen, &
2382 nEns, vLocalizeLengthScale, &
2383 modulationFactor_r4, beSilent_opt )
2384 !
2385 !:Purpose: compute modulation factor needed to multiply ensemble
2386 ! perturbation to get the modulated perturbation:
2387 ! (Nens*nLambda/(Nens - 1))^1/2 * Lambda^1/2
2388 !
2389 implicit none
2390
2391 ! Arguments:
2392 type(struct_vco), pointer, intent(in) :: vco
2393 integer, intent(in) :: eigenVectorLevelIndex
2394 integer, intent(in) :: eigenVectorColumnIndex
2395 integer, intent(in) :: numRetainedEigen
2396 integer, intent(in) :: nEns
2397 real(8), intent(in) :: vLocalizeLengthScale
2398 real(4), intent(out) :: modulationFactor_r4
2399 logical, optional, intent(in) :: beSilent_opt
2400
2401 ! Locals:
2402 integer :: levIndex1, levIndex2, eigenIndex
2403 integer :: nLev, nLev_M, nLev_depth, matrixRank
2404 real(8) :: zr, zcorr, pSurfRef
2405 real(8) :: tolerance
2406 real(8), pointer :: pressureProfile(:)
2407 real(8), allocatable, save :: eigenValues(:)
2408 real(8), allocatable, save :: eigenVectors(:,:)
2409 real(8), allocatable, save :: verticalLocalizationMat(:,:)
2410 real(8), allocatable, save :: verticalLocalizationMatLowRank(:,:)
2411 real(4), allocatable, save :: modulationFactorArray_r4(:,:)
2412 logical :: beSilent
2413
2414 logical, save :: firstCall = .true.
2415
2416 if ( present(beSilent_opt) ) then
2417 beSilent = beSilent_opt
2418 else
2419 beSilent = .false.
2420 end if
2421
2422 ! Compute vertical localization matrix and its eigenValues/Vectors on first call
2423 if ( firstCall ) then
2424 firstCall = .false.
2425 if ( mmpi_myid == 0 .and. .not. beSilent ) then
2426 write(*,*) 'getModulationFactor: computing eigenValues/Vectors'
2427 end if
2428
2429 nLev_M = vco%nLev_M
2430 nLev_depth = vco%nlev_depth
2431 nLev = max(nLev_M,nLev_depth)
2432
2433 allocate(eigenValues(nLev))
2434 allocate(eigenVectors(nLev,nLev))
2435 allocate(verticalLocalizationMat(nLev,nLev))
2436 allocate(verticalLocalizationMatLowRank(nLev,nLev))
2437 allocate(modulationFactorArray_r4(numRetainedEigen,nLev))
2438 verticalLocalizationMatLowRank(:,:) = 0.0d0
2439
2440 pSurfRef = 101000.D0
2441 call czp_fetch1DLevels(vco, pSurfRef, profM_opt=pressureProfile)
2442
2443 call lfn_Setup(LocFunctionWanted='FifthOrder')
2444
2445 ! Calculate 5'th order function
2446 do levIndex1 = 1, nLev
2447 do levIndex2 = 1, nLev
2448 zr = abs(log(pressureProfile(levIndex2)) - log(pressureProfile(levIndex1)))
2449 zcorr = lfn_response(zr,vLocalizeLengthScale)
2450 verticalLocalizationMat(levIndex1,levIndex2) = zcorr
2451 end do
2452 end do
2453
2454 ! Compute eigenValues/Vectors of vertical localization matrix
2455 tolerance = 1.0D-50
2456 call utl_eigenDecomp(verticalLocalizationMat, eigenValues, eigenVectors, &
2457 tolerance, matrixRank)
2458 if ( matrixRank < numRetainedEigen ) then
2459 write(*,*) 'matrixRank=', matrixRank
2460 call utl_abort('getModulationFactor: verticalLocalizationMat is rank deficient=')
2461 end if
2462
2463 ! Compute low-ranked vertical localization matrix
2464 do levIndex1 = 1, nLev
2465 do levIndex2 = 1, nLev
2466 do eigenIndex = 1, numRetainedEigen
2467 verticalLocalizationMatLowRank(levIndex1,levIndex2) = verticalLocalizationMatLowRank(levIndex1,levIndex2) + &
2468 eigenVectors(levIndex1,eigenIndex) * &
2469 eigenVectors(levIndex2,eigenIndex) * &
2470 eigenValues(eigenIndex)
2471 end do
2472 end do
2473 end do
2474
2475 ! now compute the 2D modulationFactor array
2476 do levIndex1 = 1, nLev
2477 do eigenIndex = 1, numRetainedEigen
2478 modulationFactorArray_r4(eigenIndex,levIndex1) = real( &
2479 1 / sqrt(verticalLocalizationMatLowRank(levIndex1,levIndex1)) * &
2480 eigenVectors(levIndex1,eigenIndex) * &
2481 eigenValues(eigenIndex) ** 0.5 * &
2482 (nEns * numRetainedEigen / (nEns - 1)) ** 0.5,4)
2483 end do
2484 end do
2485
2486 if ( mmpi_myid == 0 .and. .not. beSilent ) then
2487 do levIndex1 = 1, numRetainedEigen
2488 write(*,*) 'getModulationFactor: eigen mode=', levIndex1, ', eigenVectors=', eigenVectors(:,levIndex1)
2489 end do
2490 write(*,*) 'getModulationFactor: eigenValues=', eigenValues(1:numRetainedEigen)
2491
2492 do levIndex1 = 1, nLev
2493 write(*,*) 'getModulationFactor: verticalLocalizationMat for lev ', levIndex1, '=', verticalLocalizationMat(levIndex1,:)
2494 write(*,*) 'getModulationFactor: verticalLocalizationMatLowRank for lev ', levIndex1, '=', verticalLocalizationMatLowRank(levIndex1,:)
2495 end do
2496 end if
2497 end if
2498
2499 modulationFactor_r4 = modulationFactorArray_r4(eigenVectorColumnIndex,eigenVectorLevelIndex)
2500
2501 end subroutine getModulationFactor
2502
2503end module enkf_mod
