1module bMatrixEnsemble_mod
2 ! MODULE bMatrixEnsemble_mod (prefix='ben' category='2. B and R matrices')
3 !
4 !:Purpose: Performs transformation from control vector to analysis increment
5 ! (and the adjoint transformation) using the spatially localized
6 ! ensemble covariance matrix. This module works for both global and
7 ! limited-area applications.
8 !
9 use midasMpi_mod
10 use fileNames_mod
11 use gridStateVector_mod
12 use gridStateVectorFileIO_mod
13 use ensembleStateVector_mod
14 use horizontalCoord_mod
15 use verticalCoord_mod
16 use timeCoord_mod
17 use localization_mod
18 use mathPhysConstants_mod
19 use gridVariableTransforms_mod
20 use utilities_mod
21 use globalSpectralTransform_mod
22 use lamSpectralTransform_mod
23 use spectralFilter_mod
24 use varNameList_mod
25 use advection_mod
26 use gridBinning_mod
27 use humidityLimits_mod
28 use lamAnalysisGridTransforms_mod
29 use calcHeightAndPressure_mod
30 implicit none
31 save
32 private
33
34 ! public procedures
35 public :: ben_Setup, ben_BSqrt, ben_BSqrtAd, ben_writeAmplitude
36 public :: ben_reduceToMPILocal, ben_reduceToMPILocal_r4, ben_expandToMPIGlobal, ben_expandToMPIGlobal_r4
37 public :: ben_getScaleFactor, ben_getnEns, ben_getPerturbation, ben_getEnsMean, ben_Finalize
38 public :: ben_setFsoLeadTime, ben_getNumStepAmplitudeAssimWindow, ben_getAmplitudeAssimWindow
39 public :: ben_getAmp3dStepIndexAssimWindow, ben_getNumLoc, ben_getLoc, ben_getNumInstance
40
41 integer, parameter :: maxNumLocalLength = 20
42 integer, parameter :: nInstanceMax = 10
43
44 type :: struct_bEns
45 logical :: initialized = .false.
46
47 real(8),allocatable :: scaleFactor_M(:), scaleFactor_T(:), scaleFactor_DP(:)
48 real(8) :: scaleFactor_SF
49
50 integer :: ni, nj, myLonBeg, myLonEnd, myLatBeg, myLatEnd
51 integer :: nLevInc_M, nLevInc_T, nLevInc_DP, nLevEns_M, nLevEns_T, nLevEns_DP
52 integer :: topLevIndex_M, topLevIndex_T, topLevIndex_DP
53 integer :: nEnsOverDimension
54 integer :: cvDim_mpilocal, cvDim_mpiglobal
55 integer :: numStep, numStepAssimWindow
56 integer :: numStepAmplitudeFSOFcst, numStepAmplitudeAssimWindow, numStepAdvectAssimWindow
57 integer :: numSubEns
58 integer,allocatable :: dateStampList(:)
59 integer,allocatable :: dateStampListAdvectedFields(:)
60
61 integer :: numIncludeAnlVar
62
63 ! FSO
64 real(8) :: fsoLeadTime = -1.0D0
65 integer :: numStepAdvectFSOFcst
66
67 ! Localizations
68 integer :: nWaveBand
69 integer :: nWaveBandForFiltering = 0
70
71 ! Ensemble perturbations
72 type(struct_ens), allocatable :: ensPerts(:)
73
74 ! Ensemble amplitude (only used in diagnostic mode)
75 type(struct_ens) :: ensAmplitudeStorage
76 character(len=4) :: varNameALFA(1)
77
78 ! Localization
79 type(struct_loc), pointer :: locStorage(:)
80
81 ! The HU LQ mess
82 logical :: gsvHUcontainsLQ
83 logical :: ensShouldNotContainLQvarName
84
85 ! Vertical grid
86 type(struct_vco), pointer :: vco_anl, vco_ens, vco_file => null()
87
88 ! Horizontal grid
89 type(struct_hco), pointer :: hco_ens ! Ensemble horizontal grid parameters
90 type(struct_hco), pointer :: hco_core ! Core grid for limited area EnVar
91 type(struct_hco), pointer :: hco_file ! Input file horizontal grid parameters
92
93 ! Amplitude parameters
94 type(struct_adv) :: adv_amplitudeFSOFcst
95 type(struct_adv), pointer :: adv_amplitudeAssimWindow
96 type(struct_adv) :: adv_ensPerts
97 type(struct_adv) :: adv_analInc
98
99 integer :: amp3dStepIndexAssimWindow
100 integer :: amp3dStepIndexFSOFcst
101
102 ! Variance smoothing
103 logical :: ensPertsNormalized
104
105 type(struct_gsv) :: statevector_ensStdDev
106
107 ! Optimization
108 logical :: useSaveAmp
109
110 ! Copy of namelist variables
111 integer :: nEns
112 real(8) :: scaleFactor(vco_maxNumLevels)
113 real(8) :: scaleFactorHumidity(vco_maxNumLevels)
114 real(8) :: advectFactorFSOFcst(vco_maxNumLevels)
115 real(8) :: advectFactorAssimWindow(vco_maxNumLevels)
116 integer :: ntrunc
117 character(len=256) :: enspathname
118 real(8) :: hLocalize(maxNumLocalLength)
119 real(8) :: vLocalize(maxNumLocalLength)
120 character(len=256) :: LocalizationType
121 integer :: waveBandPeaks(maxNumLocalLength)
122 integer :: waveBandIndexSelected
123 logical :: ensDiagnostic
124 logical :: advDiagnostic
125 character(len=2) :: ctrlVarHumidity
126 character(len=32) :: advectTypeAssimWindow
127 character(len=32) :: advectStartTimeIndexAssimWindow
128 logical :: advectAmplitudeFSOFcst
129 logical :: advectAmplitudeAssimWindow = .false.
130 logical :: advectEnsPertAnlInc = .false.
131 logical :: removeSubEnsMeans
132 logical :: keepAmplitude
133 character(len=4) :: IncludeAnlVar(vnl_numvarmax)
134 logical :: ensContainsFullField
135 character(len=24) :: varianceSmoothing
136 real(8) :: footprintRadius
137 real(8) :: footprintTopoThreshold
138 logical :: useCmatrixOnly
139 integer :: ensDateOfValidity
140 character(len=20) :: transformVarKindCH
141 real(8) :: huMinValue
142 character(len=12) :: hInterpolationDegree
143 end type struct_bEns
144
145 integer :: nInstance = 0 ! The number of Bens instances
146
147 type(struct_bEns) :: bEns(nInstanceMax)
148
149 type(struct_ens), target :: ensAmplitudeSave(nInstanceMax) ! Save this to allow early allocation
150 ! for better efficiency
151
152 character(len=15) :: ben_mode
153
154 character(len=4), parameter :: varNameALFAatmMM(1) = (/ 'ALFA' /)
155 character(len=4), parameter :: varNameALFAatmTH(1) = (/ 'ALFT' /)
156 character(len=4), parameter :: varNameALFAsfc(1) = (/ 'ALFS' /)
157 character(len=4), parameter :: varNameALFAocean(1) = (/ 'ALFO' /)
158
159 logical, parameter :: verbose = .false. ! Control parameter for the level of listing output
160
161 integer, external :: get_max_rss, omp_get_thread_num
162
163CONTAINS
164
165 !--------------------------------------------------------------------------
166 ! ben_setup
167 !--------------------------------------------------------------------------
168 subroutine ben_setup(hco_anl_in, hco_core_in, vco_anl_in, cvDimPerInstance, &
169 mode_opt)
170 !
171 !:Purpose: To configure the ensemble B matrix
172 !
173 implicit none
174
175 ! Arguments:
176 type(struct_hco), pointer, intent(in) :: hco_anl_in
177 type(struct_hco), pointer, intent(in) :: hco_core_in
178 type(struct_vco), pointer, intent(in) :: vco_anl_in
179 character(len=*), optional, intent(in) :: mode_opt
180 integer, allocatable, intent(out) :: cvDimPerInstance(:)
181
182 ! Locals:
183 integer :: fnom, fclos, ierr
184 integer :: cvDimStorage(nInstanceMax)
185 integer :: nulnam = 0
186 ! Namelist variables
187 integer :: nEns ! number of ensemble members
188 real(8) :: scaleFactor(vco_maxNumLevels) ! level-dependent scaling of variances for all variables
189 real(8) :: scaleFactorHumidity(vco_maxNumLevels) ! level-dependent scaling of variances for humidity
190 real(8) :: advectFactorFSOFcst(vco_maxNumLevels) ! level-dependent scaling of winds used to advect localization
191 real(8) :: advectFactorAssimWindow(vco_maxNumLevels) ! level-dependent scaling of winds used to advect localization
192 integer :: ntrunc ! spectral truncation used for horizontal localization function
193 character(len=256) :: enspathname ! path where ensemble members are located (usually ./ensemble)
194 real(8) :: hLocalize(maxNumLocalLength) ! horiz. localization length scale for each waveband (in km)
195 real(8) :: vLocalize(maxNumLocalLength) ! vert. localization length scale for each waveband (in scale heights)
196 character(len=256) :: localizationType ! "LevelDependent", "ScaleDependent" or "ScaleDependentWithSpectralLoc"
197 integer :: waveBandPeaks(maxNumLocalLength) ! total wavenumber corresponding to peak of each waveband for SDL
198 integer :: waveBandIndexSelected ! for multiple NAMBEN blocks, waveband index of this block
199 logical :: ensDiagnostic ! when `.true.` write diagnostic info related to ens. to files
200 logical :: advDiagnostic ! when `.true.` write diagnostic info related to advection to files
201 character(len=2) :: ctrlVarHumidity ! name of humidity variable used for ensemble perturbations (LQ or HU)
202 character(len=32) :: advectTypeAssimWindow ! what is advected in assim. window: "amplitude" or "ensPertAnlInc"
203 character(len=32) :: advectStartTimeIndexAssimWindow ! time index where advection originates from "first" or "middle"
204 logical :: advectAmplitudeFSOFcst = .false. ! activate advection of ens. member amplitudes for FSOI calculation
205 logical :: advectAmplitudeAssimWindow = .false. ! activate advection of ens. member amplitudes for 4D-EnVar
206 logical :: advectEnsPertAnlInc = .false. ! activate advection of ens. pert. and anl. incr. for 4D-EnVar
207 logical :: removeSubEnsMeans ! remove mean of each subsensemble defined by "subEnsembleIndex.txt"
208 logical :: keepAmplitude ! activate storage of ens. amplitudes in instance of struct_ens
209 character(len=4) :: includeAnlVar(vnl_numvarmax) ! list of state variables for this ensemble B matrix; use all if blank
210 logical :: ensContainsFullField ! indicates full fields and not perturbations are in the ens. files
211 character(len=24) :: varianceSmoothing ! "none", "horizMean", "footprint" or "footprintLandSeaTopo"
212 real(8) :: footprintRadius ! parameter for variance smoothing (in meters)
213 real(8) :: footprintTopoThreshold ! parameter for variance smoothing (in meters)
214 logical :: useCmatrixOnly ! activate normalization of ens. perturbations by ens. stddev
215 integer :: ensDateOfValidity ! when set to -1, date in ens. files is ignored (only for 3D ens.)
216 real(8) :: huMinValue ! minimum humidity value imposed on ensemble members
217 character(len=12) :: hInterpolationDegree ! select degree of horizontal interpolation (if needed)
218 character(len=20) :: transformVarKindCH ! name of transform performed on chemistry-related variables in ens.
219
220 ! Namelist
221 NAMELIST /NAMBEN/nEns, scaleFactor, scaleFactorHumidity, ntrunc, enspathname, &
222 hLocalize, vLocalize, LocalizationType, waveBandPeaks, ensDiagnostic, advDiagnostic, &
223 ctrlVarHumidity, advectFactorFSOFcst, advectFactorAssimWindow, removeSubEnsMeans, &
224 keepAmplitude, advectTypeAssimWindow, advectStartTimeIndexAssimWindow, IncludeAnlVar,&
225 ensContainsFullField, varianceSmoothing, footprintRadius, footprintTopoThreshold, &
226 useCmatrixOnly, waveBandIndexSelected, ensDateOfValidity, transformVarKindCH, &
227 huMinValue, hInterpolationDegree
228
229 if (verbose) write(*,*) 'Entering ben_Setup'
230
231 call utl_tmg_start(54,'----B_ENS_Setup')
232
233 !- Set the module mode
234 if ( present(mode_opt) ) then
235 if ( trim(mode_opt) == 'Analysis' .or. trim(mode_opt) == 'BackgroundCheck') then
236 ben_mode = trim(mode_opt)
237 if (mmpi_myid == 0) write(*,*)
238 if (mmpi_myid == 0) write(*,*) 'ben_setup: Mode activated = ', trim(ben_mode)
239 else
240 write(*,*)
241 write(*,*) 'mode = ', trim(mode_opt)
242 call utl_abort('ben_setup: unknown mode')
243 end if
244 else
245 ben_mode = 'Analysis'
246 if (mmpi_myid == 0) write(*,*)
247 if (mmpi_myid == 0) write(*,*) 'ben_setup: Analysis mode activated (by default)'
248 end if
249
250 !- Open the namelist and loop through it
251 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
252
253 instanceLoop: do
254
255 !- Set the default values for the namalist parameters
256 scaleFactor(:) = 0.0d0
257 scaleFactorHumidity(:)= 1.0d0
258 nEns = 10
259 ntrunc = 30
260 enspathname = 'ensemble'
261 localizationType = 'LevelDependent'
262 waveBandPeaks(:) = -1.0d0
263 waveBandIndexSelected = -1
264 ensDiagnostic = .false.
265 advDiagnostic = .false.
266 hLocalize(:) = -1.0d0
267 hLocalize(1) = 2800.0d0
268 vLocalize(:) = -1.0d0
269 vLocalize(1) = 2.0d0
270 ctrlVarHumidity = 'LQ'
271 advectFactorFSOFcst(:)= 0.0D0
272 advectTypeAssimWindow = 'amplitude'
273 advectStartTimeIndexAssimWindow = 'first'
274 advectFactorAssimWindow(:) = 0.0D0
275 removeSubEnsMeans = .false.
276 keepAmplitude = .false.
277 ensContainsFullField = .true.
278 includeAnlVar(:) = ''
279 varianceSmoothing = 'none'
280 footprintRadius = 250.0d3 ! 250km
281 footprintTopoThreshold = 200.0d0 ! 200 m
282 useCmatrixOnly = .false.
283 ensDateOfValidity = MPC_missingValue_INT ! i.e. undefined
284 transformVarKindCH = ''
285 huMinValue = MPC_missingValue_R8 ! i.e. undefined
286 hInterpolationDegree = 'LINEAR' ! or 'CUBIC' or 'NEAREST'
287
288 !- Read the namelist
289 read(nulnam,nml=namben,iostat=ierr)
290 if (ierr /= 0) then
291 if (nInstance >= 1) then
292 exit instanceLoop
293 else
294 call utl_abort('ben_setup: Error reading the first instance namelist')
295 end if
296 end if
297
298 if (mmpi_myid == 0) write(*,nml=namben)
299
300 ! We have found a valid instance
301 nInstance = nInstance + 1
302
303 !- Adjust some namelist-dependent variables
304
305 ! If zero weight, skip rest of setup
306 if ( sum(scaleFactor(:)) == 0.0d0 ) then
307 if (mmpi_myid == 0) write(*,*) 'ben_setup: scaleFactor=0, skipping rest of setup and exit instance loop'
308 cvDimStorage(nInstance) = 0
309 bEns(nInstance)%initialized = .false.
310 exit instanceLoop
311 end if
312
313 if (nInstance > nInstanceMax) then
314 call utl_abort('ben_setup: the number of instance exceed the maximum currently allowed')
315 end if
316
317 if (trim(ben_mode) == 'BackgroundCheck' .and. nInstance > 1) then
318 call utl_abort('ben_setup: the background check mode is not compatible with multiple instance')
319 end if
320
321 if ( (huMinValue == MPC_missingValue_R8) .and. &
322 gsv_varExist(varName='HU') .and. &
323 (ctrlVarHumidity == 'LQ') ) then
324 call utl_abort('ben_setup: the value of huMinValue must be specified in namelist NAMBEN')
325 end if
326
327 !- Transfer the info to the structure
328 bEns(nInstance)%nEns = nEns
329 bEns(nInstance)%scaleFactor(:) = scaleFactor(:)
330 bEns(nInstance)%scaleFactorHumidity(:) = scaleFactorHumidity(:)
331 bEns(nInstance)%nTrunc = nTrunc
332 bEns(nInstance)%ensPathName = ensPathName
333 bEns(nInstance)%hLocalize(:) = hLocalize(:)
334 bEns(nInstance)%vLocalize(:) = vLocalize(:)
335 bEns(nInstance)%localizationType = localizationType
336 bEns(nInstance)%waveBandPeaks(:) = waveBandPeaks(:)
337 bEns(nInstance)%waveBandIndexSelected = waveBandIndexSelected
338 bEns(nInstance)%ensDiagnostic = ensDiagnostic
339 bEns(nInstance)%advDiagnostic = advDiagnostic
340 bEns(nInstance)%ctrlVarHumidity = ctrlVarHumidity
341 bEns(nInstance)%advectTypeAssimWindow = advectTypeAssimWindow
342 bEns(nInstance)%advectStartTimeIndexAssimWindow = advectStartTimeIndexAssimWindow
343 bEns(nInstance)%advectFactorFSOFcst(:) = advectFactorFSOFcst(:)
344 bEns(nInstance)%advectFactorAssimWindow(:) = advectFactorAssimWindow(:)
345 bEns(nInstance)%advectAmplitudeFSOFcst = advectAmplitudeFSOFcst
346 bEns(nInstance)%advectAmplitudeAssimWindow = advectAmplitudeAssimWindow
347 bEns(nInstance)%advectEnsPertAnlInc = advectEnsPertAnlInc
348 bEns(nInstance)%removeSubEnsMeans = removeSubEnsMeans
349 bEns(nInstance)%keepAmplitude = keepAmplitude
350 bEns(nInstance)%includeAnlVar(:) = includeAnlVar(:)
351 bEns(nInstance)%ensContainsFullField = ensContainsFullField
352 bEns(nInstance)%varianceSmoothing = varianceSmoothing
353 bEns(nInstance)%footprintRadius = footprintRadius
354 bEns(nInstance)%footprintTopoThreshold = footprintTopoThreshold
355 bEns(nInstance)%useCmatrixOnly = useCmatrixOnly
356 bEns(nInstance)%ensDateOfValidity = ensDateOfValidity
357 bEns(nInstance)%transformVarKindCH = transformVarKindCH
358 bEns(nInstance)%huMinValue = huMinValue
359 bEns(nInstance)%hInterpolationDegree = hInterpolationDegree
360
361 bEns(nInstance)%hco_ens => hco_anl_in
362 bEns(nInstance)%hco_core => hco_core_in
363 bEns(nInstance)%vco_anl => vco_anl_in
364
365 !- Setup the LAM analysis grid metrics
366 call lgt_setupFromHCO( hco_anl_in, hco_core_in ) ! IN
367
368 !- Set the instance
369 call ben_setupOneInstance(nInstance, & ! IN
370 cvDimStorage(nInstance)) ! OUT
371
372 bEns(nInstance)%initialized = .true.
373
374 end do instanceLoop
375
376 !- Close the namelist
377 ierr = fclos(nulnam)
378
379 !- Set the output control variable dimensions array
380 allocate(cvDimPerInstance(nInstance))
381 cvDimPerInstance(:) = cvDimStorage(1:nInstance)
382
383 call utl_tmg_stop(54)
384
385 end subroutine ben_setup
386
387 !--------------------------------------------------------------------------
388 ! ben_setupOneInstance
389 !--------------------------------------------------------------------------
390 subroutine ben_setupOneInstance(instanceIndex, cvDim)
391 !
392 !:Purpose: To configure a single instance of the ensemble B matrix
393 !
394 implicit none
395
396 ! Arguments:
397 integer, intent(in) :: instanceIndex
398 integer, intent(out) :: cvDim
399
400 ! Locals:
401 type(struct_gsv) :: statevector_ensMean4D, statevector_oneEnsPert4D
402 type(struct_gbi) :: gbi_horizontalMean, gbi_landSeaTopo
403 real(8) :: pSurfRef, delT_hour
404 real(8), allocatable :: advectFactorFSOFcst_M(:),advectFactorAssimWindow_M(:)
405 real(8),pointer :: vertLocationEns(:), vertLocationFile(:), vertLocationInc(:)
406 real(4), pointer :: bin2d(:,:,:)
407 real(8), pointer :: HeightSfc(:,:)
408 integer :: lonPerPE, latPerPE, lonPerPEmax, latPerPEmax
409 integer :: myMemberBeg, myMemberEnd, myMemberCount, maxMyMemberCount
410 integer :: levIndex, jvar, ierr
411 integer :: waveBandIndex, stepIndex
412 character(len=256) :: ensFileName
413 integer :: dateStampFSO, ensDateStampOfValidity, idate, itime, newdate
414 logical :: EnsTopMatchesAnlTop, useAnlLevelsOnly
415 character(len=32) :: direction, directionEnsPerts, directionAnlInc
416
417 if (verbose) write(*,*) 'Entering ben_SetupOneInstance'
418
419 write(*,*) 'ben_setupOneInstance: enspathname = ', trim(bEns(instanceIndex)%ensPathName)
420
421 !
422 !- 1. B matrix configuration
423 !
424
425 !- 1.1 Number of time step bins
426 bEns(instanceIndex)%numStep = tim_nstepobsinc
427 if (bEns(instanceIndex)%numStep /= 1.and.bEns(instanceIndex)%numStep /= 3.and.bEns(instanceIndex)%numStep /= 5.and.bEns(instanceIndex)%numStep /= 7) then
428 call utl_abort('ben_setupOneInstance: Invalid value for numStep (choose 1 or 3 or 5 or 7)!')
429 end if
430
431 !- 1.2 FSO-related options
432 bEns(instanceIndex)%numStepAssimWindow = bEns(instanceIndex)%numStep
433 if (bEns(instanceIndex)%fsoLeadTime > 0.0D0) then
434 bEns(instanceIndex)%numStep = bEns(instanceIndex)%numStep + 1
435 call incdatr(dateStampFSO, tim_getDatestamp(), bEns(instanceIndex)%fsoLeadTime)
436 end if
437
438 allocate(bEns(instanceIndex)%dateStampList(bEns(instanceIndex)%numStep))
439 if (bEns(instanceIndex)%fsoLeadTime > 0.0D0) then
440 call tim_getstamplist(bEns(instanceIndex)%dateStampList,bEns(instanceIndex)%numStep-1,tim_getDatestamp())
441 bEns(instanceIndex)%dateStampList(bEns(instanceIndex)%numStep) = dateStampFSO
442 else
443 if (bEns(instanceIndex)%ensDateOfValidity == MPC_missingValue_INT) then
444 call tim_getstamplist(bEns(instanceIndex)%dateStampList,bEns(instanceIndex)%numStep,tim_getDatestamp())
445 else
446 if (bEns(instanceIndex)%numStep == 1) then
447 if (bEns(instanceIndex)%ensDateOfValidity == -1) then
448 ensDateStampOfValidity = bEns(instanceIndex)%ensDateOfValidity
449 else
450 idate = bEns(instanceIndex)%ensDateOfValidity/100
451 itime = (bEns(instanceIndex)%ensDateOfValidity-idate*100)*1000000
452 ierr = newdate(ensDateStampOfValidity, idate, itime, 3)
453 end if
454 bEns(instanceIndex)%dateStampList(:) = ensDateStampOfValidity
455 else
456 call utl_abort('ben_setupOneInstance: A single date of validity cannot be specified for numStep > 1')
457 end if
458 end if
459 end if
460
461 !- 1.3 Horizontal grid
462 bEns(instanceIndex)%ni = bEns(instanceIndex)%hco_ens%ni
463 bEns(instanceIndex)%nj = bEns(instanceIndex)%hco_ens%nj
464 if (bEns(instanceIndex)%hco_ens%global) then
465 if (mmpi_myid == 0) write(*,*)
466 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: GLOBAL mode activated'
467 else
468 if (mmpi_myid == 0) write(*,*)
469 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: LAM mode activated'
470 end if
471
472 !- 1.4 Vertical levels
473 if ( mmpi_myid == 0 ) then
474 call fln_ensfileName(ensFileName, bEns(instanceIndex)%ensPathName, memberIndex_opt=1)
475 call vco_SetupFromFile(bEns(instanceIndex)%vco_file, ensFileName)
476 end if
477 call vco_mpiBcast(bEns(instanceIndex)%vco_file)
478
479 !- Do we need to read all the vertical levels from the ensemble?
480 useAnlLevelsOnly = vco_subsetOrNot(bEns(instanceIndex)%vco_anl, bEns(instanceIndex)%vco_file)
481 if ( useAnlLevelsOnly ) then
482 write(*,*)
483 write(*,*) 'ben_setupOneInstance: only the analysis levels will be read in the ensemble '
484 bEns(instanceIndex)%vco_ens => bEns(instanceIndex)%vco_anl ! the ensemble target grid is the analysis grid
485 call vco_deallocate(bEns(instanceIndex)%vco_file)
486 bEns(instanceIndex)%vco_file => bEns(instanceIndex)%vco_anl ! only the analysis levels will be read in the ensemble
487 EnsTopMatchesAnlTop = .true.
488 else
489 write(*,*)
490 write(*,*) 'ben_setupOneInstance: all the vertical levels will be read in the ensemble '
491 if ( bEns(instanceIndex)%vco_anl%nLev_M > 0 .and. bEns(instanceIndex)%vco_anl%vgridPresent ) then
492 pSurfRef = 101000.D0
493 call czp_fetch1DLevels(bEns(instanceIndex)%vco_anl, pSurfRef, &
494 profM_opt=vertLocationInc)
495 call czp_fetch1DLevels(bEns(instanceIndex)%vco_file, pSurfRef, &
496 profM_opt=vertLocationFile)
497
498 do levIndex = 1, bEns(instanceIndex)%vco_anl%nLev_M
499 vertLocationInc(levIndex) = log(vertLocationInc(levIndex))
500 end do
501 do levIndex = 1, bEns(instanceIndex)%vco_file%nLev_M
502 vertLocationFile(levIndex) = log(vertLocationFile(levIndex))
503 end do
504
505 EnsTopMatchesAnlTop = abs( vertLocationFile(1) - vertLocationInc(1) ) < 0.1d0
506 write(*,*) 'ben_setupOneInstance: EnsTopMatchesAnlTop, presEns, presInc = ', &
507 EnsTopMatchesAnlTop, vertLocationFile(1), vertLocationInc(1)
508 deallocate(vertLocationFile)
509 deallocate(vertLocationInc)
510 else
511 ! not sure what this mean when no MM levels
512 write(*,*) 'ben_setupOneInstance: nLev_M = ', bEns(instanceIndex)%vco_anl%nLev_M
513 write(*,*) 'ben_setupOneInstance: vgridPresent = ', bEns(instanceIndex)%vco_anl%vgridPresent
514 EnsTopMatchesAnlTop = .true.
515 end if
516
517 if ( EnsTopMatchesAnlTop ) then
518 if ( mmpi_myid == 0 ) write(*,*) 'ben_setupOneInstance: top level of ensemble member and analysis grid match'
519 bEns(instanceIndex)%vco_ens => bEns(instanceIndex)%vco_anl ! IMPORTANT: top levels DO match, therefore safe
520 ! to force members to be on analysis vertical levels
521 else
522 if ( mmpi_myid == 0 ) write(*,*) 'ben_setupOneInstance: top level of ensemble member and analysis grid are different, therefore'
523 if ( mmpi_myid == 0 ) write(*,*) ' assume member is already be on correct levels - NO CHECKING IS DONE'
524 bEns(instanceIndex)%vco_ens => bEns(instanceIndex)%vco_file ! IMPORTANT: top levels do not match, therefore must
525 ! assume file is already on correct vertical levels
526 end if
527 end if
528
529 if (bEns(instanceIndex)%vco_anl%Vcode /= bEns(instanceIndex)%vco_ens%Vcode) then
530 write(*,*) 'ben_setupOneInstance: vco_anl%Vcode = ', bEns(instanceIndex)%vco_anl%Vcode, ', vco_ens%Vcode = ', bEns(instanceIndex)%vco_ens%Vcode
531 call utl_abort('ben_setupOneInstance: vertical levels of ensemble not compatible with analysis grid')
532 end if
533 bEns(instanceIndex)%nLevEns_M = bEns(instanceIndex)%vco_ens%nLev_M
534 bEns(instanceIndex)%nLevEns_T = bEns(instanceIndex)%vco_ens%nLev_T
535 bEns(instanceIndex)%nLevEns_DP = bEns(instanceIndex)%vco_ens%nLev_Depth
536 bEns(instanceIndex)%nLevInc_M = bEns(instanceIndex)%vco_anl%nLev_M
537 bEns(instanceIndex)%nLevInc_T = bEns(instanceIndex)%vco_anl%nLev_T
538 bEns(instanceIndex)%nLevInc_DP = bEns(instanceIndex)%vco_anl%nLev_Depth
539 bEns(instanceIndex)%topLevIndex_M = bEns(instanceIndex)%nLevInc_M - bEns(instanceIndex)%nLevEns_M+1
540 bEns(instanceIndex)%topLevIndex_T = bEns(instanceIndex)%nLevInc_T - bEns(instanceIndex)%nLevEns_T+1
541 bEns(instanceIndex)%topLevIndex_DP = bEns(instanceIndex)%nLevInc_DP - bEns(instanceIndex)%nLevEns_DP+1
542
543 if (bEns(instanceIndex)%vco_anl%Vcode == 5002) then
544 if ( (bEns(instanceIndex)%nLevEns_T /= (bEns(instanceIndex)%nLevEns_M+1)) .and. (bEns(instanceIndex)%nLevEns_T /= 1 .or. bEns(instanceIndex)%nLevEns_M /= 1) ) then
545 write(*,*) 'ben_setupOneInstance: nLevEns_T, nLevEns_M = ',bEns(instanceIndex)%nLevEns_T,bEns(instanceIndex)%nLevEns_M
546 call utl_abort('ben_setupOneInstance: Vcode=5002, nLevEns_T must equal nLevEns_M+1!')
547 end if
548 else if (bEns(instanceIndex)%vco_anl%Vcode == 5005) then
549 if ( bEns(instanceIndex)%nLevEns_T /= bEns(instanceIndex)%nLevEns_M .and. &
550 bEns(instanceIndex)%nLevEns_T /= 0 .and. &
551 bEns(instanceIndex)%nLevEns_M /= 0 ) then
552 write(*,*) 'ben_setup: nLevEns_T, nLevEns_M = ',bEns(instanceIndex)%nLevEns_T,bEns(instanceIndex)%nLevEns_M
553 call utl_abort('ben_setupOneInstance: Vcode=5005, nLevEns_T must equal nLevEns_M!')
554 end if
555 else if (bEns(instanceIndex)%vco_anl%Vcode == 0) then
556 if ( bEns(instanceIndex)%nLevEns_T /= 0 .and. bEns(instanceIndex)%nLevEns_M /= 0 ) then
557 write(*,*) 'ben_setup: nLevEns_T, nLevEns_M = ',bEns(instanceIndex)%nLevEns_T, bEns(instanceIndex)%nLevEns_M
558 call utl_abort('ben_setupOneInstance: surface-only case (Vcode=0), bEns(instanceIndex)%nLevEns_T and nLevEns_M must equal 0!')
559 end if
560 else
561 write(*,*) 'vco_anl%Vcode = ',bEns(instanceIndex)%vco_anl%Vcode
562 call utl_abort('ben_setupOneInstance: unknown vertical coordinate type!')
563 end if
564
565 if (bEns(instanceIndex)%nLevEns_M.gt.bEns(instanceIndex)%nLevInc_M) then
566 call utl_abort('ben_setupOneInstance: ensemble has more levels than increment - not allowed!')
567 end if
568
569 if (bEns(instanceIndex)%nLevEns_M.lt.bEns(instanceIndex)%nLevInc_M) then
570 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: ensemble has less levels than increment'
571 if (mmpi_myid == 0) write(*,*) ' some levels near top will have zero increment'
572 end if
573
574 !- 1.5 Bmatrix Weight
575 if (bEns(instanceIndex)%vco_anl%Vcode == 5002 .or. bEns(instanceIndex)%vco_anl%Vcode == 5005) then
576 if (bEns(instanceIndex)%nLevEns_M > 0) then
577 ! Multi-level or momentum-level-only analysis
578 bEns(instanceIndex)%varNameALFA(:) = varNameALFAatmMM(:)
579 else
580 ! Thermo-level-only analysis
581 bEns(instanceIndex)%varNameALFA(:) = varNameALFAatmTH(:)
582 end if
583 allocate(bEns(instanceIndex)%scaleFactor_M(bEns(instanceIndex)%nLevEns_M))
584 allocate(bEns(instanceIndex)%scaleFactor_T(bEns(instanceIndex)%nLevEns_T))
585 do levIndex = 1, bEns(instanceIndex)%nLevEns_T
586 if (bEns(instanceIndex)%scaleFactor(levIndex) > 0.0d0) then
587 bEns(instanceIndex)%scaleFactor(levIndex) = sqrt(bEns(instanceIndex)%scaleFactor(levIndex))
588 else
589 bEns(instanceIndex)%scaleFactor(levIndex) = 0.0d0
590 end if
591 end do
592 bEns(instanceIndex)%scaleFactor_T(1:bEns(instanceIndex)%nLevEns_T) = bEns(instanceIndex)%scaleFactor(1:bEns(instanceIndex)%nLevEns_T)
593 if (bEns(instanceIndex)%vco_anl%Vcode == 5002) then
594 bEns(instanceIndex)%scaleFactor_M(1:bEns(instanceIndex)%nLevEns_M) = bEns(instanceIndex)%scaleFactor(2:(bEns(instanceIndex)%nLevEns_M+1))
595 else
596 bEns(instanceIndex)%scaleFactor_M(1:bEns(instanceIndex)%nLevEns_M) = bEns(instanceIndex)%scaleFactor(1:bEns(instanceIndex)%nLevEns_M)
597 end if
598
599 do levIndex = 1, bEns(instanceIndex)%nLevEns_T
600 if (bEns(instanceIndex)%scaleFactorHumidity(levIndex) > 0.0d0) then
601 bEns(instanceIndex)%scaleFactorHumidity(levIndex) = sqrt(bEns(instanceIndex)%scaleFactorHumidity(levIndex))
602 else
603 bEns(instanceIndex)%scaleFactorHumidity(levIndex) = 0.0d0
604 end if
605 end do
606
607 bEns(instanceIndex)%scaleFactor_SF = bEns(instanceIndex)%scaleFactor_T(bEns(instanceIndex)%nLevEns_T)
608
609 else if (bEns(instanceIndex)%nLevEns_DP > 0) then
610 ! Ocean variables on depth levels
611 write(*,*) 'nlev_Depth=', bEns(instanceIndex)%nLevEns_DP
612 bEns(instanceIndex)%varNameALFA(:) = varNameALFAocean(:)
613 allocate(bEns(instanceIndex)%scaleFactor_DP(bEns(instanceIndex)%nLevEns_DP))
614 do levIndex = 1, bEns(instanceIndex)%nLevEns_DP
615 if (bEns(instanceIndex)%scaleFactor(levIndex) > 0.0d0) then
616 bEns(instanceIndex)%scaleFactor(levIndex) = sqrt(bEns(instanceIndex)%scaleFactor(levIndex))
617 else
618 bEns(instanceIndex)%scaleFactor(levIndex) = 0.0d0
619 end if
620 end do
621 bEns(instanceIndex)%scaleFactor_DP(1:bEns(instanceIndex)%nLevEns_DP) = bEns(instanceIndex)%scaleFactor(1:bEns(instanceIndex)%nLevEns_DP)
622 else
623 ! 2D surface variables
624 bEns(instanceIndex)%varNameALFA(:) = varNameALFAsfc(:)
625 if (bEns(instanceIndex)%scaleFactor(1) > 0.0d0) then
626 bEns(instanceIndex)%scaleFactor_SF = sqrt(bEns(instanceIndex)%scaleFactor(1))
627 else
628 call utl_abort('ben_setupOneInstance: with vCode == 0, the scale factor should never be equal to 0')
629 end if
630 end if
631
632 !- 1.5 Domain Partionning
633 call mmpi_setup_latbands(bEns(instanceIndex)%nj, latPerPE, latPerPEmax, bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd)
634 call mmpi_setup_lonbands(bEns(instanceIndex)%ni, lonPerPE, lonPerPEmax, bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd)
635
636 !- 1.6 Localization
637 if ( trim(ben_mode) == 'Analysis' ) then
638
639 call mmpi_setup_levels(bEns(instanceIndex)%nEns,myMemberBeg,myMemberEnd,myMemberCount)
640 call rpn_comm_allreduce(myMemberCount, maxMyMemberCount, &
641 1,"MPI_INTEGER","MPI_MAX","GRID",ierr)
642 bEns(instanceIndex)%nEnsOverDimension = mmpi_npex * maxMyMemberCount
643
644 select case(trim(bEns(instanceIndex)%localizationType))
645 case('LevelDependent')
646 if (mmpi_myid == 0) write(*,*)
647 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: Level-Dependent (Standard) localization will be used'
648 bEns(instanceIndex)%nWaveBand = 1
649
650 case('ScaleDependent','ScaleDependentWithSpectralLoc')
651 if (mmpi_myid == 0) write(*,*)
652 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent') then
653 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: Scale-Dependent localization (SDL) will be used'
654 bEns(instanceIndex)%nWaveBand = count(bEns(instanceIndex)%waveBandPeaks >= 0)
655 bEns(instanceIndex)%nWaveBandForFiltering = bEns(instanceIndex)%nWaveBand
656 else
657 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: Scale-Dependent localization with Spectral localization (SDLwSL) will be used'
658 bEns(instanceIndex)%nWaveBand = 1
659 bEns(instanceIndex)%nWaveBandForFiltering = count(bEns(instanceIndex)%waveBandPeaks >= 0)
660 end if
661
662 if ( bEns(instanceIndex)%nWaveBandForFiltering <= 1 ) then
663 call utl_abort('ben_setupOneInstance: nWaveBandForFiltering <= 1')
664 end if
665 ! You must provide nWaveBand wavenumbers in decreasing order
666 ! e.g. For a 3 wave bands decomposition...
667 ! wavenumber #1 = where the response function for wave band 1 (hgh res) reaches 1
668 ! and stays at 1 for higher wavenumbers
669 ! wavenumber #2 = where the response function for wave band 2 reaches 1
670 ! wavenumber #3 = where the response function for wave band 3 (low res) reaches 1
671 ! and stays at 1 for lower wavenumbers
672 ! See FilterResponseFunction for further info...
673
674 ! Make sure that the wavenumbers are in the correct (decreasing) order
675 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBandForFiltering-1
676 if ( bEns(instanceIndex)%waveBandPeaks(waveBandIndex)-bEns(instanceIndex)%waveBandPeaks(waveBandIndex+1) <= 0 ) then
677 call utl_abort('ben_setupOneInstance: waveBandPeaks are not in decreasing wavenumber order')
678 end if
679 end do
680
681 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent') then
682 ! Make sure that we have valid localization length scales for each wave bands
683 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBand
684 if ( bEns(instanceIndex)%hLocalize(waveBandIndex) <= 0.0d0 ) then
685 call utl_abort('ben_setupOneInstance: Invalid HORIZONTAL localization length scale')
686 end if
687 if ( bEns(instanceIndex)%vLocalize(waveBandIndex) <= 0.0d0 .and. (bEns(instanceIndex)%nLevInc_M > 1 .or. bEns(instanceIndex)%nLevInc_T > 1) ) then
688 call utl_abort('ben_setupOneInstance: Invalid VERTICAL localization length scale')
689 end if
690 end do
691
692 ! Make sure the truncation is compatible with the waveBandPeaks
693 if ( bEns(instanceIndex)%nTrunc < bEns(instanceIndex)%waveBandPeaks(1) ) then
694 call utl_abort('ben_setupOneInstance: The truncation is not compatible with the your scale-dependent localization')
695 end if
696
697 else ! ScaleDependentWithSpectralLoc
698 ! Do we have a valid selected waveBand index?
699 if (bEns(instanceIndex)%waveBandIndexSelected < 1 .or. &
700 bEns(instanceIndex)%waveBandIndexSelected > count(bEns(instanceIndex)%waveBandPeaks >= 0)) then
701 write(*,*) 'ben_setupOneInstance: waveBandIndexSelected = ', bEns(instanceIndex)%waveBandIndexSelected
702 write(*,*) 'ben_setupOneInstance: number of waveBands = ', count(bEns(instanceIndex)%waveBandPeaks >= 0)
703 call utl_abort('ben_setupOneInstance: The selected waveBand index is not valid')
704 end if
705
706 ! Make sure we have only ONE localization length scales for the selected waveBand index
707 if ( bEns(instanceIndex)%hLocalize(2) > 0.0d0 .or. bEns(instanceIndex)%vLocalize(2) > 0.0d0 ) then
708 call utl_abort('ben_setupOneInstance: only a single localization length scale must be provided with SDLwSL')
709 end if
710 end if
711
712 case default
713 call utl_abort('ben_setupOneInstance: Invalid mode for LocalizationType')
714 end select
715
716 ! Setup the localization
717 if ( bEns(instanceIndex)%vco_anl%Vcode == 5002 .or. bEns(instanceIndex)%vco_anl%Vcode == 5005 ) then
718 pSurfRef = 101000.D0
719 call czp_fetch1DLevels(bEns(instanceIndex)%vco_anl, pSurfRef, &
720 profM_opt=vertLocationInc)
721 allocate(vertLocationEns(bEns(instanceIndex)%nLevEns_M))
722 do levIndex = 1, bEns(instanceIndex)%nLevEns_M
723 vertLocationEns(levIndex) = log(vertLocationInc(levIndex+bEns(instanceIndex)%topLevIndex_M-1))
724 end do
725 deallocate(vertLocationInc)
726 else if ( bEns(instanceIndex)%vco_anl%nLev_depth > 0 ) then
727 allocate(vertLocationEns(bEns(instanceIndex)%vco_anl%nLev_depth))
728 vertLocationEns(:) = bEns(instanceIndex)%vco_anl%depths(:)
729 else
730 pSurfRef = 101000.D0
731 allocate(vertLocationEns(1))
732 vertLocationEns(:) = pSurfRef
733 end if
734
735 allocate(bEns(instanceIndex)%locStorage(bEns(instanceIndex)%nWaveBand))
736 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBand
737 call loc_setup(bEns(instanceIndex)%locStorage(waveBandIndex), bEns(instanceIndex)%cvDim_mpilocal, & ! OUT
738 bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_ens, bEns(instanceIndex)%nEns, & ! IN
739 vertLocationEns, bEns(instanceIndex)%nTrunc, 'spectral', & ! IN
740 bEns(instanceIndex)%localizationType, bEns(instanceIndex)%hLocalize(waveBandIndex), & ! IN
741 bEns(instanceIndex)%hLocalize(waveBandIndex+1), bEns(instanceIndex)%vLocalize(waveBandIndex)) ! IN
742 end do
743
744 cvDim = bEns(instanceIndex)%cvDim_mpilocal
745 deallocate(vertLocationEns)
746
747 end if
748
749 !- 1.7 Control variables
750 if ( bEns(instanceIndex)%ctrlVarHumidity == 'LQ' ) then
751 write(*,*)
752 write(*,*) 'ben_setupOneInstance: Humidity control variable = ', bEns(instanceIndex)%ctrlVarHumidity
753 bEns(instanceIndex)%gsvHUcontainsLQ = .true.
754 else if ( bEns(instanceIndex)%ctrlVarHumidity == 'HU' ) then
755 write(*,*)
756 write(*,*) 'ben_setupOneInstance: Humidity control variable = ', bEns(instanceIndex)%ctrlVarHumidity
757 bEns(instanceIndex)%gsvHUcontainsLQ = .false.
758 else
759 write(*,*)
760 write(*,*) 'Unknown humidity control variable'
761 write(*,*) 'Should be LQ or LU, found = ', bEns(instanceIndex)%ctrlVarHumidity
762 call utl_abort('ben_setupOneInstance')
763 end if
764
765 !
766 !- 2. Read/Process the Ensemble
767 !
768
769 !- 2.1 Identify set of variables for which ensembles are required
770 do jvar = 1, vnl_numvarmax
771 if (trim(bEns(instanceIndex)%includeAnlVar(jvar)) == '') exit
772 if (.not.gsv_varExist(varName = trim(bEns(instanceIndex)%includeAnlVar(jvar)))) then
773 write(*,*) 'ben_setupOneInstance: This variable is not a member of ANLVAR: ', trim(bEns(instanceIndex)%includeAnlVar(jvar))
774 call utl_abort('ben_setupOneInstance: Invalid variable in includeAnlVar')
775 else
776 bEns(instanceIndex)%numIncludeAnlVar = bEns(instanceIndex)%numIncludeAnlVar+1
777 end if
778 end do
779 if (bEns(instanceIndex)%numIncludeAnlVar == 0) then
780 do jvar = 1, vnl_numvarmax
781 if (.not.gsv_varExist(varName = trim(vnl_varNamelist(jvar)))) cycle
782 bEns(instanceIndex)%numIncludeAnlVar = bEns(instanceIndex)%numIncludeAnlVar+1
783 bEns(instanceIndex)%includeAnlVar(bEns(instanceIndex)%numIncludeAnlVar) = vnl_varNamelist(jvar)
784 end do
785 end if
786
787 if (bEns(instanceIndex)%numIncludeAnlVar == 0) call utl_abort('ben_setupOneInstance: Ensembles not being requested for any variable')
788
789 bEns(instanceIndex)%ensShouldNotContainLQvarName=.false.
790 if (bEns(instanceIndex)%ctrlVarHumidity == 'LQ' .and. .not. bEns(instanceIndex)%ensContainsFullField) then
791 ! In this particular case, we must force readEnsemble to contains the LQ varName
792 ! to be able to read LQ perturbations
793 do jvar = 1, bEns(instanceIndex)%numIncludeAnlVar
794 if (bEns(instanceIndex)%includeAnlVar(jvar) == 'LQ') then
795 call utl_abort('ben_setup: LQ must not be present in ANLVAR in this case')
796 end if
797 if (bEns(instanceIndex)%includeAnlVar(jvar) == 'HU') then
798 bEns(instanceIndex)%includeAnlVar(jvar) = 'LQ'
799 bEns(instanceIndex)%ensShouldNotContainLQvarName=.true.
800 end if
801 end do
802 end if
803
804 !- 2.2 Read the ensemble data
805 call setupEnsemble(instanceIndex)
806
807 if ( trim(ben_mode) /= 'Analysis' ) then
808 cvDim = 9999 ! Dummy value > 0 to indicate to the background check (s/r ose_compute_HBHT_ensemble)
809 return
810 end if
811
812 !- 2.3 Convert into a C matrix
813 if (bEns(instanceIndex)%useCmatrixOnly .or. trim(bEns(instanceIndex)%varianceSmoothing) /= 'none') then
814 call ens_computeStdDev(bEns(instanceIndex)%ensPerts(1), containsScaledPerts_opt=.true.)
815 call ens_normalize(bEns(instanceIndex)%ensPerts(1))
816 bEns(instanceIndex)%ensPertsNormalized = .true.
817 else
818 bEns(instanceIndex)%ensPertsNormalized = .false.
819 end if
820
821 !- 2.4 Variance smoothing
822 if (trim(bEns(instanceIndex)%varianceSmoothing) /= 'none' .and. .not. bEns(instanceIndex)%useCmatrixOnly) then
823 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: variance smoothing will be performed'
824
825 call ens_copyEnsStdDev(bEns(instanceIndex)%ensPerts(1), bEns(instanceIndex)%statevector_ensStdDev)
826 if (bEns(instanceIndex)%ensDiagnostic) then
827 call gio_writeToFile(bEns(instanceIndex)%statevector_ensStdDev,'./ens_stddev.fst', & ! IN
828 'STDDEV_RAW', HUcontainsLQ_opt=bEns(instanceIndex)%gsvHUcontainsLQ) ! IN
829 end if
830
831 call gsv_power(bEns(instanceIndex)%statevector_ensStdDev, 2.d0 ) ! StdDev -> Variance
832
833 if (trim(bEns(instanceIndex)%varianceSmoothing) == 'horizMean') then
834 if (mmpi_myid == 0) write(*,*) 'ben_setup: variance smoothing type = horizMean'
835 call gbi_setup(gbi_horizontalMean, 'HorizontalMean', &
836 bEns(instanceIndex)%statevector_ensStdDev, bEns(instanceIndex)%hco_core)
837 call gbi_mean(gbi_horizontalMean, bEns(instanceIndex)%statevector_ensStdDev, & ! IN
838 bEns(instanceIndex)%statevector_ensStdDev) ! OUT
839 call gbi_deallocate(gbi_horizontalMean)
840 else if (trim(bEns(instanceIndex)%varianceSmoothing) == 'footprint') then
841 call gsv_smoothHorizontal(bEns(instanceIndex)%statevector_ensStdDev, & ! INOUT
842 bEns(instanceIndex)%footprintRadius) ! IN
843 else if (trim(bEns(instanceIndex)%varianceSmoothing) == 'footprintLandSeaTopo') then
844 call gbi_setup(gbi_landSeaTopo, 'landSeaTopo', bEns(instanceIndex)%statevector_ensStdDev, &
845 bEns(instanceIndex)%hco_core, &
846 mpi_distribution_opt='None', writeBinsToFile_opt=bEns(instanceIndex)%ensDiagnostic)
847 call gsv_getField(gbi_landSeaTopo%statevector_bin2d,bin2d)
848 HeightSfc => gsv_getHeightSfc(gbi_landSeaTopo%statevector_bin2d)
849 call gsv_smoothHorizontal(bEns(instanceIndex)%statevector_ensStdDev, & ! INOUT
850 bEns(instanceIndex)%footprintRadius, binInteger_opt=bin2d, binReal_opt=HeightSfc,& ! IN
851 binRealThreshold_opt=bEns(instanceIndex)%footprintTopoThreshold) ! IN
852 call gbi_deallocate(gbi_landSeaTopo)
853 else
854 call utl_abort('ben_setupOneInstance: Invalid variance smoothing type = '//trim(bEns(instanceIndex)%varianceSmoothing))
855 end if
856
857 call gsv_power(bEns(instanceIndex)%statevector_ensStdDev, 0.5d0) ! Variance -> StdDev
858
859 if (bEns(instanceIndex)%ensDiagnostic) then
860 call gio_writeToFile(bEns(instanceIndex)%statevector_ensStdDev,'./ens_stddev.fst',& ! IN
861 'STDDEV_SMOOT', HUcontainsLQ_opt=bEns(instanceIndex)%gsvHUcontainsLQ) ! IN
862 end if
863 end if
864
865 !- 2.5 Pre-compute everything for advection in FSO mode
866 if (bEns(instanceIndex)%fsoLeadTime > 0.0D0) then
867 bEns(instanceIndex)%amp3dStepIndexFSOFcst = 1
868 if ( sum(bEns(instanceIndex)%advectFactorFSOFcst(:)) == 0.0D0 .or. bEns(instanceIndex)%numStep == 1) then
869 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: advection not activated for FSO'
870 bEns(instanceIndex)%advectAmplitudeFSOFcst = .false.
871 bEns(instanceIndex)%numStepAmplitudeFSOFcst = 1
872 else
873 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: advection activated in FSO mode'
874 bEns(instanceIndex)%advectAmplitudeFSOFcst = .true.
875 bEns(instanceIndex)%numStepAmplitudeFSOFcst = 2
876 bEns(instanceIndex)%numStepAdvectFSOFcst = nint(bEns(instanceIndex)%fsoLeadTime/6.0D0) + 1
877 allocate(advectFactorFSOFcst_M(bEns(instanceIndex)%vco_ens%nLev_M))
878 advectFactorFSOFcst_M(:) = bEns(instanceIndex)%advectFactorFSOFcst(1:bEns(instanceIndex)%vco_ens%nLev_M)
879 allocate(bEns(instanceIndex)%dateStampListAdvectedFields(bEns(instanceIndex)%numStepAmplitudeFSOFcst))
880 bEns(instanceIndex)%dateStampListAdvectedFields(1) = tim_getDatestamp()
881 bEns(instanceIndex)%dateStampListAdvectedFields(2) = bEns(instanceIndex)%dateStampList(bEns(instanceIndex)%numStep)
882 delT_hour = bEns(instanceIndex)%fsoLeadTime/real(bEns(instanceIndex)%numStepAdvectFSOFcst-1,8) ! time between winds
883 call utl_tmg_start(55,'------B_ENS_SetupAdvecFSO')
884 call adv_setup( bEns(instanceIndex)%adv_amplitudeFSOFcst, & ! OUT
885 'fromFirstTimeIndex', bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_ens, & ! IN
886 bEns(instanceIndex)%numStepAmplitudeFSOFcst, bEns(instanceIndex)%dateStampListAdvectedFields, & ! IN
887 bEns(instanceIndex)%numStepAdvectFSOFcst, delT_hour, advectFactorFSOFcst_M, & ! IN
888 'MMLevsOnly', & ! IN
889 steeringFlowFilename_opt=trim(bEns(instanceIndex)%ensPathName)//'/forecast_for_advection' ) ! IN
890 call utl_tmg_stop(55)
891 deallocate(advectFactorFSOFcst_M)
892 end if
893 else
894 bEns(instanceIndex)%numStepAmplitudeFSOFcst = 0
895 end if
896
897 !- 2.6 Pre-compute everything for advection in ANALYSIS mode
898 if ( sum(bEns(instanceIndex)%advectFactorAssimWindow(:)) == 0.0D0 .or. bEns(instanceIndex)%numStep == 1) then
899 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: advection not activated in ANALYSIS mode'
900
901 bEns(instanceIndex)%advectAmplitudeAssimWindow = .false.
902 bEns(instanceIndex)%numStepAmplitudeAssimWindow = 1
903 bEns(instanceIndex)%amp3dStepIndexAssimWindow = 1
904
905 else
906 if (mmpi_myid == 0) write(*,*) 'ben_setupOneInstance: advection activated in ANALYSIS mode'
907
908 delT_hour = tim_dstepobsinc
909 allocate(advectFactorAssimWindow_M(bEns(instanceIndex)%vco_ens%nLev_M))
910 advectFactorAssimWindow_M(:) = bEns(instanceIndex)%advectFactorAssimWindow(1:bEns(instanceIndex)%vco_ens%nLev_M)
911 allocate(bEns(instanceIndex)%dateStampListAdvectedFields(bEns(instanceIndex)%numStep))
912 bEns(instanceIndex)%dateStampListAdvectedFields(:) = bEns(instanceIndex)%dateStampList(:)
913 call gsv_allocate(statevector_ensMean4D, bEns(instanceIndex)%numStep, bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_ens, &
914 datestampList_opt=bEns(instanceIndex)%dateStampListAdvectedFields, &
915 mpi_local_opt=.true., &
916 allocHeight_opt=.false., allocPressure_opt=.false.)
917 call ens_copyEnsMean(bEns(instanceIndex)%ensPerts(1), & ! IN
918 statevector_ensMean4D ) ! OUT
919
920 call utl_tmg_start(56,'------B_ENS_SetupAdvecAnl')
921
922 select case(trim(bEns(instanceIndex)%advectTypeAssimWindow))
923 case ('amplitude')
924 if (mmpi_myid == 0) write(*,*) ' amplitude fields will be advected'
925 bEns(instanceIndex)%advectAmplitudeAssimWindow = .true.
926 bEns(instanceIndex)%numStepAmplitudeAssimWindow = bEns(instanceIndex)%numStep
927 bEns(instanceIndex)%numStepAdvectAssimWindow = bEns(instanceIndex)%numStep
928
929 select case(trim(bEns(instanceIndex)%advectStartTimeIndexAssimWindow))
930 case ('first')
931 direction='fromFirstTimeIndex'
932 bEns(instanceIndex)%amp3dStepIndexAssimWindow = 1
933 case ('middle')
934 direction='fromMiddleTimeIndex'
935 bEns(instanceIndex)%amp3dStepIndexAssimWindow = (bEns(instanceIndex)%numStepAmplitudeAssimWindow+1)/2
936 case default
937 write(*,*)
938 write(*,*) 'Unsupported starting timeIndex : ', trim(bEns(instanceIndex)%advectStartTimeIndexAssimWindow)
939 call utl_abort('ben_setupOneInstance')
940 end select
941 call adv_setup( bEns(instanceIndex)%adv_amplitudeAssimWindow, & ! OUT
942 direction, bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_ens, & ! IN
943 bEns(instanceIndex)%numStepAmplitudeAssimWindow, bEns(instanceIndex)%dateStampListAdvectedFields, & ! IN
944 bEns(instanceIndex)%numStepAdvectAssimWindow, delT_hour, advectFactorAssimWindow_M, & ! IN
945 'MMLevsOnly', statevector_steeringFlow_opt = statevector_ensMean4D ) ! IN
946
947 case('ensPertAnlInc')
948 if (mmpi_myid == 0) write(*,*) ' ensPerts and AnalInc will be advected'
949
950 if (.not. EnsTopMatchesAnlTop) then
951 call utl_abort('ben_setupOneInstance: for advectTypeAssimWindow=ensPertAnlInc, ensTop and anlTop must match!')
952 end if
953
954 bEns(instanceIndex)%advectEnsPertAnlInc = .true.
955 bEns(instanceIndex)%amp3dStepIndexAssimWindow = 1
956 bEns(instanceIndex)%numStepAmplitudeAssimWindow = 1
957 bEns(instanceIndex)%numStepAdvectAssimWindow = bEns(instanceIndex)%numStep
958
959 select case(trim(bEns(instanceIndex)%advectStartTimeIndexAssimWindow))
960 case ('first')
961 directionEnsPerts='towardFirstTimeIndex'
962 directionAnlInc ='towardFirstTimeIndexInverse'
963 case ('middle')
964 directionEnsPerts='towardMiddleTimeIndex'
965 directionAnlInc ='towardMiddleTimeIndexInverse'
966 case default
967 write(*,*)
968 write(*,*) 'Unsupported starting timeIndex : ', trim(bEns(instanceIndex)%advectStartTimeIndexAssimWindow)
969 call utl_abort('ben_setupOneInstance')
970 end select
971
972 call adv_setup( bEns(instanceIndex)%adv_ensPerts, & ! OUT
973 directionEnsPerts, bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_ens, & ! IN
974 bEns(instanceIndex)%numStepAdvectAssimWindow, bEns(instanceIndex)%dateStampListAdvectedFields, & ! IN
975 bEns(instanceIndex)%numStepAdvectAssimWindow, delT_hour, advectFactorAssimWindow_M, & ! IN
976 'allLevs', statevector_steeringFlow_opt=statevector_ensMean4D ) ! IN
977
978 call adv_setup( bEns(instanceIndex)%adv_analInc, & ! OUT
979 directionAnlInc, bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_ens, & ! IN
980 bEns(instanceIndex)%numStepAdvectAssimWindow, bEns(instanceIndex)%dateStampListAdvectedFields, & ! IN
981 bEns(instanceIndex)%numStepAdvectAssimWindow, delT_hour, advectFactorAssimWindow_M, & ! IN
982 'allLevs', statevector_steeringFlow_opt=statevector_ensMean4D ) ! IN
983
984 case default
985 write(*,*)
986 write(*,*) 'Unsupported advectTypeAssimWindow : ', trim(bEns(instanceIndex)%advectTypeAssimWindow)
987 call utl_abort('ben_setupOneInstance')
988 end select
989
990 call utl_tmg_stop(56)
991
992 deallocate(advectFactorAssimWindow_M)
993
994 !- If wanted, write the ensemble mean
995 if (bEns(instanceIndex)%advDiagnostic) then
996 do stepIndex = 1, tim_nstepobsinc
997 call gio_writeToFile(statevector_ensMean4D,'./ens_mean.fst','ENSMEAN4D', & ! IN
998 stepIndex_opt=stepIndex, HUcontainsLQ_opt=bEns(instanceIndex)%gsvHUcontainsLQ ) ! IN
999 end do
1000 end if
1001
1002 call gsv_deallocate(statevector_ensMean4D)
1003
1004 end if
1005
1006 !- 2.7 Ensemble perturbations advection
1007 if ( bEns(instanceIndex)%advectEnsPertAnlInc ) then
1008
1009 !- If wanted, write the original ensemble perturbations for member #1
1010 if (bEns(instanceIndex)%advDiagnostic) then
1011 call ens_copyMember(bEns(instanceIndex)%ensPerts(1), statevector_oneEnsPert4D, 1)
1012 do stepIndex = 1, tim_nstepobsinc
1013 call gio_writeToFile(statevector_oneEnsPert4D,'./ens_pert1.fst','ORIGINAL', & ! IN
1014 stepIndex_opt=stepIndex, HUcontainsLQ_opt=bEns(instanceIndex)%gsvHUcontainsLQ ) ! IN
1015 end do
1016 end if
1017
1018 !- Do the advection of all the members
1019 call adv_ensemble_tl( bEns(instanceIndex)%ensPerts(1), & ! INOUT
1020 bEns(instanceIndex)%adv_ensPerts, bEns(instanceIndex)%nEns ) ! IN
1021
1022 !- If wanted, write the advected ensemble perturbations for member #1
1023 if (bEns(instanceIndex)%advDiagnostic) then
1024 call ens_copyMember(bEns(instanceIndex)%ensPerts(1), statevector_oneEnsPert4D, 1)
1025 do stepIndex = 1, tim_nstepobsinc
1026 call gio_writeToFile(statevector_oneEnsPert4D,'./ens_pert1_advected.fst','ADVECTED', & ! IN
1027 stepIndex_opt=stepIndex,HUcontainsLQ_opt=bEns(instanceIndex)%gsvHUcontainsLQ ) ! IN
1028 end do
1029 end if
1030
1031 end if
1032
1033 !- 2.8 Compute and write Std. Dev.
1034 if (bEns(instanceIndex)%ensDiagnostic) call ensembleDiagnostic(instanceIndex,'FullPerturbations')
1035
1036 !- 2.9 Partitioned the ensemble perturbations into wave bands
1037 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent' .or. &
1038 trim(bEns(instanceIndex)%localizationType) == 'ScaleDependentWithSpectralLoc') then
1039 call ensembleScaleDecomposition(instanceIndex)
1040 if (bEns(instanceIndex)%ensDiagnostic) call ensembleDiagnostic(instanceIndex,'WaveBandPerturbations')
1041 end if
1042
1043 ! Allocate main ensemble Amplitude arrays for improved efficiency - do not know why!!!
1044 if (bEns(instanceIndex)%numStepAmplitudeAssimWindow > 1 .or. bEns(instanceIndex)%numStepAmplitudeFSOFcst > 1 ) then
1045 write(*,*) 'ben_setupOneInstance: WARNING: due to advection being activated, cannot currently used saved ensemble amplitudes!'
1046 bEns(instanceIndex)%useSaveAmp = .false.
1047 else
1048 write(*,*) 'ben_setupOneInstance: using saved memory for ensemble amplitudes for improved efficiency'
1049 bEns(instanceIndex)%useSaveAmp = .true.
1050 call ens_allocate(ensAmplitudeSave(instanceIndex), &
1051 bEns(instanceIndex)%nEnsOverDimension, &
1052 bEns(instanceIndex)%numStepAmplitudeAssimWindow, &
1053 bEns(instanceIndex)%hco_ens, &
1054 bEns(instanceIndex)%vco_ens, &
1055 bEns(instanceIndex)%dateStampList, &
1056 hco_core_opt=bEns(instanceIndex)%hco_core, &
1057 varNames_opt=bEns(instanceIndex)%varNameALFA, dataKind_opt=8)
1058 call ens_zero(ensAmplitudeSave(instanceIndex))
1059 end if
1060
1061 !- 2.10 Setup en ensGridStateVector to store the amplitude fields (for writing)
1062 if (bEns(instanceIndex)%keepAmplitude) then
1063 write(*,*)
1064 write(*,*) 'ben_setupOneInstance: ensAmplitude fields will be store for potential write to file'
1065 call ens_allocate(bEns(instanceIndex)%ensAmplitudeStorage, bEns(instanceIndex)%nEns, &
1066 bEns(instanceIndex)%numStepAmplitudeAssimWindow, &
1067 bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_ens, &
1068 bEns(instanceIndex)%dateStampList, &
1069 hco_core_opt=bEns(instanceIndex)%hco_core, &
1070 varNames_opt=bEns(instanceIndex)%varNameALFA, dataKind_opt=8)
1071 end if
1072
1073 end subroutine ben_setupOneInstance
1074
1075 !--------------------------------------------------------------------------
1076 ! ben_finalize
1077 !--------------------------------------------------------------------------
1078 subroutine ben_finalize()
1079 implicit none
1080
1081 ! Locals:
1082 integer :: waveBandIndex
1083 integer :: instanceIndex
1084
1085 if (verbose) write(*,*) 'Entering ben_Finalize'
1086
1087 do instanceIndex = 1, nInstance
1088 if (bEns(instanceIndex)%initialized) then
1089 write(*,*) 'ben_finalize: deallocating B_ensemble arrays for instance #', instanceIndex
1090 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBand
1091 call ens_deallocate(bEns(instanceIndex)%ensPerts(waveBandIndex))
1092 call loc_finalize(bEns(instanceIndex)%locStorage(waveBandIndex))
1093 end do
1094 deallocate(bEns(instanceIndex)%ensPerts)
1095 if (bEns(instanceIndex)%keepAmplitude) call ens_deallocate(bEns(instanceIndex)%ensAmplitudeStorage)
1096 end if
1097 end do
1098
1099 end subroutine ben_finalize
1100
1101 !--------------------------------------------------------------------------
1102 ! ben_getScaleFactor
1103 !--------------------------------------------------------------------------
1104 subroutine ben_getScaleFactor(scaleFactor_out, instanceIndex_opt)
1105 implicit none
1106
1107 ! Arguments:
1108 real(8), intent(inout) :: scaleFactor_out(:)
1109 integer, optional, intent(in) :: instanceIndex_opt
1110
1111 ! Locals:
1112 integer :: levIndex, instanceIndex
1113
1114 if (verbose) write(*,*) 'Entering ben_getScaleFactor'
1115
1116 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
1117
1118 ! return value of 0 above highest level of ensemble
1119 do levIndex = 1, (bEns(instanceIndex)%topLevIndex_T - 1)
1120 scaleFactor_out(levIndex) = 0.0d0
1121 end do
1122 ! return scale factor for thermo levels
1123 do levIndex = bEns(instanceIndex)%topLevIndex_T, bEns(instanceIndex)%nLevInc_T
1124 scaleFactor_out(levIndex) = bEns(instanceIndex)%scaleFactor_T(levIndex-bEns(instanceIndex)%topLevIndex_T+1)
1125 end do
1126
1127 end subroutine ben_getScaleFactor
1128
1129 !--------------------------------------------------------------------------
1130 ! ben_setInstanceIndex
1131 !--------------------------------------------------------------------------
1132 function ben_setInstanceIndex(instanceIndex_opt) result(instanceIndex)
1133 !
1134 !:Purpose: To return the appropriate instance index
1135 !
1136 implicit none
1137
1138 ! Arguments:
1139 integer, optional, intent(in) :: instanceIndex_opt
1140 ! Result:
1141 integer :: instanceIndex
1142
1143 if (present(instanceIndex_opt)) then
1144 instanceIndex = instanceIndex_opt
1145 else
1146 instanceIndex = 1
1147 end if
1148
1149 end function ben_setInstanceIndex
1150
1151 !--------------------------------------------------------------------------
1152 ! ben_getnEns
1153 !--------------------------------------------------------------------------
1154 integer function ben_getnEns(instanceIndex_opt)
1155 !
1156 !:Purpose: To return the number of ensemble members
1157 !
1158 implicit none
1159
1160 ! Arguments:
1161 integer, optional, intent(in) :: instanceIndex_opt
1162
1163 ! Locals:
1164 integer :: instanceIndex
1165
1166 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
1167
1168 ben_getnEns = bEns(instanceIndex)%nEns
1169
1170 end function ben_getnEns
1171
1172 !--------------------------------------------------------------------------
1173 ! setupEnsemble
1174 !--------------------------------------------------------------------------
1175 subroutine setupEnsemble(instanceIndex)
1176 implicit none
1177
1178 ! Arguments:
1179 integer, intent(in) :: instanceIndex
1180
1181 ! Locals:
1182 real(4), pointer :: ptr4d_r4(:,:,:,:)
1183 real(8) :: multFactor
1184 integer :: stepIndex,levIndex,lev,waveBandIndex,memberIndex,varIndex
1185 logical :: makeBiPeriodic
1186 character(len=4) :: varName
1187 character(len=30) :: transform
1188
1189 write(*,*) 'setupEnsemble: Start'
1190 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1191
1192 !- 1. Memory allocation
1193 allocate(bEns(instanceIndex)%ensPerts(bEns(instanceIndex)%nWaveBand))
1194 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBand
1195 call ens_allocate(bEns(instanceIndex)%ensPerts(waveBandIndex), &
1196 bEns(instanceIndex)%nEns, bEns(instanceIndex)%numStep, &
1197 bEns(instanceIndex)%hco_ens, &
1198 bEns(instanceIndex)%vco_ens, bEns(instanceIndex)%dateStampList, &
1199 hco_core_opt = bEns(instanceIndex)%hco_core, &
1200 varNames_opt = bEns(instanceIndex)%includeAnlVar(1:bEns(instanceIndex)%numIncludeAnlVar), &
1201 hInterpolateDegree_opt = bEns(instanceIndex)%hInterpolationDegree)
1202 end do
1203
1204 !- 2. Read ensemble
1205 makeBiPeriodic = (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent' .or. &
1206 trim(bEns(instanceIndex)%localizationType) == 'ScaleDependentWithSpectralLoc' .or. &
1207 trim(bEns(instanceIndex)%varianceSmoothing) /= 'none')
1208 call ens_readEnsemble(bEns(instanceIndex)%ensPerts(1), bEns(instanceIndex)%ensPathName, makeBiPeriodic, &
1209 vco_file_opt = bEns(instanceIndex)%vco_file, &
1210 varNames_opt = bEns(instanceIndex)%includeAnlVar(1:bEns(instanceIndex)%numIncludeAnlVar), &
1211 containsFullField_opt=bEns(instanceIndex)%ensContainsFullField)
1212
1213 if ( bEns(instanceIndex)%ctrlVarHumidity == 'LQ' .and. ens_varExist(bEns(instanceIndex)%ensPerts(1),'HU') .and. &
1214 bEns(instanceIndex)%ensContainsFullField ) then
1215 call gvt_transform(bEns(instanceIndex)%ensPerts(1),'HUtoLQ',huMinValue_opt=bEns(instanceIndex)%huMinValue)
1216 else if ( bEns(instanceIndex)%ctrlVarHumidity == 'HU' .and. ens_varExist(bEns(instanceIndex)%ensPerts(1),'HU') .and. &
1217 bEns(instanceIndex)%ensContainsFullField .and. bEns(instanceIndex)%huMinValue /= MPC_missingValue_R8 ) then
1218 call qlim_setMin(bEns(instanceIndex)%ensPerts(1), bEns(instanceIndex)%huMinValue)
1219 else if ( trim(bEns(instanceIndex)%transformVarKindCH) /= '' ) then
1220 do varIndex = 1, bEns(instanceIndex)%numIncludeAnlVar
1221 if ( vnl_varKindFromVarname(bEns(instanceIndex)%includeAnlVar(varIndex)) /= 'CH' ) cycle
1222
1223 transform = trim(bens(instanceIndex)%transformVarKindCH)//'CH'
1224 call gvt_transform( bEns(instanceIndex)%ensPerts(1), trim(transform), &
1225 varName_opt=bEns(instanceIndex)%includeAnlVar(varIndex) )
1226 end do
1227 end if
1228
1229 !- 3. From ensemble FORECASTS to ensemble PERTURBATIONS
1230
1231 !- 3.1 remove mean
1232 call ens_computeMean( bEns(instanceIndex)%ensPerts(1), bEns(instanceIndex)%removeSubEnsMeans, numSubEns_opt=bEns(instanceIndex)%numSubEns )
1233 call ens_removeMean( bEns(instanceIndex)%ensPerts(1) )
1234
1235 !- 3.2 normalize and apply scale factors
1236 !$OMP PARALLEL DO PRIVATE (levIndex,varName,lev,ptr4d_r4,stepIndex,memberIndex,multFactor)
1237 do levIndex = 1, ens_getNumK(bEns(instanceIndex)%ensPerts(1))
1238 varName = ens_getVarNameFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
1239 lev = ens_getLevFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
1240
1241 if ( .not. ens_varExist(bEns(instanceIndex)%ensPerts(1), varName) ) cycle
1242
1243 ptr4d_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(1),levIndex)
1244
1245 do stepIndex = 1, bEns(instanceIndex)%numStep
1246 do memberIndex = 1, bEns(instanceIndex)%nEns
1247
1248 if ( vnl_varLevelFromVarname(varName) == 'MM' ) then
1249 multFactor = bEns(instanceIndex)%scaleFactor_M(lev)
1250 else if ( vnl_varLevelFromVarname(varName) == 'TH' ) then
1251 multFactor = bEns(instanceIndex)%scaleFactor_T(lev)
1252 else if ( vnl_varLevelFromVarname(varName) == 'SF' ) then
1253 multFactor = bEns(instanceIndex)%scaleFactor_SF
1254 else if ( vnl_varLevelFromVarname(varName) == 'DP' ) then
1255 multFactor = bEns(instanceIndex)%scaleFactor_DP(lev)
1256 else if ( vnl_varLevelFromVarname(varName) == 'SS' ) then
1257 multFactor = bEns(instanceIndex)%scaleFactor_DP(1)
1258 else
1259 write(*,*) 'varName = ', varName, ', varLevel = ', vnl_varLevelFromVarname(varName)
1260 call utl_abort('setupEnsemble: unknown varLevel')
1261 end if
1262
1263 multFactor = multFactor/sqrt(1.0d0*dble(bEns(instanceIndex)%nEns-bEns(instanceIndex)%numSubEns))
1264
1265 if (trim(varName) == 'HU') then
1266 multFactor = multFactor*bEns(instanceIndex)%scaleFactorHumidity(lev)
1267 end if
1268
1269 ptr4d_r4(memberIndex,stepIndex,:,:) = real( real(ptr4d_r4(memberIndex,stepIndex,:,:),8)*multFactor, 4 )
1270
1271 end do ! memberIndex
1272 end do ! stepIndex
1273
1274 end do ! levIndex
1275 !$OMP END PARALLEL DO
1276
1277 write(*,*) 'ben_setupEnsemble: finished adjusting ensemble members...'
1278
1279 end subroutine setupEnsemble
1280
1281 !--------------------------------------------------------------------------
1282 ! ben_getPerturbation
1283 !--------------------------------------------------------------------------
1284 subroutine ben_getPerturbation(statevector, memberIndexWanted, &
1285 upwardExtrapolationMethod, waveBandIndexWanted_opt, &
1286 undoNormalization_opt, instanceIndex_opt)
1287 implicit none
1288
1289 ! Arguments:
1290 type(struct_gsv), intent(inout) :: statevector
1291 integer, intent(in) :: memberIndexWanted
1292 character(len=*), intent(in) :: upwardExtrapolationMethod
1293 integer, optional, intent(in) :: waveBandIndexWanted_opt
1294 logical, optional, intent(in) :: undoNormalization_opt
1295 integer, optional, intent(in) :: instanceIndex_opt
1296
1297 ! Locals:
1298 integer :: instanceIndex
1299 real(8), pointer :: ptr4d_r8(:,:,:,:)
1300 real(4), pointer :: ensOneLev_r4(:,:,:,:)
1301 real(8) :: dnens2, scaleFactor_MT
1302 logical :: undoNormalization
1303 integer :: waveBandIndex
1304 integer :: lonIndex,latIndex,stepIndex,levIndex,lev,levInc,topLevOffset
1305 character(len=4) :: varName
1306
1307 if (verbose) write(*,*) 'Entering ben_getPerturbation'
1308
1309 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
1310
1311 if ( trim(upwardExtrapolationMethod) /= "ConstantValue" ) then
1312 call utl_abort('ben_getPerturbation : Invalid value for upwardExtrapolationMethod')
1313 end if
1314
1315 if ( present(waveBandIndexWanted_opt) ) then
1316 waveBandIndex = waveBandIndexWanted_opt
1317 else
1318 waveBandIndex = 1
1319 end if
1320
1321 ! set default value for optional argument undoNormalization
1322 if ( present(undoNormalization_opt) ) then
1323 undoNormalization = undoNormalization_opt
1324 else
1325 undoNormalization = .false.
1326 end if
1327
1328 do levIndex = 1, ens_getNumK(bEns(instanceIndex)%ensPerts(1))
1329 varName = ens_getVarNameFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
1330 lev = ens_getLevFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
1331
1332 if (varName == 'LQ' .and. bEns(instanceIndex)%ensShouldNotContainLQvarName) then
1333 call gsv_getField(statevector, ptr4d_r8, 'HU')
1334 else
1335 call gsv_getField(statevector, ptr4d_r8, varName)
1336 end if
1337 ensOneLev_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(waveBandIndex),levIndex)
1338
1339 !$OMP PARALLEL DO PRIVATE(stepIndex,topLevOffset,scaleFactor_MT,levInc,dnens2,latIndex,lonIndex)
1340 do stepIndex = 1, bEns(instanceIndex)%numStep
1341
1342 if ( vnl_varLevelFromVarname(varName) == 'MM' ) then
1343 topLevOffset = bEns(instanceIndex)%topLevIndex_M - 1
1344 scaleFactor_MT = bEns(instanceIndex)%scaleFactor_M(lev)
1345 else if ( vnl_varLevelFromVarname(varName) == 'TH' ) then
1346 topLevOffset = bEns(instanceIndex)%topLevIndex_T - 1
1347 scaleFactor_MT = bEns(instanceIndex)%scaleFactor_T(lev)
1348 else ! SF
1349 topLevOffset = 0
1350 scaleFactor_MT = bEns(instanceIndex)%scaleFactor_SF
1351 end if
1352
1353 levInc = lev + topLevOffset
1354
1355 ! undo the normalization (optional)
1356 if (undoNormalization) then
1357 if (scaleFactor_MT > 0.0d0) then
1358 dnens2 = sqrt(1.0d0*dble(bEns(instanceIndex)%nEns-1))/scaleFactor_MT
1359 else
1360 if (stepIndex == 1) then
1361 write(*,*) 'scalefactor not positive, cannot undo normalization!'
1362 write(*,*) varName,scaleFactor_MT,lev
1363 end if
1364 dnens2 = 0.0d0
1365 end if
1366 if (varName == 'HU ') then
1367 if (bEns(instanceIndex)%scaleFactorHumidity(lev).gt.0.0d0) then
1368 dnens2 = dnens2/bEns(instanceIndex)%scaleFactorHumidity(lev)
1369 else
1370 if (stepIndex == 1) then
1371 write(*,*) 'Humidity scalefactor not positive, cannot undo normalization!'
1372 write(*,*) varName,bEns(instanceIndex)%scaleFactorHumidity(lev),lev
1373 end if
1374 end if
1375 end if
1376 else
1377 dnens2 = 1.0d0
1378 end if
1379
1380 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1381 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1382 ptr4d_r8(lonIndex,latIndex,levInc,stepIndex) = &
1383 dnens2*dble(ensOneLev_r4(memberIndexWanted,stepIndex,lonIndex,latIndex))
1384 end do
1385 end do
1386
1387 if ( topLevOffset > 0 .and. lev == 1) then
1388 ! Fill the gap between the ensemble lid and the analysis lid
1389
1390 ! undo the normalization (optional)
1391 if (undoNormalization) then
1392 if (bEns(instanceIndex)%scaleFactor(1) > 0.0d0) then
1393 dnens2 = sqrt(1.0d0*dble(bEns(instanceIndex)%nEns-1))/bEns(instanceIndex)%scaleFactor(1)
1394 else
1395 if (stepIndex == 1) then
1396 write(*,*) 'scalefactor(top) not positive, cannot undo normalization!'
1397 write(*,*) varName,bEns(instanceIndex)%scaleFactor(1)
1398 end if
1399 dnens2 = 0.0d0
1400 end if
1401 if (varName == 'HU ') then
1402 if (bEns(instanceIndex)%scaleFactorHumidity(1) > 0.0d0) then
1403 dnens2 = dnens2/bEns(instanceIndex)%scaleFactorHumidity(1)
1404 else
1405 if (stepIndex == 1) then
1406 write(*,*) 'Humidity scalefactor(top) not positive, cannot undo normalization!'
1407 write(*,*) varName,bEns(instanceIndex)%scaleFactorHumidity(1)
1408 end if
1409 end if
1410 end if
1411 else
1412 dnens2 = 1.0d0
1413 end if
1414
1415 do levInc = 1, topLevOffset
1416 ! using a constant value
1417 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1418 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1419 ptr4d_r8(lonIndex,latIndex,levInc,stepIndex) = dnens2 * &
1420 dble(ensOneLev_r4(memberIndexWanted,stepIndex,lonIndex,latIndex))
1421 end do
1422 end do
1423 end do
1424
1425 end if ! topLevOffset > 0
1426
1427 end do ! stepIndex
1428 !$OMP END PARALLEL DO
1429
1430 end do ! levIndex
1431
1432 end subroutine ben_getPerturbation
1433
1434 !--------------------------------------------------------------------------
1435 ! ben_getEnsMean
1436 !--------------------------------------------------------------------------
1437 subroutine ben_getEnsMean(statevector, upwardExtrapolationMethod, &
1438 instanceIndex_opt)
1439 implicit none
1440
1441 ! Arguments:
1442 type(struct_gsv), intent(inout) :: statevector
1443 character(len=*), intent(in) :: upwardExtrapolationMethod
1444 integer, optional, intent(in) :: instanceIndex_opt
1445
1446 ! Locals:
1447 real(8), pointer :: ptr4d_out(:,:,:,:)
1448 real(8), pointer :: ensOneLev_mean(:,:,:)
1449 integer :: instanceIndex
1450 integer :: lonIndex,latIndex,stepIndex,levIndex,lev,levInc,topLevOffset
1451 character(len=4) :: varName
1452
1453 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
1454
1455 if (.not. bEns(instanceIndex)%initialized) then
1456 if (mmpi_myid == 0) write(*,*) 'ben_getEnsMean: bMatrixEnsemble not initialized, returning zero vector'
1457 call gsv_zero(statevector)
1458 return
1459 end if
1460
1461 if ( trim(upwardExtrapolationMethod) /= "ConstantValue" ) then
1462 call utl_abort('ben_getEnsMean : Invalid value for upwardExtrapolationMethod')
1463 end if
1464
1465 do levIndex = 1, ens_getNumK(bEns(instanceIndex)%ensPerts(1))
1466 varName = ens_getVarNameFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
1467 lev = ens_getLevFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
1468
1469 if (varName == 'LQ' .and. bEns(instanceIndex)%ensShouldNotContainLQvarName) then
1470 call gsv_getField(statevector, ptr4d_out, 'HU')
1471 else
1472 call gsv_getField(statevector, ptr4d_out, varName)
1473 end if
1474 ensOneLev_mean => ens_getOneLevMean_r8(bEns(instanceIndex)%ensPerts(1), 1, levIndex)
1475
1476 !$OMP PARALLEL DO PRIVATE(stepIndex,topLevOffset,levInc,latIndex,lonIndex)
1477 do stepIndex = 1, bEns(instanceIndex)%numStep
1478
1479 if ( vnl_varLevelFromVarname(varName) == 'MM' ) then
1480 topLevOffset = bEns(instanceIndex)%topLevIndex_M - 1
1481 else if ( vnl_varLevelFromVarname(varName) == 'TH' ) then
1482 topLevOffset = bEns(instanceIndex)%topLevIndex_T - 1
1483 else ! SF
1484 topLevOffset = 0
1485 end if
1486
1487 levInc = lev + topLevOffset
1488
1489 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1490 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1491 ptr4d_out(lonIndex,latIndex,levInc,stepIndex) = ensOneLev_mean(stepIndex,lonIndex,latIndex)
1492 end do
1493 end do
1494
1495 if ( topLevOffset > 0 .and. lev == 1 ) then
1496 ! Fill the gap between the ensemble lid and the analysis lid
1497
1498 do levInc = 1, topLevOffset
1499 ! using a constant value
1500 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1501 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1502 ptr4d_out(lonIndex,latIndex,levInc,stepIndex) = ensOneLev_mean(stepIndex,lonIndex,latIndex)
1503 end do
1504 end do
1505 end do
1506
1507 end if ! topLevOffset > 0
1508
1509 end do ! stepIndex
1510 !$OMP END PARALLEL DO
1511
1512 end do ! levIndex
1513
1514 end subroutine ben_getEnsMean
1515
1516 !--------------------------------------------------------------------------
1517 ! ensembleScaleDecomposition
1518 !--------------------------------------------------------------------------
1519 subroutine ensembleScaleDecomposition(instanceIndex)
1520 implicit none
1521
1522 ! Arguments:
1523 integer, intent(in) :: instanceIndex
1524
1525 ! Locals:
1526 integer :: waveBandIndex, memberindex, stepIndex, levIndex, latIndex, lonIndex
1527 integer :: ila_filter, p, nla_filter, nphase_filter
1528 real(8), allocatable :: ResponseFunction(:,:)
1529 real(8), allocatable :: bandSum(:,:)
1530 real(8) :: totwvnb_r8
1531 real(8), allocatable :: ensPertSP(:,:,:)
1532 real(8), allocatable :: ensPertSPfiltered(:,:,:)
1533 real(8), allocatable :: ensPertGD(:,:,:)
1534 real(4), pointer :: ptr4d_r4(:,:,:,:)
1535 integer, allocatable :: nIndex_vec(:)
1536 integer :: nTrunc
1537 integer :: gstFilterID, mIndex, nIndex, mymBeg, mymEnd, mynBeg, mynEnd, mymSkip, mynSkip
1538 integer :: mymCount, mynCount
1539 integer :: waveBandIndexStart, waveBandIndexEnd
1540 integer :: waveBandIndexLoopStart, waveBandIndexLoopEnd, waveBandIndexLoopDirection
1541 type(struct_lst) :: lst_ben_filter ! Spectral transform Parameters for filtering
1542 character(len=19) :: kind
1543
1544 !
1545 ! --> For SDL
1546 !
1547 ! --- Ensemble Perturbation Data at the Start ---
1548 ! ensPerts(1 ,:) contains the full perturbations
1549 ! ensPerts(2:nWaveBand,:) already allocated but empty
1550 !
1551 ! --- Ensemble Perturbation Data at the End ---
1552 ! ensPerts(nWaveBand,:) contains the largest scales
1553 ! ...
1554 ! ensPerts(1 ,:) contains the smallest scales
1555 !
1556
1557 !
1558 ! --> For SDLwSL
1559 !
1560 ! --- Ensemble Perturbation Data at the Start ---
1561 ! ensPerts(1,:) contains the full perturbations
1562 !
1563 ! --- Ensemble Perturbation Data at the End ---
1564 ! ensPerts(1,:) contains the selected scales
1565 !
1566
1567 if ( mmpi_myid == 0 ) then
1568 write(*,*)
1569 write(*,*) 'Scale decomposition of the ensemble perturbations'
1570 write(*,*) ' number of WaveBands for filtering = ', bEns(instanceIndex)%nWaveBandForFiltering
1571 write(*,*) ' WaveBand Peaks (total wavenumber)...'
1572 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBandForFiltering
1573 write(*,*) waveBandIndex, bEns(instanceIndex)%waveBandPeaks(waveBandIndex)
1574 end do
1575 end if
1576
1577 !
1578 !- Setup a spectral transform for filtering (nk = nEnsOverDimension)
1579 !
1580 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent') then
1581 nTrunc = bEns(instanceIndex)%waveBandPeaks(1)
1582 else ! ScaleDependentWithSpectralLoc
1583 if (bEns(instanceIndex)%waveBandIndexSelected == 1) then
1584 nTrunc = -1 ! no truncation needed to extract the smallest scales
1585 else
1586 nTrunc = bEns(instanceIndex)%waveBandPeaks(bEns(instanceIndex)%waveBandIndexSelected-1)
1587 end if
1588 end if
1589
1590 if (bEns(instanceIndex)%hco_ens%global) then
1591 ! Global mode
1592 gstFilterID = gst_setup(bEns(instanceIndex)%ni,bEns(instanceIndex)%nj,nTrunc,bEns(instanceIndex)%nEnsOverDimension)
1593 if (mmpi_myid == 0) write(*,*) 'ben : returned value of gstFilterID = ',gstFilterID
1594
1595 nla_filter = gst_getNla(gstFilterID)
1596 nphase_filter = 2
1597
1598 allocate(nIndex_vec(nla_filter))
1599 call mmpi_setup_m(nTrunc,mymBeg,mymEnd,mymSkip,mymCount)
1600 call mmpi_setup_n(nTrunc,mynBeg,mynEnd,mynSkip,mynCount)
1601 ila_filter = 0
1602 do mIndex = mymBeg, mymEnd, mymSkip
1603 do nIndex = mynBeg, mynEnd, mynSkip
1604 if (mIndex.le.nIndex) then
1605 ila_filter = ila_filter + 1
1606 nIndex_vec(ila_filter) = nIndex
1607 end if
1608 end do
1609 end do
1610
1611 else
1612 ! LAM mode
1613 call lst_Setup(lst_ben_filter, & ! OUT
1614 bEns(instanceIndex)%ni, bEns(instanceIndex)%nj, bEns(instanceIndex)%hco_ens%dlon, nTrunc, & ! IN
1615 'LatLonMN', maxlevels_opt=bEns(instanceIndex)%nEnsOverDimension, gridDataOrder_opt='kij' ) ! IN
1616
1617 nla_filter = lst_ben_filter%nla
1618 nphase_filter = lst_ben_filter%nphase
1619 end if
1620
1621 !
1622 !- 1. Scale decomposition for every wave band except for wave band #1
1623 !
1624 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent') then
1625 waveBandIndexStart = 2 ! Skip the smallest scales
1626 waveBandIndexEnd = bEns(instanceIndex)%nWaveBand
1627 waveBandIndexLoopStart = waveBandIndexEnd ! Start with the largest scales
1628 waveBandIndexLoopEnd = waveBandIndexStart
1629 waveBandIndexLoopDirection = -1
1630 else ! ScaleDependentWithSpectralLoc
1631 waveBandIndexStart = bEns(instanceIndex)%waveBandIndexSelected
1632 waveBandIndexEnd = bEns(instanceIndex)%waveBandIndexSelected
1633 waveBandIndexLoopStart = waveBandIndexStart
1634 waveBandIndexLoopEnd = waveBandIndexEnd
1635 waveBandIndexLoopDirection = 1
1636 end if
1637
1638 allocate(ResponseFunction(nla_filter,waveBandIndexStart:waveBandIndexEnd))
1639 allocate(ensPertSP(nla_filter,nphase_filter,bEns(instanceIndex)%nEnsOverDimension))
1640 allocate(ensPertSPfiltered(nla_filter,nphase_filter,bEns(instanceIndex)%nEnsOverDimension))
1641 allocate(ensPertGD(bEns(instanceIndex)%nEnsOverDimension,bEns(instanceIndex)%myLonBeg:bEns(instanceIndex)%myLonEnd,bEns(instanceIndex)%myLatBeg:bEns(instanceIndex)%myLatEnd))
1642
1643 ensPertSP (:,:,:) = 0.0d0
1644 ensPertSPfiltered(:,:,:) = 0.0d0
1645
1646 !- 1.1 Pre-compute the response function
1647 do waveBandIndex = waveBandIndexLoopStart, waveBandIndexLoopEnd, waveBandIndexLoopDirection
1648 do ila_filter = 1, nla_filter
1649 if (bEns(instanceIndex)%hco_ens%global) then
1650 totwvnb_r8 = real(nIndex_vec(ila_filter),8)
1651 else
1652 totwvnb_r8 = lst_ben_filter%k_r8(ila_filter)
1653 end if
1654 ResponseFunction(ila_filter,waveBandIndex) = spf_FilterResponseFunction(totwvnb_r8,waveBandIndex, bEns(instanceIndex)%waveBandPeaks, &
1655 bEns(instanceIndex)%nWaveBandForFiltering)
1656 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependentWithSpectralLoc') then
1657 ResponseFunction(ila_filter,waveBandIndex) = sqrt(ResponseFunction(ila_filter,waveBandIndex))
1658 end if
1659 write(*,*) totwvnb_r8, ResponseFunction(ila_filter,waveBandIndex)
1660 end do
1661 end do
1662 if (bEns(instanceIndex)%hco_ens%global) deallocate(nIndex_vec)
1663
1664 do stepIndex = 1, bEns(instanceIndex)%numStep ! Loop on ensemble time bin
1665 do levIndex = 1, ens_getNumK(bEns(instanceIndex)%ensPerts(1)) ! Loop on variables and vertical levels
1666
1667 ptr4d_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(1),levIndex)
1668
1669 !- 1.2 GridPoint space -> Spectral Space
1670 !$OMP PARALLEL DO PRIVATE (latIndex)
1671 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1672 ensPertGD(:,:,latIndex) = 0.0d0
1673 end do
1674 !$OMP END PARALLEL DO
1675 !$OMP PARALLEL DO PRIVATE (memberIndex,latIndex,lonIndex)
1676 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1677 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1678 do memberIndex = 1, bEns(instanceIndex)%nEns
1679 ensPertGD(memberIndex,lonIndex,latIndex) = dble(ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex))
1680 end do
1681 end do
1682 end do
1683 !$OMP END PARALLEL DO
1684 if (bEns(instanceIndex)%hco_ens%global) then
1685 ! Global Mode
1686 call gst_setID(gstFilterID) ! IN
1687 call gst_reespe_kij(ensPertSP, & ! OUT
1688 ensPertGD) ! IN
1689 else
1690 ! LAM mode
1691 kind = 'GridPointToSpectral'
1692 call lst_VarTransform(lst_ben_filter, & ! IN
1693 ensPertSP, & ! OUT
1694 ensPertGD, & ! IN
1695 kind, bEns(instanceIndex)%nEnsOverDimension ) ! IN
1696 end if
1697
1698 !- 1.3 Filtering and transformation back to grid point space
1699 do waveBandIndex = waveBandIndexLoopStart, waveBandIndexLoopEnd, waveBandIndexLoopDirection
1700
1701 ! Filtering
1702 !$OMP PARALLEL DO PRIVATE (memberIndex,p,ila_filter)
1703 do memberIndex = 1, bEns(instanceIndex)%nEns
1704 do p = 1, nphase_filter
1705 do ila_filter = 1, nla_filter
1706 ensPertSPfiltered(ila_filter,p,memberIndex) = &
1707 ensPertSP(ila_filter,p,memberIndex) * ResponseFunction(ila_filter,waveBandIndex)
1708 end do
1709 end do
1710 end do
1711 !$OMP END PARALLEL DO
1712
1713 ! Spectral Space -> GridPoint space
1714 if (bEns(instanceIndex)%hco_ens%global) then
1715 ! Global Mode
1716 call gst_setID(gstFilterID) ! IN
1717 call gst_speree_kij(ensPertSPfiltered, & ! IN
1718 ensPertGD) ! OUT
1719 else
1720 ! LAM mode
1721 kind = 'SpectralToGridPoint'
1722 call lst_VarTransform(lst_ben_filter, & ! IN
1723 ensPertSPfiltered, & ! IN
1724 ensPertGD, & ! OUT
1725 kind, bEns(instanceIndex)%nEnsOverDimension ) ! IN
1726 end if
1727
1728 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent') then
1729 ptr4d_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(waveBandIndex),levIndex)
1730 else
1731 ptr4d_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(1),levIndex)
1732 end if
1733 !$OMP PARALLEL DO PRIVATE (memberIndex,latIndex,lonIndex)
1734 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1735 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1736 do memberIndex = 1, bEns(instanceIndex)%nEns
1737 ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex) = sngl(ensPertGD(memberIndex,lonIndex,latIndex))
1738 end do
1739 end do
1740 end do
1741 !$OMP END PARALLEL DO
1742
1743 end do ! waveBandIndex
1744 end do ! time bins
1745 end do ! variables&levels
1746
1747 deallocate(ensPertGD)
1748 deallocate(ResponseFunction)
1749 deallocate(ensPertSP)
1750 deallocate(ensPertSPfiltered)
1751
1752 !
1753 !- 2. For SDL only, Isolate the smallest scales in waveBandIndex = 1 by difference in grid point space
1754 !
1755 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent') then
1756
1757 allocate(bandSum(bEns(instanceIndex)%myLonBeg:bEns(instanceIndex)%myLonEnd,bEns(instanceIndex)%myLatBeg:bEns(instanceIndex)%myLatEnd))
1758 do stepIndex = 1, bEns(instanceIndex)%numStep
1759 !$OMP PARALLEL DO PRIVATE (memberIndex,levIndex,latIndex,lonIndex,waveBandIndex,bandsum,ptr4d_r4)
1760 do levIndex = 1, ens_getNumK(bEns(instanceIndex)%ensPerts(1))
1761 do memberIndex = 1, bEns(instanceIndex)%nEns
1762 bandSum(:,:) = 0.d0
1763 do waveBandIndex = 2, bEns(instanceIndex)%nWaveBand
1764 ptr4d_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(waveBandIndex),levIndex)
1765 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1766 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1767 bandSum(lonIndex,latIndex) = bandSum(lonIndex,latIndex) + dble(ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex))
1768 end do
1769 end do
1770 end do
1771 ptr4d_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(1),levIndex)
1772 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
1773 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
1774 ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex) = sngl(dble(ptr4d_r4(memberIndex,stepIndex,lonIndex,latIndex)) - bandSum(lonIndex,latIndex))
1775 end do
1776 end do
1777 end do
1778 end do
1779 !$OMP END PARALLEL DO
1780 end do
1781 deallocate(bandSum)
1782
1783 end if
1784
1785 end subroutine ensembleScaleDecomposition
1786
1787 !--------------------------------------------------------------------------
1788 ! ben_reduceToMPILocal
1789 !--------------------------------------------------------------------------
1790 subroutine ben_reduceToMPILocal(cv_mpilocal,cv_mpiglobal,instanceIndex)
1791 implicit none
1792
1793 ! Arguments:
1794 integer, intent(in) :: instanceIndex
1795 real(8), intent(out) :: cv_mpilocal(bEns(instanceIndex)%cvDim_mpilocal)
1796 real(8), intent(in) :: cv_mpiglobal(:)
1797
1798 if (verbose) write(*,*) 'Entering ben_reduceToMPILocal'
1799
1800 call loc_reduceToMPILocal(bEns(instanceIndex)%locStorage(1),cv_mpilocal,cv_mpiglobal)
1801
1802 end subroutine ben_reduceToMPILocal
1803
1804 !--------------------------------------------------------------------------
1805 ! ben_reduceToMPILocal_r4
1806 !--------------------------------------------------------------------------
1807 subroutine ben_reduceToMPILocal_r4(cv_mpilocal,cv_mpiglobal,instanceIndex)
1808 implicit none
1809
1810 ! Arguments:
1811 integer, intent(in) :: instanceIndex
1812 real(4), intent(out) :: cv_mpilocal(bEns(instanceIndex)%cvDim_mpilocal)
1813 real(4), intent(in) :: cv_mpiglobal(:)
1814
1815 if (verbose) write(*,*) 'Entering reduceToMPILocal_r4'
1816
1817 call loc_reduceToMPILocal_r4(bEns(instanceIndex)%locStorage(1),cv_mpilocal,cv_mpiglobal) ! IN
1818
1819 end subroutine ben_reduceToMPILocal_r4
1820
1821 !--------------------------------------------------------------------------
1822 ! ben_expandToMPIGlobal
1823 !--------------------------------------------------------------------------
1824 subroutine ben_expandToMPIGlobal(cv_mpilocal,cv_mpiglobal,instanceIndex)
1825 implicit none
1826
1827 ! Arguments:
1828 integer, intent(in) :: instanceIndex
1829 real(8), intent(in) :: cv_mpilocal(bEns(instanceIndex)%cvDim_mpilocal)
1830 real(8), intent(out) :: cv_mpiglobal(:)
1831
1832 if (verbose) write(*,*) 'Entering ben_expandToMPIGlobal'
1833
1834 call loc_expandToMPIGlobal(bEns(instanceIndex)%locStorage(1), cv_mpilocal, & ! IN
1835 cv_mpiglobal) ! OUT
1836
1837 end subroutine ben_expandToMPIGlobal
1838
1839 !--------------------------------------------------------------------------
1840 ! ben_expandToMPIGlobal_r4
1841 !--------------------------------------------------------------------------
1842 subroutine ben_expandToMPIGlobal_r4(cv_mpilocal,cv_mpiglobal,instanceIndex)
1843 implicit none
1844
1845 ! Arguments:
1846 integer, intent(in) :: instanceIndex
1847 real(4), intent(in) :: cv_mpilocal(bEns(instanceIndex)%cvDim_mpilocal)
1848 real(4), intent(out) :: cv_mpiglobal(:)
1849
1850 if (verbose) write(*,*) 'Entering ben_expandToMPIGlobal_r4'
1851
1852 call loc_expandToMPIGlobal_r4(bEns(instanceIndex)%locStorage(1), cv_mpilocal, & ! IN
1853 cv_mpiglobal) ! OUT
1854
1855 end subroutine ben_expandToMPIGlobal_r4
1856
1857 !--------------------------------------------------------------------------
1858 ! ben_BSqrt
1859 !--------------------------------------------------------------------------
1860 subroutine ben_BSqrt(instanceIndex, controlVector_in, statevector, &
1861 useFSOFcst_opt, stateVectorRef_opt)
1862 implicit none
1863
1864 ! Arguments:
1865 integer, intent(in) :: instanceIndex
1866 real(8), intent(in) :: controlVector_in(bEns(instanceIndex)%cvDim_mpilocal)
1867 type(struct_gsv), intent(inout) :: statevector
1868 logical, optional, intent(in) :: useFSOFcst_opt
1869 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
1870
1871 ! Locals:
1872 type(struct_ens), target :: ensAmplitude
1873 type(struct_ens), pointer :: ensAmplitude_ptr
1874 integer :: ierr, waveBandIndex
1875 integer :: numStepAmplitude, amp3dStepIndex
1876 logical :: immediateReturn
1877 logical :: useFSOFcst
1878
1879 if (mmpi_doBarrier) call rpn_comm_barrier('GRID',ierr)
1880
1881 !
1882 !- 1. Tests
1883 !
1884 if (.not. bEns(instanceIndex)%initialized) then
1885 if (mmpi_myid == 0) write(*,*) 'ben_bsqrt: bMatrixEnsemble not initialized'
1886 return
1887 end if
1888
1889 if (sum(bEns(instanceIndex)%scaleFactor) == 0.0d0) then
1890 if (mmpi_myid == 0) write(*,*) 'ben_bsqrt: scaleFactor=0, skipping bSqrt'
1891 return
1892 end if
1893
1894 ! only check controlVector on proc 0, since may be zero-length on some procs
1895 if (mmpi_myid == 0) then
1896 immediateReturn = .false.
1897 if (maxval(controlVector_in) == 0.0d0 .and. minval(controlVector_in) == 0.0d0) then
1898 write(*,*) 'ben_bsqrt: controlVector=0, skipping bSqrt'
1899 immediateReturn = .true.
1900 end if
1901 end if
1902 call rpn_comm_bcast(immediateReturn, 1, 'MPI_LOGICAL', 0, 'GRID', ierr)
1903 if (immediateReturn) return
1904
1905 if (mmpi_myid == 0) write(*,*) 'ben_bsqrt: starting'
1906 if (mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1907
1908 if (present(useFSOFcst_opt)) then
1909 useFSOFcst = useFSOFcst_opt
1910 else
1911 useFSOFcst = .false.
1912 end if
1913
1914 !
1915 !- 2. Compute the analysis increment from Bens
1916 !
1917 if (verbose) write(*,*) 'ben_bsqrt: allocating ensAmplitude'
1918 if (useFSOFcst) then
1919 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeFSOFcst
1920 amp3dStepIndex = bEns(instanceIndex)%amp3dStepIndexFSOFcst
1921 else
1922 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeAssimWindow
1923 amp3dStepIndex = bEns(instanceIndex)%amp3dStepIndexAssimWindow
1924 end if
1925 if (bEns(instanceIndex)%useSaveAmp) then
1926 ensAmplitude_ptr => ensAmplitudeSave(instanceIndex)
1927 else
1928 call ens_allocate(ensAmplitude, bEns(instanceIndex)%nEnsOverDimension, numStepAmplitude, &
1929 bEns(instanceIndex)%hco_ens, &
1930 bEns(instanceIndex)%vco_ens, bEns(instanceIndex)%dateStampList, &
1931 hco_core_opt=bEns(instanceIndex)%hco_core, &
1932 varNames_opt=bEns(instanceIndex)%varNameALFA, dataKind_opt=8)
1933 ensAmplitude_ptr => ensAmplitude
1934 end if
1935 call gsv_zero(statevector)
1936
1937 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBand ! Loop on WaveBand (for ScaleDependent Localization)
1938
1939 ! 2.1 Compute the ensemble amplitudes
1940 call utl_tmg_start(60,'------LocSpectral_TL')
1941 call loc_Lsqrt(bEns(instanceIndex)%locStorage(waveBandIndex),controlVector_in, & ! IN
1942 ensAmplitude_ptr, & ! OUT
1943 amp3dStepIndex) ! IN
1944 call utl_tmg_stop(60)
1945
1946 ! 2.2 Advect the amplitudes
1947 if (bEns(instanceIndex)%advectAmplitudeFSOFcst .and. useFSOFcst) then
1948 call adv_ensemble_tl( ensAmplitude_ptr, & ! INOUT
1949 bEns(instanceIndex)%adv_amplitudeFSOFcst, bEns(instanceIndex)%nEns ) ! IN
1950 else if (bEns(instanceIndex)%advectAmplitudeAssimWindow .and. .not. useFSOFcst) then
1951 call adv_ensemble_tl( ensAmplitude_ptr, & ! INOUT
1952 bEns(instanceIndex)%adv_amplitudeAssimWindow, bEns(instanceIndex)%nEns ) ! IN
1953 end if
1954
1955 if ( bEns(instanceIndex)%keepAmplitude .and. waveBandIndex == 1 ) then
1956 call ens_copy(ensAmplitude_ptr, bEns(instanceIndex)%ensAmplitudeStorage)
1957 end if
1958
1959 ! 2.3 Compute increment by multiplying amplitudes by member perturbations
1960 call addEnsMember( ensAmplitude_ptr, statevector, & ! INOUT
1961 instanceIndex, waveBandIndex, useFSOFcst ) ! IN
1962
1963 end do ! Loop on WaveBand
1964
1965 if (.not. bEns(instanceIndex)%useSaveAmp) call ens_deallocate(ensAmplitude)
1966
1967 ! 2.4 Apply the Std. Dev (if needed)
1968 if (bEns(instanceIndex)%ensPertsNormalized .and. .not. bEns(instanceIndex)%useCmatrixOnly) then
1969 call gsv_schurProduct(bEns(instanceIndex)%statevector_ensStdDev, & ! IN
1970 statevector) ! INOUT
1971 end if
1972
1973 ! 2.5 Advect Increments
1974 if ( bEns(instanceIndex)%advectEnsPertAnlInc ) then
1975 call adv_statevector_tl( statevector, & ! INOUT
1976 bEns(instanceIndex)%adv_analInc ) ! IN
1977 end if
1978
1979 !
1980 !- 3. Variable transforms
1981 !
1982 if ( bEns(instanceIndex)%ctrlVarHumidity == 'LQ' .and. gsv_varExist(varName='HU') ) then
1983 call gvt_transform( statevector, & ! INOUT
1984 'LQtoHU_tlm', & ! IN
1985 stateVectorRef_opt=stateVectorRef_opt ) ! IN
1986 end if
1987
1988 if (mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
1989 if (mmpi_myid == 0) write(*,*) 'ben_bsqrt: done'
1990
1991 end subroutine ben_BSqrt
1992
1993 !--------------------------------------------------------------------------
1994 ! ben_BSqrtAd
1995 !--------------------------------------------------------------------------
1996 subroutine ben_BSqrtAd(instanceIndex, statevector, controlVector_out, &
1997 useFSOFcst_opt, stateVectorRef_opt)
1998 implicit none
1999
2000 ! Arguments:
2001 integer, intent(in) :: instanceIndex
2002 real(8) , intent(inout) :: controlVector_out(bEns(instanceIndex)%cvDim_mpilocal)
2003 type(struct_gsv), intent(inout) :: statevector
2004 logical, optional, intent(in) :: useFSOFcst_opt
2005 type(struct_gsv), optional, intent(in) :: statevectorRef_opt
2006
2007 ! Locals:
2008 type(struct_ens), target :: ensAmplitude
2009 type(struct_ens), pointer :: ensAmplitude_ptr
2010 integer :: ierr, waveBandIndex
2011 integer :: numStepAmplitude,amp3dStepIndex
2012 logical :: useFSOFcst
2013
2014 !
2015 !- 1. Tests
2016 !
2017 if (mmpi_doBarrier) call rpn_comm_barrier('GRID',ierr)
2018
2019 if (.not. bEns(instanceIndex)%initialized) then
2020 if (mmpi_myid == 0) write(*,*) 'ben_bsqrtad: bMatrixEnsemble not initialized'
2021 return
2022 end if
2023
2024 if (sum(bEns(instanceIndex)%scaleFactor) == 0.0d0) then
2025 if (mmpi_myid == 0) write(*,*) 'ben_bsqrtad: scaleFactor=0, skipping bSqrtAd'
2026 return
2027 end if
2028
2029 if (mmpi_myid == 0) write(*,*) 'ben_bsqrtad: starting'
2030 if (mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2031
2032 if (present(useFSOFcst_opt)) then
2033 useFSOFcst = useFSOFcst_opt
2034 else
2035 useFSOFcst = .false.
2036 end if
2037
2038 !
2039 !- 3. Variable transforms
2040 !
2041 if ( bEns(instanceIndex)%ctrlVarHumidity == 'LQ' .and. gsv_varExist(varName='HU') ) then
2042 call gvt_transform( statevector, & ! INOUT
2043 'LQtoHU_ad', & ! IN
2044 stateVectorRef_opt=stateVectorRef_opt ) ! IN
2045 end if
2046
2047 !
2048 !- 2. Compute the analysis increment from Bens
2049 !
2050
2051 ! 2.5 Advect Increments
2052 if ( bEns(instanceIndex)%advectEnsPertAnlInc ) then
2053 call adv_statevector_ad( statevector, & ! INOUT
2054 bEns(instanceIndex)%adv_analInc ) ! IN
2055 end if
2056
2057 ! 2.4 Apply the Std. Dev (if needed)
2058 if (bEns(instanceIndex)%ensPertsNormalized .and. .not. bEns(instanceIndex)%useCmatrixOnly) then
2059 call gsv_schurProduct(bEns(instanceIndex)%statevector_ensStdDev, & ! IN
2060 statevector) ! INOUT
2061 end if
2062
2063 if (verbose) write(*,*) 'ben_bsqrtAd: allocating ensAmplitude'
2064 if (useFSOFcst) then
2065 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeFSOFcst
2066 amp3dStepIndex = bEns(instanceIndex)%amp3dStepIndexFSOFcst
2067 else
2068 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeAssimWindow
2069 amp3dStepIndex = bEns(instanceIndex)%amp3dStepIndexAssimWindow
2070 end if
2071 if (bEns(instanceIndex)%useSaveAmp) then
2072 ensAmplitude_ptr => ensAmplitudeSave(instanceIndex)
2073 else
2074 call ens_allocate(ensAmplitude, bEns(instanceIndex)%nEnsOverDimension, numStepAmplitude, &
2075 bEns(instanceIndex)%hco_ens, &
2076 bEns(instanceIndex)%vco_ens, bEns(instanceIndex)%dateStampList, &
2077 hco_core_opt=bEns(instanceIndex)%hco_core, &
2078 varNames_opt=bEns(instanceIndex)%varNameALFA, dataKind_opt=8)
2079 ensAmplitude_ptr => ensAmplitude
2080 end if
2081
2082 do waveBandIndex = 1, bEns(instanceIndex)%nWaveBand ! Loop on WaveBand (for ScaleDependent Localization)
2083
2084 ! 2.3 Compute increment by multiplying amplitudes by member perturbations
2085 call addEnsMemberAd( statevector, ensAmplitude_ptr, & ! INOUT
2086 instanceIndex, waveBandIndex, useFSOFcst) ! IN
2087
2088 ! 2.2 Advect the amplitudes
2089 if (bEns(instanceIndex)%advectAmplitudeFSOFcst .and. useFSOFcst) then
2090 call adv_ensemble_ad( ensAmplitude_ptr, & ! INOUT
2091 bEns(instanceIndex)%adv_amplitudeFSOFcst, bEns(instanceIndex)%nEns ) ! IN
2092 else if (bEns(instanceIndex)%advectAmplitudeAssimWindow .and. .not. useFSOFcst) then
2093 call adv_ensemble_ad( ensAmplitude_ptr, & ! INOUT
2094 bEns(instanceIndex)%adv_amplitudeAssimWindow, bEns(instanceIndex)%nEns ) ! IN
2095 end if
2096
2097 ! 2.1 Compute the ensemble amplitudes
2098 call utl_tmg_start(64,'------LocSpectral_AD')
2099 call loc_LsqrtAd(bEns(instanceIndex)%locStorage(waveBandIndex),ensAmplitude_ptr, & ! IN
2100 controlVector_out, & ! OUT
2101 amp3dStepIndex) ! IN
2102 call utl_tmg_stop(64)
2103
2104 end do ! Loop on WaveBand
2105
2106 if (.not. bEns(instanceIndex)%useSaveAmp) call ens_deallocate(ensAmplitude)
2107
2108 if (mmpi_myid == 0) write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2109 if (mmpi_myid == 0) write(*,*) 'ben_bsqrtAd: done'
2110
2111 end subroutine ben_BSqrtAd
2112
2113 !--------------------------------------------------------------------------
2114 ! addEnsMember
2115 !--------------------------------------------------------------------------
2116 subroutine addEnsMember(ensAmplitude, statevector_out, &
2117 instanceIndex, waveBandIndex, useFSOFcst_opt)
2118 implicit none
2119
2120 ! Arguments:
2121 type(struct_ens), intent(in) :: ensAmplitude
2122 type(struct_gsv), intent(inout) :: statevector_out
2123 integer, intent(in) :: instanceIndex
2124 integer, intent(in) :: waveBandIndex
2125 logical, optional, intent(in) :: useFSOFcst_opt
2126
2127 ! Locals:
2128 real(8), pointer :: ensAmplitude_oneLev(:,:,:,:)
2129 real(8), pointer :: ensAmplitude_oneLevM1(:,:,:,:), ensAmplitude_oneLevP1(:,:,:,:)
2130 real(8), allocatable, target :: ensAmplitude_MT(:,:,:,:)
2131 real(8), pointer :: ensAmplitude_MT_ptr(:,:,:,:)
2132 real(4), pointer :: increment_out_r4(:,:,:,:)
2133 real(8), pointer :: increment_out_r8(:,:,:,:)
2134 real(8), allocatable :: increment_out2(:,:,:)
2135 real(4), pointer :: ensMemberAll_r4(:,:,:,:)
2136 integer :: lev, lev2, levIndex, stepIndex, stepIndex_amp, latIndex, lonIndex, topLevOffset, memberIndex
2137 character(len=4) :: varName
2138 logical :: useFSOFcst
2139 integer :: stepIndex2, stepBeg, stepEnd, numStepAmplitude
2140
2141 if (verbose) write(*,*) 'Entering ben_addEnsMember'
2142
2143 call utl_tmg_start(58,'------AddMem_TL')
2144
2145 if (present(useFSOFcst_opt)) then
2146 useFSOFcst = useFSOFcst_opt
2147 else
2148 useFSOFcst = .false.
2149 end if
2150 if (useFSOFcst .and. bEns(instanceIndex)%fsoLeadTime > 0.0d0) then
2151 stepBeg = bEns(instanceIndex)%numStep
2152 stepEnd = stepBeg
2153 if (mmpi_myid == 0) write(*,*) 'ben_bsqrtad: using forecast ensemble stored at timestep ',stepEnd
2154 else
2155 stepBeg = 1
2156 stepEnd = bEns(instanceIndex)%numStepAssimWindow
2157 end if
2158 if (useFSOFcst) then
2159 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeFSOFcst
2160 else
2161 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeAssimWindow
2162 end if
2163
2164 allocate(ensAmplitude_MT(bEns(instanceIndex)%nEns,numStepAmplitude,bEns(instanceIndex)%myLonBeg:bEns(instanceIndex)%myLonEnd,bEns(instanceIndex)%myLatBeg:bEns(instanceIndex)%myLatEnd))
2165 allocate(increment_out2(bEns(instanceIndex)%numStep,bEns(instanceIndex)%myLonBeg:bEns(instanceIndex)%myLonEnd,bEns(instanceIndex)%myLatBeg:bEns(instanceIndex)%myLatEnd))
2166
2167 do levIndex = 1, ens_getNumK(bEns(instanceIndex)%ensPerts(waveBandIndex))
2168
2169 lev = ens_getLevFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
2170 varName = ens_getVarNameFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
2171
2172 !$OMP PARALLEL DO PRIVATE (latIndex)
2173 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
2174 increment_out2(:,:,latIndex) = 0.0d0
2175 end do
2176 !$OMP END PARALLEL DO
2177
2178 if ( vnl_varLevelFromVarname(varName) /= 'SF' .and. &
2179 vnl_varLevelFromVarname(varName) == vnl_varLevelFromVarname(bEns(instanceIndex)%varNameALFA(1)) ) then
2180 ! The non-surface variable varName is on the same levels than the amplitude field
2181 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2182 ensAmplitude_MT_ptr(1:,1:,bEns(instanceIndex)%myLonBeg:,bEns(instanceIndex)%myLatBeg:) => ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,:,:)
2183
2184 else if (vnl_varLevelFromVarname(varName) == 'TH') then
2185 ! The non-surface variable varName is on TH levels whereas the amplitude field is on MM levels
2186
2187 if (bEns(instanceIndex)%vco_anl%Vcode == 5002) then
2188
2189 if (lev == 1) then
2190 ! use top momentum level amplitudes for top thermo level
2191 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2192 ensAmplitude_MT_ptr(1:,1:,bEns(instanceIndex)%myLonBeg:,bEns(instanceIndex)%myLatBeg:) => ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,:,:)
2193 else if (lev == bEns(instanceIndex)%nLevEns_T) then
2194 ! use surface momentum level amplitudes for surface thermo level
2195 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,bEns(instanceIndex)%nLevEns_M)
2196 ensAmplitude_MT_ptr(1:,1:,bEns(instanceIndex)%myLonBeg:,bEns(instanceIndex)%myLatBeg:) => ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,:,:)
2197 else
2198 ! for other levels, interpolate momentum weights to get thermo amplitudes
2199 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2200 ensAmplitude_oneLevM1 => ens_getOneLev_r8(ensAmplitude,lev-1)
2201 !$OMP PARALLEL DO PRIVATE (latIndex)
2202 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
2203 ensAmplitude_MT(:,:,:,latIndex) = 0.5d0*( ensAmplitude_oneLevM1(1:bEns(instanceIndex)%nEns,:,:,latIndex) + &
2204 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,:,latIndex) )
2205 end do
2206 !$OMP END PARALLEL DO
2207 ensAmplitude_MT_ptr(1:,1:,bEns(instanceIndex)%myLonBeg:,bEns(instanceIndex)%myLatBeg:) => ensAmplitude_MT(:,:,:,:)
2208 end if
2209
2210 else if (bEns(instanceIndex)%vco_anl%Vcode == 5005) then
2211
2212 if (lev == bEns(instanceIndex)%nLevEns_T) then
2213 ! use surface momentum level amplitudes for surface thermo level
2214 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,bEns(instanceIndex)%nLevEns_M)
2215 ensAmplitude_MT_ptr(1:,1:,bEns(instanceIndex)%myLonBeg:,bEns(instanceIndex)%myLatBeg:) => ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,:,:)
2216 else
2217 ! for other levels, interpolate momentum weights to get thermo amplitudes
2218 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2219 ensAmplitude_oneLevP1 => ens_getOneLev_r8(ensAmplitude,lev+1)
2220 !$OMP PARALLEL DO PRIVATE (latIndex)
2221 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
2222 ensAmplitude_MT(:,:,:,latIndex) = 0.5d0*( ensAmplitude_oneLevP1(1:bEns(instanceIndex)%nEns,:,:,latIndex) + &
2223 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,:,latIndex) )
2224 end do
2225 !$OMP END PARALLEL DO
2226 ensAmplitude_MT_ptr(1:,1:,bEns(instanceIndex)%myLonBeg:,bEns(instanceIndex)%myLatBeg:) => ensAmplitude_MT(:,:,:,:)
2227 end if
2228
2229 else
2230 call utl_abort('ben_addEnsMember: incompatible vcode')
2231 end if
2232
2233 else if (vnl_varLevelFromVarname(varName) == 'SF') then
2234
2235 ! Surface variable cases (atmosphere or surface only)
2236 if (bEns(instanceIndex)%vco_anl%Vcode == 5002 .or. bEns(instanceIndex)%vco_anl%Vcode == 5005) then
2237 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,bEns(instanceIndex)%nLevEns_M)
2238 else ! vco_anl%Vcode == 0
2239 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,1)
2240 end if
2241 ensAmplitude_MT_ptr(1:,1:,bEns(instanceIndex)%myLonBeg:,bEns(instanceIndex)%myLatBeg:) => ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,:,:)
2242
2243 else if (vnl_varLevelFromVarname(varName) == 'SS') then
2244
2245 ! Surface variable cases (ocean)
2246 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,1)
2247
2248 else
2249
2250 write(*,*) 'variable name = ', varName
2251 write(*,*) 'varLevel = ', vnl_varLevelFromVarname(varName)
2252 call utl_abort('ben_addEnsMember: unknown value of varLevel')
2253
2254 end if
2255
2256 call utl_tmg_start(59,'--------AddMemInner_TL')
2257
2258 ensMemberAll_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(waveBandIndex),levIndex)
2259 !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex,stepIndex,stepIndex2,stepIndex_amp,memberIndex)
2260 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
2261 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
2262 do stepIndex = StepBeg, StepEnd
2263 stepIndex2 = stepIndex - stepBeg + 1
2264 if (bEns(instanceIndex)%advectAmplitudeFSOFcst .and. useFSOFcst) then
2265 stepIndex_amp = 2
2266 else if (bEns(instanceIndex)%advectAmplitudeAssimWindow .and. .not. useFSOFcst) then
2267 stepIndex_amp = stepIndex
2268 else
2269 stepIndex_amp = 1
2270 end if
2271 do memberIndex = 1, bEns(instanceIndex)%nEns
2272 increment_out2(stepIndex2,lonIndex,latIndex) = increment_out2(stepIndex2,lonIndex,latIndex) + &
2273 ensAmplitude_MT_ptr(memberIndex,stepIndex_amp,lonIndex,latIndex) * &
2274 dble(ensMemberAll_r4(memberIndex,stepIndex,lonIndex,latIndex))
2275 end do ! memberIndex
2276 end do ! stepIndex
2277 end do ! lonIndex
2278 end do ! latIndex
2279 !$OMP END PARALLEL DO
2280
2281 call utl_tmg_stop(59)
2282
2283 ! compute increment level from amplitude/member level
2284 if (vnl_varLevelFromVarname(varName) == 'SF') then
2285 topLevOffset = 1
2286 else if (vnl_varLevelFromVarname(varName) == 'MM') then
2287 topLevOffset = bEns(instanceIndex)%topLevIndex_M
2288 else if (vnl_varLevelFromVarname(varName) == 'TH') then
2289 topLevOffset = bEns(instanceIndex)%topLevIndex_T
2290 else if (vnl_varLevelFromVarname(varName) == 'SS') then
2291 topLevOffset = 1
2292 else if (vnl_varLevelFromVarname(varName) == 'DP') then
2293 topLevOffset = 1
2294 else
2295 call utl_abort('ben_addEnsMember: unknown value of varLevel')
2296 end if
2297 lev2 = lev - 1 + topLevOffset
2298
2299 if (varName == 'LQ' .and. bEns(instanceIndex)%ensShouldNotContainLQvarName) then
2300 if (gsv_getDataKind(statevector_out) == 4) then
2301 call gsv_getField(statevector_out, increment_out_r4, 'HU')
2302 else
2303 call gsv_getField(statevector_out, increment_out_r8, 'HU')
2304 end if
2305 else
2306 if (gsv_getDataKind(statevector_out) == 4) then
2307 call gsv_getField(statevector_out, increment_out_r4, varName)
2308 else
2309 call gsv_getField(statevector_out, increment_out_r8, varName)
2310 end if
2311 end if
2312 !$OMP PARALLEL DO PRIVATE (stepIndex, stepIndex2)
2313 do stepIndex = StepBeg, StepEnd
2314 stepIndex2 = stepIndex - StepBeg + 1
2315 if (gsv_getDataKind(statevector_out) == 4) then
2316 increment_out_r4(:,:,lev2,stepIndex2) = increment_out_r4(:,:,lev2,stepIndex2) + increment_out2(stepIndex2,:,:)
2317 else
2318 increment_out_r8(:,:,lev2,stepIndex2) = increment_out_r8(:,:,lev2,stepIndex2) + increment_out2(stepIndex2,:,:)
2319 end if
2320 end do
2321 !$OMP END PARALLEL DO
2322
2323 end do ! levIndex
2324
2325 deallocate(ensAmplitude_MT)
2326 deallocate(increment_out2)
2327
2328 call utl_tmg_stop(58)
2329
2330 end subroutine addEnsMember
2331
2332 !--------------------------------------------------------------------------
2333 ! addEnsMemberAd
2334 !--------------------------------------------------------------------------
2335 subroutine addEnsMemberAd(statevector_in, ensAmplitude, &
2336 instanceIndex, waveBandIndex, useFSOFcst_opt)
2337 implicit none
2338
2339 ! Arguments:
2340 type(struct_ens), intent(inout) :: ensAmplitude
2341 type(struct_gsv), intent(inout) :: statevector_in
2342 integer, intent(in) :: instanceIndex
2343 integer, intent(in) :: waveBandIndex
2344 logical, optional, intent(in) :: useFSOFcst_opt
2345
2346 ! Locals:
2347 real(8), pointer :: ensAmplitude_oneLev(:,:,:,:)
2348 real(8), pointer :: ensAmplitude_oneLevM1(:,:,:,:), ensAmplitude_oneLevP1(:,:,:,:)
2349 real(8), allocatable :: ensAmplitude_MT(:,:)
2350 real(4), pointer :: increment_in_r4(:,:,:,:)
2351 real(8), pointer :: increment_in_r8(:,:,:,:)
2352 real(8), allocatable :: increment_in2(:,:,:)
2353 real(4), pointer :: ensMemberAll_r4(:,:,:,:)
2354 integer :: levIndex, lev, lev2, stepIndex, stepIndex_amp, latIndex, lonIndex, topLevOffset, memberIndex
2355 character(len=4) :: varName
2356 character(len=4) :: varLevel, varLevelAlfa
2357 integer :: stepBeg, stepEnd, stepIndex2, numStepAmplitude
2358 logical :: useFSOFcst
2359
2360 if (verbose) write(*,*) 'Entering ben_addEnsMemberAd'
2361
2362 call utl_tmg_start(62,'------AddMem_AD')
2363
2364 if (present(useFSOFcst_opt)) then
2365 useFSOFcst = useFSOFcst_opt
2366 else
2367 useFSOFcst = .false.
2368 end if
2369
2370 if (useFSOFcst .and. bEns(instanceIndex)%fsoLeadTime > 0.0d0) then
2371 stepBeg = bEns(instanceIndex)%numStep
2372 stepEnd = stepBeg
2373 if (mmpi_myid == 0) write(*,*) 'ben_addEnsMemberAd: using forecast ensemble stored at timestep ',stepEnd
2374 else
2375 stepBeg = 1
2376 stepEnd = bEns(instanceIndex)%numStepAssimWindow
2377 end if
2378 if (useFSOFcst) then
2379 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeFSOFcst
2380 else
2381 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeAssimWindow
2382 end if
2383
2384 allocate(ensAmplitude_MT(bEns(instanceIndex)%nEns,numStepAmplitude))
2385 allocate(increment_in2(bEns(instanceIndex)%numStep,bEns(instanceIndex)%myLonBeg:bEns(instanceIndex)%myLonEnd,bEns(instanceIndex)%myLatBeg:bEns(instanceIndex)%myLatEnd))
2386
2387 ! set output ensemble Amplitude to zero
2388 !$OMP PARALLEL DO PRIVATE (levIndex, ensAmplitude_oneLev)
2389 do levIndex = 1, ens_getNumLev(ensAmplitude,vnl_varLevelFromVarname(bEns(instanceIndex)%varNameALFA(1)))
2390 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,levIndex)
2391 ensAmplitude_oneLev(:,:,:,:) = 0.0d0
2392 end do
2393 !$OMP END PARALLEL DO
2394
2395 do levIndex = 1, ens_getNumK(bEns(instanceIndex)%ensPerts(waveBandIndex))
2396
2397 lev = ens_getLevFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
2398 varName = ens_getVarNameFromK(bEns(instanceIndex)%ensPerts(1),levIndex)
2399 varLevel = vnl_varLevelFromVarname(varName)
2400 varLevelAlfa = vnl_varLevelFromVarname(bEns(instanceIndex)%varNameALFA(1))
2401
2402 ! compute increment level from amplitude/member level
2403 if (varLevel == 'SF') then
2404 topLevOffset = 1
2405 else if (varLevel == 'MM') then
2406 topLevOffset = bEns(instanceIndex)%topLevIndex_M
2407 else if (varLevel == 'TH') then
2408 topLevOffset = bEns(instanceIndex)%topLevIndex_T
2409 else
2410 call utl_abort('ben_addEnsMemberAd: unknown value of varLevel')
2411 end if
2412 lev2 = lev - 1 + topLevOffset
2413
2414 if (varName == 'LQ' .and. bEns(instanceIndex)%ensShouldNotContainLQvarName) then
2415 if (gsv_getDataKind(statevector_in) == 4) then
2416 call gsv_getField(statevector_in, increment_in_r4, 'HU')
2417 else
2418 call gsv_getField(statevector_in, increment_in_r8, 'HU')
2419 end if
2420 else
2421 if (gsv_getDataKind(statevector_in) == 4) then
2422 call gsv_getField(statevector_in, increment_in_r4, varName)
2423 else
2424 call gsv_getField(statevector_in, increment_in_r8, varName)
2425 end if
2426 end if
2427 !$OMP PARALLEL DO PRIVATE (stepIndex, stepIndex2)
2428 do stepIndex = stepBeg, stepEnd
2429 stepIndex2 = stepIndex - stepBeg + 1
2430 if (gsv_getDataKind(statevector_in) == 4) then
2431 increment_in2(stepIndex2,:,:) = increment_in_r4(:,:,lev2,stepIndex2)
2432 else
2433 increment_in2(stepIndex2,:,:) = increment_in_r8(:,:,lev2,stepIndex2)
2434 end if
2435 end do
2436 !$OMP END PARALLEL DO
2437
2438 ensMemberAll_r4 => ens_getOneLev_r4(bEns(instanceIndex)%ensPerts(waveBandIndex),levIndex)
2439 !$OMP PARALLEL DO PRIVATE (latIndex,lonIndex,stepIndex, stepIndex2, stepIndex_amp, &
2440 memberIndex,ensAmplitude_oneLev, ensAmplitude_oneLevM1, &
2441 ensAmplitude_oneLevP1, ensAmplitude_MT)
2442 do latIndex = bEns(instanceIndex)%myLatBeg, bEns(instanceIndex)%myLatEnd
2443 do lonIndex = bEns(instanceIndex)%myLonBeg, bEns(instanceIndex)%myLonEnd
2444
2445 if (omp_get_thread_num() == 0) call utl_tmg_start(63,'--------AddMemInner_AD')
2446 ensAmplitude_MT(:,:) = 0.0d0
2447 do stepIndex = stepBeg, stepEnd
2448 stepIndex2 = stepIndex - stepBeg + 1
2449 if (bEns(instanceIndex)%advectAmplitudeFSOFcst .and. useFSOFcst) then
2450 stepIndex_amp = 2
2451 else if (bEns(instanceIndex)%advectAmplitudeAssimWindow .and. .not. useFSOFcst) then
2452 stepIndex_amp = stepIndex
2453 else
2454 stepIndex_amp = 1
2455 end if
2456 do memberIndex = 1, bEns(instanceIndex)%nEns
2457 ensAmplitude_MT(memberIndex,stepIndex_amp) = ensAmplitude_MT(memberIndex,stepIndex_amp) + &
2458 increment_in2(stepIndex2,lonIndex,latIndex) * dble(ensMemberAll_r4(memberIndex,stepIndex,lonIndex,latIndex))
2459 end do ! memberIndex
2460 end do ! stepIndex
2461 if (omp_get_thread_num() == 0) call utl_tmg_stop(63)
2462
2463 ! transform thermo/momentum level amplitude sensitivites appropriately
2464
2465 if ( varLevel /= 'SF' .and. &
2466 varLevel == varLevelAlfa ) then
2467 ! The non-surface variable varName is on the same levels than the amplitude field
2468
2469 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2470 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2471 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + ensAmplitude_MT(:,:)
2472
2473 else if (varLevel == 'TH') then
2474 ! The non-surface variable varName is on TH levels whereas the amplitude field is on MM levels
2475
2476 if (bEns(instanceIndex)%vco_anl%Vcode == 5002) then
2477
2478 if (lev == 1) then
2479 ! use top momentum level amplitudes for top thermo level
2480 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2481 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2482 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + ensAmplitude_MT(:,:)
2483 else if (lev == bEns(instanceIndex)%nLevEns_T) then
2484 ! use surface momentum level amplitudes for surface thermo level
2485 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,bEns(instanceIndex)%nLevEns_M)
2486 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2487 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + ensAmplitude_MT(:,:)
2488 else
2489 ! for other levels, interpolate momentum weights to get thermo amplitudes
2490 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2491 ensAmplitude_oneLevM1 => ens_getOneLev_r8(ensAmplitude,lev-1)
2492 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2493 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + 0.5d0*ensAmplitude_MT(:,:)
2494 ensAmplitude_oneLevM1(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2495 ensAmplitude_oneLevM1(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + 0.5d0*ensAmplitude_MT(:,:)
2496 end if
2497
2498 else if (bEns(instanceIndex)%vco_anl%Vcode == 5005) then
2499
2500 if (lev == bEns(instanceIndex)%nLevEns_T) then
2501 ! use surface momentum level amplitudes for surface thermo level
2502 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,bEns(instanceIndex)%nLevEns_M)
2503 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2504 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + ensAmplitude_MT(:,:)
2505 else
2506 ! for other levels, interpolate momentum weights to get thermo amplitudes
2507 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,lev)
2508 ensAmplitude_oneLevP1 => ens_getOneLev_r8(ensAmplitude,lev+1)
2509 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2510 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + 0.5d0*ensAmplitude_MT(:,:)
2511 ensAmplitude_oneLevP1(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2512 ensAmplitude_oneLevP1(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + 0.5d0*ensAmplitude_MT(:,:)
2513 end if
2514
2515 else
2516 call utl_abort('ben_addEnsMemberAd: incompatible vcode')
2517 end if
2518
2519 else if (varLevel == 'SF') then
2520 ! Surface variable cases
2521
2522 if (bEns(instanceIndex)%vco_anl%Vcode == 5002 .or. bEns(instanceIndex)%vco_anl%Vcode == 5005) then
2523 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,bEns(instanceIndex)%nLevEns_M)
2524 else ! vco_anl%Vcode == 0
2525 ensAmplitude_oneLev => ens_getOneLev_r8(ensAmplitude,1)
2526 end if
2527 ensAmplitude_oneLev (1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) = &
2528 ensAmplitude_oneLev(1:bEns(instanceIndex)%nEns,:,lonIndex,latIndex) + ensAmplitude_MT(:,:)
2529
2530 else
2531 call utl_abort('ben_addEnsMemberAd: unknown value of varLevel')
2532 end if
2533
2534 end do ! lonIndex
2535 end do ! latIndex
2536 !$OMP END PARALLEL DO
2537
2538 end do ! levIndex
2539
2540 deallocate(ensAmplitude_MT)
2541 deallocate(increment_in2)
2542
2543 call utl_tmg_stop(62)
2544
2545 end subroutine addEnsMemberAd
2546
2547 !--------------------------------------------------------------------------
2548 ! ensembleDiagnostic
2549 !--------------------------------------------------------------------------
2550 subroutine ensembleDiagnostic(instanceIndex, mode)
2551 implicit none
2552
2553 ! Arguments:
2554 integer, intent(in) :: instanceIndex
2555 character(len=*), intent(in) :: mode
2556
2557 ! Locals:
2558 type(struct_gsv) :: statevector, statevector_temp
2559 integer :: nWaveBandToDiagnose, waveBandIndex, memberIndex
2560 real(8) :: dnens2
2561 character(len=48):: fileName
2562 character(len=12):: etiket
2563 character(len=2) :: waveBandNumber
2564 character(len=2) :: instanceNumber
2565
2566 if ( trim(mode) == 'FullPerturbations') then
2567 nWaveBandToDiagnose = 1
2568 else if ( trim(mode) == 'WaveBandPerturbations' ) then
2569 if (trim(bEns(instanceIndex)%localizationType) == 'ScaleDependent') then
2570 nWaveBandToDiagnose = bEns(instanceIndex)%nWaveBand
2571 else
2572 nWaveBandToDiagnose = 1
2573 end if
2574 else
2575 write(*,*)
2576 write(*,*) 'mode = ', trim(mode)
2577 call utl_abort('EnsembleDiagnostic: unknown mode')
2578 end if
2579
2580 if ( mmpi_myid == 0 ) write(*,*)
2581 if ( mmpi_myid == 0 ) write(*,*) 'EnsembleDiagnostic in mode: ', mode
2582
2583 !
2584 !- Write each wave band for a selected member
2585 !
2586 if ( mmpi_myid == 0 ) write(*,*) ' writing perturbations for member 001'
2587 memberIndex = 1
2588 dnens2 = sqrt(1.0d0*dble(bEns(instanceIndex)%nEns-1))
2589 do waveBandIndex = 1, nWaveBandToDiagnose
2590 if ( mmpi_myid == 0 ) write(*,*) ' waveBandIndex = ', waveBandIndex
2591 call gsv_allocate(statevector, tim_nstepobsinc, bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_anl, &
2592 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
2593 allocHeight_opt=.false., allocPressure_opt=.false.)
2594 call ben_getPerturbation( statevector, & ! OUT
2595 memberIndex, & ! IN
2596 'ConstantValue', waveBandIndex, instanceIndex_opt = instanceIndex) ! IN
2597 if ( trim(mode) == 'FullPerturbations') then
2598 etiket = 'PERT001_FULL'
2599 else
2600 write(waveBandNumber,'(I2.2)') waveBandIndex
2601 etiket = 'PERT001_WB' // trim(waveBandNumber)
2602 end if
2603 write(instanceNumber,'(I2.2)') instanceIndex
2604 fileName = './ens_pert001_i' // trim(instanceNumber) // '.fst'
2605
2606 call gio_writeToFile(statevector,fileName,etiket, & ! IN
2607 scaleFactor_opt=dnens2, & ! IN
2608 HUcontainsLQ_opt=bEns(instanceIndex)%gsvHUcontainsLQ ) ! IN
2609 call gsv_deallocate(statevector)
2610 end do
2611
2612 !
2613 !- Compute the standard deviations for each wave band
2614 !
2615 if ( mmpi_myid == 0 ) write(*,*) ' computing Std.Dev.'
2616 call gsv_allocate(statevector_temp, tim_nstepobsinc, bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_anl, &
2617 mpi_local_opt=.true., dataKind_opt=ens_getDataKind(bEns(instanceIndex)%ensPerts(1)), &
2618 allocHeight_opt=.false., allocPressure_opt=.false.)
2619
2620 do waveBandIndex = 1, nWaveBandToDiagnose
2621 if ( mmpi_myid == 0 ) write(*,*) ' waveBandIndex = ', waveBandIndex
2622 call gsv_allocate(statevector, tim_nstepobsinc, bEns(instanceIndex)%hco_ens, bEns(instanceIndex)%vco_anl, &
2623 datestamp_opt=tim_getDatestamp(), mpi_local_opt=.true., &
2624 dataKind_opt=ens_getDataKind(bEns(instanceIndex)%ensPerts(1)), &
2625 allocHeight_opt=.false., allocPressure_opt=.false.)
2626 call gsv_zero(statevector)
2627 do memberIndex = 1, bEns(instanceIndex)%nEns
2628 !- Get normalized perturbations
2629 call ens_copyMember(bEns(instanceIndex)%ensPerts(waveBandIndex), & ! IN
2630 statevector_temp, & ! OUT
2631 memberIndex) ! IN
2632
2633 !- Square
2634 call gsv_power(statevector_temp, & ! INOUT
2635 2.d0) ! IN
2636 !- Sum square values, result in statevector
2637 call gsv_add(statevector_temp, & ! IN
2638 statevector) ! INOUT
2639 end do
2640
2641 !- Convert to StdDev
2642 call gsv_power(statevector, & ! INOUT
2643 0.5d0) ! IN
2644
2645 !- Write to file
2646 if ( trim(mode) == 'FullPerturbations') then
2647 etiket = 'STDDEV_FULL'
2648 else
2649 write(waveBandNumber,'(I2.2)') waveBandIndex
2650 etiket = 'STDDEV_WB' // trim(waveBandNumber)
2651 end if
2652 write(instanceNumber,'(I2.2)') instanceIndex
2653 fileName = './ens_stddev_i' // trim(instanceNumber) // '.fst'
2654
2655 call gio_writeToFile(statevector,fileName,etiket, & ! IN
2656 HUcontainsLQ_opt=bEns(instanceIndex)%gsvHUcontainsLQ) ! IN
2657 call gsv_deallocate(statevector)
2658 end do
2659
2660 call gsv_deallocate(statevector_temp)
2661
2662 end subroutine ensembleDiagnostic
2663
2664 !--------------------------------------------------------------------------
2665 ! ben_writeAmplitude
2666 !--------------------------------------------------------------------------
2667 subroutine ben_writeAmplitude(ensPathName, ensFileNamePrefix, ip3, instanceIndex_opt)
2668 implicit none
2669
2670 ! Arguments:
2671 character(len=*), intent(in) :: ensPathName
2672 character(len=*), intent(in) :: ensFileNamePrefix
2673 integer, intent(in) :: ip3
2674 integer, optional, intent(in) :: instanceIndex_opt
2675
2676 ! Locals:
2677 integer :: instanceIndex
2678
2679 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
2680
2681 if (bEns(instanceIndex)%initialized .and. bEns(instanceIndex)%keepAmplitude) then
2682 if ( mmpi_myid == 0 ) write(*,*)
2683 if ( mmpi_myid == 0 ) write(*,*) 'bmatrixEnsemble_mod: Writing the amplitude field'
2684 call ens_writeEnsemble(bEns(instanceIndex)%ensAmplitudeStorage, ensPathName, ensFileNamePrefix, &
2685 'FROM_BENS', 'R',varNames_opt=bEns(instanceIndex)%varNameALFA, ip3_opt=ip3)
2686 end if
2687
2688 end subroutine ben_writeAmplitude
2689
2690 !--------------------------------------------------------------------------
2691 ! ben_setFsoLeadTime
2692 !--------------------------------------------------------------------------
2693 subroutine ben_setFsoLeadTime(fsoLeadTime_in, instanceIndex_opt)
2694 implicit none
2695
2696 ! Arguments:
2697 real(8), intent(in) :: fsoLeadTime_in
2698 integer, optional, intent(in) :: instanceIndex_opt
2699
2700 ! Locals:
2701 integer :: instanceIndex
2702
2703 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
2704
2705 bEns(instanceIndex)%fsoLeadTime = fsoLeadTime_in
2706
2707 end subroutine ben_setFsoLeadTime
2708
2709 !--------------------------------------------------------------------------
2710 ! ben_getNumStepAmplitudeAssimWindow
2711 !--------------------------------------------------------------------------
2712 function ben_getNumStepAmplitudeAssimWindow(instanceIndex_opt) result(numStepAmplitude)
2713 implicit none
2714
2715 ! Arguments:
2716 integer, optional, intent(in) :: instanceIndex_opt
2717 ! Result:
2718 integer numStepAmplitude
2719
2720 ! Locals:
2721 integer :: instanceIndex
2722
2723 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
2724
2725 numStepAmplitude = bEns(instanceIndex)%numStepAmplitudeAssimWindow
2726
2727 end function ben_getNumStepAmplitudeAssimWindow
2728
2729 !--------------------------------------------------------------------------
2730 ! ben_getAmplitudeAssimWindow
2731 !--------------------------------------------------------------------------
2732 function ben_getAmplitudeAssimWindow(instanceIndex_opt) result(adv_amplitude)
2733 implicit none
2734
2735 ! Arguments:
2736 integer, optional, intent(in) :: instanceIndex_opt
2737 ! Result:
2738 type(struct_adv), pointer :: adv_amplitude
2739
2740 ! Locals:
2741 integer :: instanceIndex
2742
2743 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
2744
2745 adv_amplitude => bEns(instanceIndex)%adv_amplitudeAssimWindow
2746
2747 end function ben_getAmplitudeAssimWindow
2748
2749 !--------------------------------------------------------------------------
2750 ! ben_getAmp3dStepIndexAssimWindow
2751 !--------------------------------------------------------------------------
2752 function ben_getAmp3dStepIndexAssimWindow(instanceIndex_opt) result(stepIndex)
2753 implicit none
2754
2755 ! Arguments:
2756 integer, optional, intent(in) :: instanceIndex_opt
2757 ! Result:
2758 integer :: stepIndex
2759
2760 ! Locals:
2761 integer :: instanceIndex
2762
2763 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
2764
2765 stepIndex = bEns(instanceIndex)%amp3dStepIndexAssimWindow
2766
2767 end function ben_getAmp3dStepIndexAssimWindow
2768
2769 !--------------------------------------------------------------------------
2770 ! ben_getNumInstance
2771 !--------------------------------------------------------------------------
2772 function ben_getNumInstance() result(numInstance)
2773 implicit none
2774
2775 ! Result:
2776 integer :: numInstance
2777
2778 numInstance = nInstance
2779
2780 end function ben_getNumInstance
2781
2782 !--------------------------------------------------------------------------
2783 ! ben_getNumLoc
2784 !--------------------------------------------------------------------------
2785 function ben_getNumLoc(instanceIndex_opt) result(numLoc)
2786 implicit none
2787
2788 ! Arguments:
2789 integer, optional, intent(in) :: instanceIndex_opt
2790 ! Result:
2791 integer :: numLoc
2792
2793 ! Locals:
2794 integer :: instanceIndex
2795
2796 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
2797
2798 numLoc = bEns(instanceIndex)%nWaveBand
2799
2800 end function ben_getNumLoc
2801
2802 !--------------------------------------------------------------------------
2803 ! ben_getLoc
2804 !--------------------------------------------------------------------------
2805 function ben_getLoc(locIndex,instanceIndex_opt) result(loc)
2806 implicit none
2807
2808 ! Arguments:
2809 integer, intent(in) :: locIndex
2810 integer, optional, intent(in) :: instanceIndex_opt
2811 ! Result:
2812 type(struct_loc), pointer :: loc
2813
2814 ! Locals:
2815 integer :: instanceIndex
2816
2817 instanceIndex = ben_setInstanceIndex(instanceIndex_opt)
2818
2819 loc => bEns(instanceIndex)%locStorage(locIndex)
2820
2821 end function ben_getLoc
2822
2823end module bMatrixEnsemble_mod
