1module ensembleStateVector_mod
2 ! MODULE ensembleStateVector_mod (prefix='ens' category='6. High-level data objects')
3 !
4 !:Purpose: Store and manipulate ensemble of state vectors and the ensemble
5 ! mean.
6 !
7 use ramDisk_mod
8 use midasMpi_mod
9 use fileNames_mod
10 use gridStateVector_mod
11 use gridStateVectorFileIO_mod
12 use interpolation_mod
13 use horizontalCoord_mod
14 use verticalCoord_mod
15 use lamAnalysisGridTransforms_mod
16 use oceanMask_mod
17 use timeCoord_mod
18 use utilities_mod
19 use varNameList_mod
20 use codePrecision_mod
21 implicit none
22 save
23 private
24
25 ! public procedures
26 public :: struct_ens, ens_isAllocated, ens_allocate, ens_deallocate, ens_zero
27 public :: ens_readEnsemble, ens_writeEnsemble, ens_copy, ens_copy4Dto3D, ens_add
28 public :: ens_getOneLevMean_r8, ens_modifyVarName
29 public :: ens_varExist, ens_getNumLev, ens_getNumMembers, ens_getNumSubEns
30 public :: ens_computeMean, ens_removeMean, ens_removeGlobalMean, ens_recenter
31 public :: ens_copyEnsMean, ens_copyToEnsMean, ens_copyMember, ens_insertMember
32 public :: ens_computeStdDev, ens_copyEnsStdDev, ens_normalize
33 public :: ens_getMask, ens_copyMaskToGsv
34 public :: ens_getOneLev_r4, ens_getOneLev_r8
35 public :: ens_getOffsetFromVarName, ens_getLevFromK, ens_getVarNameFromK
36 public :: ens_getNumK, ens_getKFromLevVarName, ens_getDataKind, ens_getPathName
37 public :: ens_getVco, ens_getHco, ens_getLatLonBounds, ens_getNumStep
38 public :: ens_varNamesList, ens_applyMaskLAM
39
40 integer,external :: get_max_rss
41
42 type :: struct_oneLev_r4
43 real(4), pointer :: onelevel(:,:,:,:) => null()
44 end type struct_oneLev_r4
45
46 type :: struct_oneLev_r8
47 real(8), pointer :: onelevel(:,:,:,:) => null()
48 end type struct_oneLev_r8
49
50 type :: struct_ens
51 private
52 logical :: allocated = .false.
53 integer :: numMembers
54 integer :: dataKind = 4 ! default value
55 integer :: fileMemberIndex1 = 1 ! first member number in ensemble set
56 type(struct_gsv) :: statevector_work
57 type(struct_hco), pointer :: hco_core
58 type(struct_oneLev_r8), allocatable :: allLev_ensMean_r8(:), allLev_ensStdDev_r8(:)
59 type(struct_oneLev_r4), allocatable :: allLev_r4(:)
60 type(struct_oneLev_r8), allocatable :: allLev_r8(:)
61 logical :: meanIsComputed = .false.
62 logical :: stdDevIsComputed = .false.
63 logical :: meanIsRemoved = .false.
64 integer, allocatable :: subEnsIndexList(:), nEnsSubEns(:)
65 integer :: numSubEns
66 character(len=256) :: enspathname
67 character(len=12) :: hInterpolateDegree='UNSPECIFIED' ! 'LINEAR' or 'CUBIC' or 'NEAREST'
68 end type struct_ens
69
70CONTAINS
71
72 !--------------------------------------------------------------------------
73 ! ens_isAllocated
74 !--------------------------------------------------------------------------
75 function ens_isAllocated(ens) result(isAllocated)
76 !
77 !:Purpose: Return true if object has been allocated
78 !
79 implicit none
80
81 ! Arguments:
82 type(struct_ens), intent(in) :: ens
83 ! Result:
84 logical :: isAllocated
85
86 isAllocated = ens%allocated
87 end function ens_isAllocated
88
89 !--------------------------------------------------------------------------
90 ! ens_allocate
91 !--------------------------------------------------------------------------
92 subroutine ens_allocate(ens, numMembers, numStep, hco_comp, vco_ens, &
93 dateStampList, hco_core_opt, varNames_opt, dataKind_opt, &
94 hInterpolateDegree_opt, fileMemberIndex1_opt)
95 !
96 !:Purpose: Allocate an ensembleStateVector object
97 !
98 implicit none
99
100 ! Arguments:
101 type(struct_ens), intent(inout) :: ens
102 integer, intent(in) :: numMembers
103 integer, intent(in) :: numStep
104 type(struct_hco), pointer, intent(in) :: hco_comp
105 type(struct_hco), pointer, optional, intent(in) :: hco_core_opt
106 type(struct_vco), pointer, intent(in) :: vco_ens
107 integer, intent(in) :: dateStampList(:)
108 character(len=*), optional, intent(in) :: varNames_opt(:)
109 integer, optional, intent(in) :: dataKind_opt
110 integer, optional, intent(in) :: fileMemberIndex1_opt
111 character(len=*), optional, intent(in) :: hInterpolateDegree_opt
112
113 ! Locals:
114 integer :: varLevIndex, lon1, lon2, lat1, lat2, k1, k2
115 character(len=4), pointer :: varNames(:)
116
117 if ( ens%allocated ) then
118 write(*,*) 'ens_allocate: this object is already allocated, deallocating first.'
119 call ens_deallocate( ens )
120 end if
121
122 if ( present(dataKind_opt) ) ens%dataKind = dataKind_opt
123
124 if ( present(fileMemberIndex1_opt) ) ens%fileMemberIndex1 = fileMemberIndex1_opt
125
126 if ( present(hInterpolateDegree_opt) ) then
127 ! set the horizontal interpolation degree
128 ens%hInterpolateDegree = trim(hInterpolateDegree_opt)
129 else
130 ens%hInterpolateDegree = 'LINEAR' ! default, for legacy purposes
131 end if
132
133 nullify(varNames)
134 if (present(varNames_opt)) then
135 allocate(varNames(size(varNames_opt)))
136 varNames(:) = varNames_opt(:)
137 else
138 call gsv_varNamesList(varNames)
139 end if
140
141 call gsv_allocate( ens%statevector_work, &
142 numStep, hco_comp, vco_ens, &
143 datestamplist_opt=dateStampList, mpi_local_opt=.true., &
144 varNames_opt=varNames, dataKind_opt=ens%dataKind, &
145 hInterpolateDegree_opt = hInterpolateDegree_opt)
146
147 lon1 = ens%statevector_work%myLonBeg
148 lon2 = ens%statevector_work%myLonEnd
149 lat1 = ens%statevector_work%myLatBeg
150 lat2 = ens%statevector_work%myLatEnd
151 k1 = ens%statevector_work%mykBeg
152 k2 = ens%statevector_work%mykEnd
153
154 if (ens%dataKind == 8) then
155 allocate( ens%allLev_r8(k1:k2) )
156 do varLevIndex = k1, k2
157 allocate( ens%allLev_r8(varLevIndex)%onelevel(numMembers,numStep,lon1:lon2,lat1:lat2) )
158 end do
159 else if (ens%dataKind == 4) then
160 allocate( ens%allLev_r4(k1:k2) )
161 do varLevIndex = k1, k2
162 allocate( ens%allLev_r4(varLevIndex)%onelevel(numMembers,numStep,lon1:lon2,lat1:lat2) )
163 end do
164 else
165 call utl_abort('ens_allocate: unknown value of datakind')
166 end if
167
168 ens%allocated = .true.
169 ens%numMembers = numMembers
170 if (present(hco_core_opt)) then
171 ens%hco_core => hco_core_opt
172 else
173 ens%hco_core => hco_comp
174 end if
175
176 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
177
178 end subroutine ens_allocate
179
180 !--------------------------------------------------------------------------
181 ! ens_allocateMean
182 !--------------------------------------------------------------------------
183 subroutine ens_allocateMean(ens)
184 !
185 !:Purpose: Allocate the ensemble mean arrays within an ensembleStateVector object
186 !
187 implicit none
188
189 ! Arguments:
190 type(struct_ens), intent(inout) :: ens
191
192 ! Locals:
193 integer :: lon1, lon2, lat1, lat2, k1, k2, varLevIndex, numStep
194
195 lon1 = ens%statevector_work%myLonBeg
196 lon2 = ens%statevector_work%myLonEnd
197 lat1 = ens%statevector_work%myLatBeg
198 lat2 = ens%statevector_work%myLatEnd
199 k1 = ens%statevector_work%mykBeg
200 k2 = ens%statevector_work%mykEnd
201 numStep = ens%statevector_work%numStep
202
203 allocate( ens%allLev_ensMean_r8(k1:k2) )
204 do varLevIndex = k1, k2
205 allocate( ens%allLev_ensMean_r8(varLevIndex)%onelevel(ens%numSubEns,numStep,lon1:lon2,lat1:lat2) )
206 ens%allLev_ensMean_r8(varLevIndex)%onelevel(:,:,:,:) = 0.0d0
207 end do
208
209 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
210
211 end subroutine ens_allocateMean
212
213 !--------------------------------------------------------------------------
214 ! ens_allocateStdDev
215 !--------------------------------------------------------------------------
216 subroutine ens_allocateStdDev(ens)
217 !
218 !:Purpose: Allocate the ensemble stddev arrays within an ensembleStateVector object
219 !
220 implicit none
221
222 ! Arguments:
223 type(struct_ens), intent(inout) :: ens
224
225 ! Locals:
226 integer :: lon1, lon2, lat1, lat2, k1, k2, varLevIndex, numStep
227
228 lon1 = ens%statevector_work%myLonBeg
229 lon2 = ens%statevector_work%myLonEnd
230 lat1 = ens%statevector_work%myLatBeg
231 lat2 = ens%statevector_work%myLatEnd
232 k1 = ens%statevector_work%mykBeg
233 k2 = ens%statevector_work%mykEnd
234 numStep = ens%statevector_work%numStep
235
236 allocate( ens%allLev_ensStdDev_r8(k1:k2) )
237 do varLevIndex = k1, k2
238 allocate( ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,numStep,lon1:lon2,lat1:lat2) )
239 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(:,:,:,:) = 0.0d0
240 end do
241
242 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
243
244 end subroutine ens_allocateStdDev
245
246 !--------------------------------------------------------------------------
247 ! ens_deallocate
248 !--------------------------------------------------------------------------
249 subroutine ens_deallocate( ens )
250 !
251 !:Purpose: Deallocate an ensembleStateVector object
252 !
253 implicit none
254
255 ! Arguments:
256 type(struct_ens), intent(inout) :: ens
257
258 ! Locals:
259 integer :: k1, k2, varLevIndex
260
261 if ( .not. ens%allocated ) return
262
263 k1 = ens%statevector_work%mykBeg
264 k2 = ens%statevector_work%mykEnd
265
266 if (ens%dataKind == 8) then
267 do varLevIndex = k1, k2
268 deallocate( ens%allLev_r8(varLevIndex)%onelevel )
269 end do
270 deallocate( ens%allLev_r8 )
271 else if (ens%dataKind == 4) then
272 do varLevIndex = k1, k2
273 deallocate( ens%allLev_r4(varLevIndex)%onelevel )
274 end do
275 deallocate( ens%allLev_r4 )
276 end if
277
278 if (ens%stdDevIsComputed) then
279 do varLevIndex = k1, k2
280 deallocate( ens%allLev_ensStdDev_r8(varLevIndex)%onelevel )
281 end do
282 deallocate( ens%allLev_ensStdDev_r8 )
283 end if
284
285 if (ens%meanIsComputed) then
286 do varLevIndex = k1, k2
287 deallocate( ens%allLev_ensMean_r8(varLevIndex)%onelevel )
288 end do
289 deallocate( ens%allLev_ensMean_r8 )
290 deallocate( ens%subEnsIndexList )
291 deallocate( ens%nEnsSubEns )
292 end if
293
294 ens%allocated = .false.
295
296 end subroutine ens_deallocate
297
298 !--------------------------------------------------------------------------
299 ! ens_modifyVarName
300 !--------------------------------------------------------------------------
301 subroutine ens_modifyVarName(ens, oldVarName, newVarName)
302 !
303 !:Purpose: Change an existing variable name within the ensemble.
304 ! This is only used when the contents of a variable are
305 ! transformed into another variable **in place**.
306 !
307 implicit none
308
309 ! Arguments:
310 type(struct_ens), intent(inout) :: ens
311 character(len=*), intent(in) :: oldVarName
312 character(len=*), intent(in) :: newVarName
313
314 call gsv_modifyVarName(ens%statevector_work,oldVarName, newVarName)
315
316 end subroutine ens_modifyVarName
317
318 !--------------------------------------------------------------------------
319 ! ens_copy
320 !--------------------------------------------------------------------------
321 subroutine ens_copy(ens_in,ens_out)
322 !
323 !:Purpose: Copy the contents of one ensembleStateVector object into another
324 !
325 implicit none
326
327 ! Arguments:
328 type(struct_ens), intent(in) :: ens_in
329 type(struct_ens), intent(inout) :: ens_out
330
331 ! Locals:
332 integer :: lon1, lon2, lat1, lat2, k1, k2
333 integer :: varLevIndex, stepIndex, latIndex, lonIndex, memberIndex
334
335 if (.not.ens_in%allocated) then
336 call utl_abort('ens_copy: ens_in not yet allocated')
337 end if
338 if (.not.ens_out%allocated) then
339 call utl_abort('ens_copy: ens_out not yet allocated')
340 end if
341
342 call gsv_copyMask(ens_in%statevector_work,ens_out%statevector_work)
343
344 lon1 = ens_out%statevector_work%myLonBeg
345 lon2 = ens_out%statevector_work%myLonEnd
346 lat1 = ens_out%statevector_work%myLatBeg
347 lat2 = ens_out%statevector_work%myLatEnd
348 k1 = ens_out%statevector_work%mykBeg
349 k2 = ens_out%statevector_work%mykEnd
350
351 if ( ens_out%dataKind == 8 .and. ens_in%dataKind == 8 ) then
352
353 !$OMP PARALLEL DO PRIVATE (varLevIndex,stepIndex,latIndex,lonIndex,memberIndex)
354 do varLevIndex = k1, k2
355 do latIndex = lat1, lat2
356 do lonIndex = lon1, lon2
357 do stepIndex = 1, ens_out%statevector_work%numStep
358 do memberIndex = 1, ens_out%numMembers
359 ens_out%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
360 ens_in %allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex)
361 end do
362 end do
363 end do
364 end do
365 end do
366 !$OMP END PARALLEL DO
367
368 else if ( ens_out%dataKind == 4 .and. ens_in%dataKind == 4 ) then
369
370 !$OMP PARALLEL DO PRIVATE (varLevIndex,stepIndex,latIndex,lonIndex,memberIndex)
371 do varLevIndex = k1, k2
372 do latIndex = lat1, lat2
373 do lonIndex = lon1, lon2
374 do stepIndex = 1, ens_out%statevector_work%numStep
375 do memberIndex = 1, ens_out%numMembers
376 ens_out%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
377 ens_in %allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex)
378 end do
379 end do
380 end do
381 end do
382 end do
383 !$OMP END PARALLEL DO
384
385 else
386 call utl_abort('ens_copy: Data type must be the same for both ensembleStatevectors')
387 end if
388
389 end subroutine ens_copy
390
391 !--------------------------------------------------------------------------
392 ! ens_copy4Dto3D
393 !--------------------------------------------------------------------------
394 subroutine ens_copy4Dto3D(ens_in,ens_out)
395 !
396 !:Purpose: Copy contents of a 4D ensemble into a 3D ensemble object by
397 ! extracting the middle time step.
398 !
399 implicit none
400
401 ! Arguments:
402 type(struct_ens), intent(in) :: ens_in
403 type(struct_ens), intent(inout) :: ens_out
404
405 ! Locals:
406 integer :: lon1, lon2, lat1, lat2, k1, k2
407 integer :: varLevIndex, latIndex, lonIndex, memberIndex
408 integer :: numStepIn, numStepOut, middleStepIndex
409
410 if (.not.ens_in%allocated) then
411 call utl_abort('ens_copy4Dto3D: ens_in not yet allocated')
412 end if
413 if (.not.ens_out%allocated) then
414 call utl_abort('ens_copy4Dto3D: ens_out not yet allocated')
415 end if
416
417 call gsv_copyMask(ens_in%statevector_work,ens_out%statevector_work)
418
419 lon1 = ens_out%statevector_work%myLonBeg
420 lon2 = ens_out%statevector_work%myLonEnd
421 lat1 = ens_out%statevector_work%myLatBeg
422 lat2 = ens_out%statevector_work%myLatEnd
423 k1 = ens_out%statevector_work%mykBeg
424 k2 = ens_out%statevector_work%mykEnd
425 numStepIn = ens_in%statevector_work%numStep
426 numStepOut = ens_out%statevector_work%numStep
427
428 if (numStepOut /= 1) call utl_abort('ens_copy4Dto3D: output ensemble must have only 1 timestep')
429 if (numStepIn == 1) then
430 write(*,*) 'ens_copy4Dto3D: WARNING: input ensemble only has 1 timestep, will simply copy.'
431 end if
432 middleStepIndex = (numStepIn + 1) / 2
433
434 if ( ens_out%dataKind == 8 .and. ens_in%dataKind == 8 ) then
435
436 !$OMP PARALLEL DO PRIVATE (varLevIndex,latIndex,lonIndex,memberIndex)
437 do varLevIndex = k1, k2
438 do latIndex = lat1, lat2
439 do lonIndex = lon1, lon2
440 do memberIndex = 1, ens_out%numMembers
441 ens_out%allLev_r8(varLevIndex)%onelevel(memberIndex,1,lonIndex,latIndex) = &
442 ens_in %allLev_r8(varLevIndex)%onelevel(memberIndex,middleStepIndex,lonIndex,latIndex)
443 end do
444 end do
445 end do
446 end do
447 !$OMP END PARALLEL DO
448
449 else if ( ens_out%dataKind == 4 .and. ens_in%dataKind == 4 ) then
450
451 !$OMP PARALLEL DO PRIVATE (varLevIndex,latIndex,lonIndex,memberIndex)
452 do varLevIndex = k1, k2
453 do latIndex = lat1, lat2
454 do lonIndex = lon1, lon2
455 do memberIndex = 1, ens_out%numMembers
456 ens_out%allLev_r4(varLevIndex)%onelevel(memberIndex,1,lonIndex,latIndex) = &
457 ens_in %allLev_r4(varLevIndex)%onelevel(memberIndex,middleStepIndex,lonIndex,latIndex)
458 end do
459 end do
460 end do
461 end do
462 !$OMP END PARALLEL DO
463
464 else
465 call utl_abort('ens_copy4Dto3D: Data type must be the same for both ensembleStatevectors')
466 end if
467
468 end subroutine ens_copy4Dto3D
469
470 !--------------------------------------------------------------------------
471 ! ens_add
472 !--------------------------------------------------------------------------
473 subroutine ens_add(ens_in, ens_inOut, scaleFactorIn_opt, scaleFactorInOut_opt)
474 !
475 !:Purpose: Add the contents of the ens_in ensemble to the ens_inOut ensemble.
476 !
477 implicit none
478
479 ! arguments
480 type(struct_ens), intent(in) :: ens_in
481 type(struct_ens), intent(inout) :: ens_inOut
482 real(8), optional, intent(in) :: scaleFactorIn_opt
483 real(8), optional, intent(in) :: scaleFactorInOut_opt
484
485 ! locals
486 integer :: lon1, lon2, lat1, lat2, k1, k2
487 integer :: varLevIndex, stepIndex, latIndex, lonIndex, memberIndex
488 real(4) :: scaleFactorIn_r4, scaleFactorInOut_r4
489 real(8) :: scaleFactorIn, scaleFactorInOut
490
491 if (.not.ens_in%allocated) then
492 call utl_abort('ens_add: ens_in not yet allocated')
493 end if
494 if (.not.ens_inOut%allocated) then
495 call utl_abort('ens_add: ens_inOut not yet allocated')
496 end if
497
498 lon1 = ens_inOut%statevector_work%myLonBeg
499 lon2 = ens_inOut%statevector_work%myLonEnd
500 lat1 = ens_inOut%statevector_work%myLatBeg
501 lat2 = ens_inOut%statevector_work%myLatEnd
502 k1 = ens_inOut%statevector_work%mykBeg
503 k2 = ens_inOut%statevector_work%mykEnd
504
505 if ( ens_inOut%dataKind == 8 .and. ens_in%dataKind == 8 ) then
506
507 if (present(scaleFactorIn_opt)) then
508 scaleFactorIn = scaleFactorIn_opt
509 else
510 scaleFactorIn = 1.0d0
511 end if
512 if (present(scaleFactorInOut_opt)) then
513 scaleFactorInOut = scaleFactorInOut_opt
514 else
515 scaleFactorInOut = 1.0d0
516 end if
517
518 !$OMP PARALLEL DO PRIVATE (varLevIndex,stepIndex,latIndex,lonIndex,memberIndex)
519 do varLevIndex = k1, k2
520 do latIndex = lat1, lat2
521 do lonIndex = lon1, lon2
522 do stepIndex = 1, ens_inOut%statevector_work%numStep
523 do memberIndex = 1, ens_inOut%numMembers
524 ens_inOut%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
525 scaleFactorInOut * &
526 ens_inOut%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) + &
527 scaleFactorIn * &
528 ens_in%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex)
529 end do
530 end do
531 end do
532 end do
533 end do
534 !$OMP END PARALLEL DO
535
536 else if ( ens_inOut%dataKind == 4 .and. ens_in%dataKind == 4 ) then
537
538 if (present(scaleFactorIn_opt)) then
539 scaleFactorIn_r4 = real(scaleFactorIn_opt,4)
540 else
541 scaleFactorIn_r4 = 1.0
542 end if
543 if (present(scaleFactorInOut_opt)) then
544 scaleFactorInOut_r4 = real(scaleFactorInOut_opt,4)
545 else
546 scaleFactorInOut_r4 = 1.0
547 end if
548
549 !$OMP PARALLEL DO PRIVATE (varLevIndex,stepIndex,latIndex,lonIndex,memberIndex)
550 do varLevIndex = k1, k2
551 do latIndex = lat1, lat2
552 do lonIndex = lon1, lon2
553 do stepIndex = 1, ens_inOut%statevector_work%numStep
554 do memberIndex = 1, ens_inOut%numMembers
555 ens_inOut%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
556 scaleFactorInOut_r4 * &
557 ens_inOut%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) + &
558 scaleFactorIn_r4 * &
559 ens_in%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex)
560 end do
561 end do
562 end do
563 end do
564 end do
565 !$OMP END PARALLEL DO
566
567 else
568 call utl_abort('ens_add: Data type must be the same for both ensembleStatevectors')
569 end if
570
571 end subroutine ens_add
572
573 !--------------------------------------------------------------------------
574 ! ens_zero
575 !--------------------------------------------------------------------------
576 subroutine ens_zero(ens)
577 !
578 !:Purpose: Set the main contents of an ensembleStateVector object to zero
579 !
580 implicit none
581
582 ! Arguments:
583 type(struct_ens), intent(inout) :: ens
584
585 ! Locals:
586 integer :: lon1, lon2, lat1, lat2, k1, k2
587 integer :: varLevIndex, stepIndex, latIndex, lonIndex, memberIndex
588
589 if (.not.ens%allocated) then
590 call utl_abort('ens_zero: ens not yet allocated')
591 end if
592
593 lon1 = ens%statevector_work%myLonBeg
594 lon2 = ens%statevector_work%myLonEnd
595 lat1 = ens%statevector_work%myLatBeg
596 lat2 = ens%statevector_work%myLatEnd
597 k1 = ens%statevector_work%mykBeg
598 k2 = ens%statevector_work%mykEnd
599
600 if ( ens%dataKind == 8 ) then
601
602 !$OMP PARALLEL DO PRIVATE (varLevIndex,stepIndex,latIndex,lonIndex,memberIndex)
603 do varLevIndex = k1, k2
604 do latIndex = lat1, lat2
605 do lonIndex = lon1, lon2
606 do stepIndex = 1, ens%statevector_work%numStep
607 do memberIndex = 1, ens%numMembers
608 ens%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = 0.d0
609 end do
610 end do
611 end do
612 end do
613 end do
614 !$OMP END PARALLEL DO
615
616 else if ( ens%dataKind == 4 ) then
617
618 !$OMP PARALLEL DO PRIVATE (varLevIndex,stepIndex,latIndex,lonIndex,memberIndex)
619 do varLevIndex = k1, k2
620 do latIndex = lat1, lat2
621 do lonIndex = lon1, lon2
622 do stepIndex = 1, ens%statevector_work%numStep
623 do memberIndex = 1, ens%numMembers
624 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = 0.0
625 end do
626 end do
627 end do
628 end do
629 end do
630 !$OMP END PARALLEL DO
631
632 end if
633
634 end subroutine ens_zero
635
636 !--------------------------------------------------------------------------
637 ! ens_copyToStateWork (private routine)
638 !--------------------------------------------------------------------------
639 subroutine ens_copyToStateWork(ens, dataType, memberIndex_opt, &
640 subEnsIndex_opt)
641 !
642 !:Purpose: Copy the selected contents of an ensembleStateVector into
643 ! the 'work' stateVector within the object. The possible
644 ! types of content are: 'member', 'mean', or 'stdDev'.
645 !
646 implicit none
647
648 ! Arguments:
649 type(struct_ens), intent(inout) :: ens
650 character(len=*), intent(in) :: dataType
651 integer, optional, intent(in) :: memberIndex_opt
652 integer, optional, intent(in) :: subEnsIndex_opt
653
654 ! Locals:
655 real(4), pointer :: ptr4d_r4(:,:,:,:)
656 real(8), pointer :: ptr4d_r8(:,:,:,:)
657 integer :: k1, k2, varLevIndex, numStep, stepIndex
658 integer :: memberIndex, subEnsIndex
659
660 k1 = ens%statevector_work%mykBeg
661 k2 = ens%statevector_work%mykEnd
662 numStep = ens%statevector_work%numStep
663
664 select case(trim(dataType))
665 case ('member')
666 if (present(memberIndex_opt)) then
667 memberIndex = memberIndex_opt
668 else
669 call utl_abort('ens_copyToStateWork: memberIndex_opt must be provided with dataType=member')
670 end if
671 case ('mean','stdDev')
672 if (present(subEnsIndex_opt)) then
673 subEnsIndex = subEnsIndex_opt
674 else
675 subEnsIndex = 1
676 end if
677 case default
678 write(*,*)
679 write(*,*) 'Unsupported dataType: ', trim(dataType)
680 write(*,*) ' please select either: member, mean or stdDev'
681 call utl_abort('ens_copyToStateWork')
682 end select
683
684 if (ens%dataKind == 8) then
685 call gsv_getField(ens%statevector_work,ptr4d_r8)
686 do stepIndex = 1, numStep
687 do varLevIndex = k1, k2
688 if (dataType == 'member') then
689 ptr4d_r8(:,:,varLevIndex,stepIndex) = &
690 ens%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,:,:)
691 else if (dataType == 'mean') then
692 ptr4d_r8(:,:,varLevIndex,stepIndex) = &
693 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:)
694 else if (dataType == 'stdDev') then
695 ptr4d_r8(:,:,varLevIndex,stepIndex) = &
696 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:)
697 end if
698 end do
699 end do
700 else if (ens%dataKind == 4) then
701 call gsv_getField(ens%statevector_work,ptr4d_r4)
702 do stepIndex = 1, numStep
703 do varLevIndex = k1, k2
704 if (dataType == 'member') then
705 ptr4d_r4(:,:,varLevIndex,stepIndex) = &
706 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,:,:)
707 else if (dataType == 'mean') then
708 ptr4d_r4(:,:,varLevIndex,stepIndex) = &
709 real(ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:),4)
710 else if (dataType == 'stdDev') then
711 ptr4d_r4(:,:,varLevIndex,stepIndex) = &
712 real(ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:),4)
713 end if
714 end do
715 end do
716 end if
717
718 end subroutine ens_copyToStateWork
719
720 !--------------------------------------------------------------------------
721 ! ens_copyFromStateWork
722 !--------------------------------------------------------------------------
723 subroutine ens_copyFromStateWork(ens, dataType, memberIndex_opt, &
724 subEnsIndex_opt)
725 !
726 !:Purpose: This is the inverse operation as ens_copyToStateWork.
727 !
728 implicit none
729
730 ! Arguments:
731 type(struct_ens), intent(inout) :: ens
732 character(len=*), intent(in) :: dataType
733 integer, optional, intent(in) :: memberIndex_opt
734 integer, optional, intent(in) :: subEnsIndex_opt
735
736 ! Locals:
737 real(4), pointer :: ptr4d_r4(:,:,:,:)
738 real(8), pointer :: ptr4d_r8(:,:,:,:)
739 integer :: k1, k2, varLevIndex, numStep, stepIndex
740 integer :: memberIndex, subEnsIndex
741
742 k1 = ens%statevector_work%mykBeg
743 k2 = ens%statevector_work%mykEnd
744 numStep = ens%statevector_work%numStep
745
746 select case(trim(dataType))
747 case ('member')
748 if (present(memberIndex_opt)) then
749 memberIndex = memberIndex_opt
750 else
751 call utl_abort('ens_copyFromStateWork: memberIndex_opt must be provided with dataType=member')
752 end if
753 case ('mean','stdDev')
754 if (present(subEnsIndex_opt)) then
755 subEnsIndex = subEnsIndex_opt
756 else
757 subEnsIndex = 1
758 end if
759 case default
760 write(*,*)
761 write(*,*) 'Unsupported dataType: ', trim(dataType)
762 write(*,*) ' please select either: member, mean or stdDev'
763 call utl_abort('ens_copyToStateWork')
764 end select
765
766 if (ens%dataKind == 8) then
767 call gsv_getField(ens%statevector_work,ptr4d_r8)
768 do stepIndex = 1, numStep
769 do varLevIndex = k1, k2
770 if (dataType == 'member') then
771 ens%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
772 ptr4d_r8(:,:,varLevIndex,stepIndex)
773 else if (dataType == 'mean') then
774 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:) = &
775 ptr4d_r8(:,:,varLevIndex,stepIndex)
776 else if (dataType == 'stdDev') then
777 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:) = &
778 ptr4d_r8(:,:,varLevIndex,stepIndex)
779 end if
780 end do
781 end do
782 else if (ens%dataKind == 4) then
783 call gsv_getField(ens%statevector_work,ptr4d_r4)
784 do stepIndex = 1, numStep
785 do varLevIndex = k1, k2
786 if (dataType == 'member') then
787 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
788 ptr4d_r4(:,:,varLevIndex,stepIndex)
789 else if (dataType == 'mean') then
790 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:) = &
791 real(ptr4d_r4(:,:,varLevIndex,stepIndex),8)
792 else if (dataType == 'stdDev') then
793 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:) = &
794 real(ptr4d_r4(:,:,varLevIndex,stepIndex),8)
795 end if
796 end do
797 end do
798 end if
799
800 end subroutine ens_copyFromStateWork
801
802 !--------------------------------------------------------------------------
803 ! ens_getOneLev_r4
804 !--------------------------------------------------------------------------
805 function ens_getOneLev_r4(ens,kIndex) result(oneLevLevel)
806 !
807 !:Purpose: Return a 4D pointer to a single level of a real4 ensemble. The
808 ! dimensions of the pointer are:
809 ! (memberIndex, stepIndex, lonIndex, latIndex)
810 !
811 implicit none
812
813 ! Arguments:
814 type(struct_ens), intent(inout) :: ens
815 integer, intent(in) :: kIndex
816 ! Result:
817 real(4), pointer :: oneLevLevel(:,:,:,:)
818
819 ! Locals:
820 integer :: lon1, lat1
821
822 lon1 = ens%statevector_work%myLonBeg
823 lat1 = ens%statevector_work%myLatBeg
824
825 oneLevLevel(1:,1:,lon1:,lat1:) => ens%allLev_r4(kIndex)%onelevel(:,:,:,:)
826
827 end function ens_getOneLev_r4
828
829 !--------------------------------------------------------------------------
830 ! ens_getOneLev_r8
831 !--------------------------------------------------------------------------
832 function ens_getOneLev_r8(ens,kIndex) result(oneLevLevel)
833 !
834 !:Purpose: Return a 4D pointer to a single level of a real8 ensemble. The
835 ! dimensions of the pointer are:
836 ! (memberIndex, stepIndex, lonIndex, latIndex)
837 !
838 implicit none
839
840 ! Arguments:
841 type(struct_ens), intent(in) :: ens
842 integer , intent(in) :: kIndex
843 ! Result:
844 real(8), pointer :: oneLevLevel(:,:,:,:)
845
846 ! Locals:
847 integer :: lon1, lat1
848
849 lon1 = ens%statevector_work%myLonBeg
850 lat1 = ens%statevector_work%myLatBeg
851
852 oneLevLevel(1:,1:,lon1:,lat1:) => ens%allLev_r8(kIndex)%onelevel(:,:,:,:)
853
854 end function ens_getOneLev_r8
855
856 !--------------------------------------------------------------------------
857 ! ens_getOneLevMean_r8
858 !--------------------------------------------------------------------------
859 function ens_getOneLevMean_r8(ens,subEnsIndex,kIndex) result(field)
860 !
861 !:Purpose: Return a 3D pointer to a single level the ensemble mean. The
862 ! dimensions of the pointer are:
863 ! (stepIndex, lonIndex, latIndex)
864 !
865 implicit none
866
867 ! Arguments:
868 type(struct_ens), intent(inout) :: ens
869 integer, intent(in) :: subEnsIndex, kIndex
870 ! Result:
871 real(8), pointer :: field(:,:,:)
872
873 ! Locals:
874 integer :: lon1, lat1
875
876 lon1 = ens%statevector_work%myLonBeg
877 lat1 = ens%statevector_work%myLatBeg
878
879 field(1:, lon1:, lat1:) => ens%allLev_ensMean_r8(kIndex)%onelevel(subEnsIndex,:,:,:)
880
881 end function ens_getOneLevMean_r8
882
883 !--------------------------------------------------------------------------
884 ! ens_copyEnsMean
885 !--------------------------------------------------------------------------
886 subroutine ens_copyEnsMean(ens, statevector, subEnsIndex_opt)
887 !
888 !:Purpose: Copy the ensemble mean into the supplied stateVector.
889 !
890 implicit none
891
892 ! Arguments:
893 type(struct_ens), intent(inout) :: ens
894 type(struct_gsv), intent(inout) :: statevector
895 integer, optional, intent(in) :: subEnsIndex_opt
896
897 ! Locals:
898 real(8), pointer :: ptr4d_r8(:,:,:,:)
899 real(4), pointer :: ptr4d_r4(:,:,:,:)
900 integer :: k1, k2, varLevIndex, stepIndex, numStep, subEnsIndex
901 character(len=4), pointer :: varNamesInEns(:)
902
903 if( present(subEnsIndex_opt) ) then
904 subEnsIndex = subEnsIndex_opt
905 else
906 subEnsIndex = 1
907 end if
908
909 k1 = ens%statevector_work%mykBeg
910 k2 = ens%statevector_work%mykEnd
911 numStep = ens%statevector_work%numStep
912
913 if (.not. gsv_isAllocated(statevector)) then
914 nullify(varNamesInEns)
915 call gsv_varNamesList(varNamesInEns,ens%statevector_work)
916 call gsv_allocate(statevector, numStep, &
917 ens%statevector_work%hco, ens%statevector_work%vco, &
918 varNames_opt=varNamesInEns, datestamp_opt=tim_getDatestamp(), &
919 mpi_local_opt=.true., dataKind_opt=8 )
920 deallocate(varNamesInEns)
921 end if
922
923 statevector%onPhysicsGrid(:) = ens%statevector_work%onPhysicsGrid
924 statevector%hco_physics => ens%statevector_work%hco_physics
925
926 if (statevector%dataKind == 8) then
927 call gsv_getField(statevector,ptr4d_r8)
928 do stepIndex = 1, numStep
929 do varLevIndex = k1, k2
930 ptr4d_r8(:,:,varLevIndex,stepIndex) = &
931 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:)
932 end do
933 end do
934 else
935 call gsv_getField(statevector,ptr4d_r4)
936 do stepIndex = 1, numStep
937 do varLevIndex = k1, k2
938 ptr4d_r4(:,:,varLevIndex,stepIndex) = &
939 real(ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:),4)
940 end do
941 end do
942 end if
943
944 end subroutine ens_copyEnsMean
945
946 !--------------------------------------------------------------------------
947 ! ens_copyToEnsMean
948 !--------------------------------------------------------------------------
949 subroutine ens_copyToEnsMean(ens, statevector, subEnsIndex_opt)
950 !
951 !:Purpose: Copy the supplied stateVector into the ensemble mean.
952 !
953 implicit none
954
955 ! Arguments:
956 type(struct_ens), intent(inout) :: ens
957 type(struct_gsv), intent(inout) :: statevector
958 integer, optional, intent(in) :: subEnsIndex_opt
959
960 ! Locals:
961 real(8), pointer :: ptr4d_r8(:,:,:,:)
962 real(4), pointer :: ptr4d_r4(:,:,:,:)
963 integer :: k1, k2, varLevIndex, stepIndex, numStep, subEnsIndex
964
965 if( present(subEnsIndex_opt) ) then
966 subEnsIndex = subEnsIndex_opt
967 else
968 subEnsIndex = 1
969 end if
970
971 k1 = ens%statevector_work%mykBeg
972 k2 = ens%statevector_work%mykEnd
973 numStep = ens%statevector_work%numStep
974
975 if (.not. gsv_isAllocated(statevector)) then
976 call utl_abort('ens_copyToEnsMean: supplied stateVector must be allocated')
977 end if
978
979 if (.not. allocated(ens%allLev_ensMean_r8)) then
980 call ens_allocateMean(ens)
981 else
982 do varLevIndex = k1, k2
983 ens%allLev_ensMean_r8(varLevIndex)%onelevel(:,:,:,:) = 0.0d0
984 end do
985 end if
986 ens%meanIsComputed = .true.
987
988 if (statevector%dataKind == 8) then
989 call gsv_getField(statevector,ptr4d_r8)
990 do stepIndex = 1, numStep
991 do varLevIndex = k1, k2
992 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:) = &
993 ptr4d_r8(:,:,varLevIndex,stepIndex)
994 end do
995 end do
996 else
997 call gsv_getField(statevector,ptr4d_r4)
998 do stepIndex = 1, numStep
999 do varLevIndex = k1, k2
1000 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,:,:) = &
1001 real(ptr4d_r4(:,:,varLevIndex,stepIndex),8)
1002 end do
1003 end do
1004 end if
1005
1006 end subroutine ens_copyToEnsMean
1007
1008 !--------------------------------------------------------------------------
1009 ! ens_copyEnsStdDev
1010 !--------------------------------------------------------------------------
1011 subroutine ens_copyEnsStdDev(ens, statevector)
1012 !
1013 !:Purpose: Copy the ensemble StdDev into the supplied stateVector.
1014 !
1015 implicit none
1016
1017 ! Arguments:
1018 type(struct_ens), intent(inout) :: ens
1019 type(struct_gsv), intent(inout) :: statevector
1020
1021 ! Locals:
1022 real(8), pointer :: ptr4d_r8(:,:,:,:)
1023 real(4), pointer :: ptr4d_r4(:,:,:,:)
1024 integer :: k1, k2, varLevIndex, stepIndex, numStep
1025 character(len=4), pointer :: varNamesInEns(:)
1026
1027 k1 = ens%statevector_work%mykBeg
1028 k2 = ens%statevector_work%mykEnd
1029 numStep = ens%statevector_work%numStep
1030
1031 if (.not. gsv_isAllocated(statevector)) then
1032 nullify(varNamesInEns)
1033 call gsv_varNamesList(varNamesInEns,ens%statevector_work)
1034 call gsv_allocate(statevector, numStep, &
1035 ens%statevector_work%hco, ens%statevector_work%vco, &
1036 varNames_opt=varNamesInEns, datestamp_opt=tim_getDatestamp(), &
1037 mpi_local_opt=.true., dataKind_opt=8 )
1038 deallocate(varNamesInEns)
1039 end if
1040
1041 if (statevector%dataKind == 8) then
1042 call gsv_getField(statevector,ptr4d_r8)
1043 do stepIndex = 1, numStep
1044 do varLevIndex = k1, k2
1045 ptr4d_r8(:,:,varLevIndex,stepIndex) = &
1046 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,:,:)
1047 end do
1048 end do
1049 else
1050 call gsv_getField(statevector,ptr4d_r4)
1051 do stepIndex = 1, numStep
1052 do varLevIndex = k1, k2
1053 ptr4d_r4(:,:,varLevIndex,stepIndex) = &
1054 real(ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,:,:),4)
1055 end do
1056 end do
1057 end if
1058
1059 end subroutine ens_copyEnsStdDev
1060
1061 !--------------------------------------------------------------------------
1062 ! ens_copyMember
1063 !--------------------------------------------------------------------------
1064 subroutine ens_copyMember(ens, statevector, memberIndex)
1065 !
1066 !:Purpose: Copy a selected ensemble member into the supplied stateVector.
1067 !
1068 implicit none
1069
1070 ! Arguments:
1071 type(struct_ens), intent(in) :: ens
1072 type(struct_gsv), intent(inout) :: statevector
1073 integer, intent(in) :: memberIndex
1074
1075 ! Locals:
1076 real(8), pointer :: ptr4d_r8(:,:,:,:)
1077 real(4), pointer :: ptr4d_r4(:,:,:,:)
1078 integer :: k1, k2, varLevIndex, stepIndex, numStep, varIndex
1079 integer :: gsvLevIndex, ensVarLevIndex, nLev
1080 character(len=4), pointer :: varNamesInEns(:)
1081 character(len=4), pointer :: varNamesInGsv(:)
1082 character(len=4) :: varName
1083 logical :: sameVariables
1084
1085 numStep = ens%statevector_work%numStep
1086
1087 nullify(varNamesInEns)
1088 call gsv_varNamesList(varNamesInEns, ens%statevector_work)
1089
1090 if (.not. gsv_isAllocated(statevector)) then
1091 call utl_abort('ens_copyMember: statevector not allocated')
1092 else
1093 nullify(varNamesInGsv)
1094 call gsv_varNamesList(varNamesInGsv, statevector)
1095 end if
1096
1097 sameVariables = .false.
1098 if (size(ens%statevector_work%varExistlist) == size(statevector%varExistlist)) then
1099 if (all(ens%statevector_work%varExistlist == statevector%varExistlist)) then
1100 sameVariables = .true.
1101 end if
1102 end if
1103
1104 if (sameVariables) then
1105
1106 k1 = ens%statevector_work%mykBeg
1107 k2 = ens%statevector_work%mykEnd
1108
1109 if (ens%dataKind == 8) then
1110 call gsv_getField(statevector,ptr4d_r8)
1111 do stepIndex = 1, numStep
1112 do varLevIndex = k1, k2
1113 ptr4d_r8(:,:,varLevIndex,stepIndex) = &
1114 ens%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,:,:)
1115 end do
1116 end do
1117 else if (ens%dataKind == 4) then
1118 if (gsv_getDataKind(statevector) == 8) then
1119 call gsv_getField(statevector,ptr4d_r8)
1120 do stepIndex = 1, numStep
1121 do varLevIndex = k1, k2
1122 ptr4d_r8(:,:,varLevIndex,stepIndex) = &
1123 real(ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,:,:),8)
1124 end do
1125 end do
1126 else
1127 call gsv_getField(statevector,ptr4d_r4)
1128 do stepIndex = 1, numStep
1129 do varLevIndex = k1, k2
1130 ptr4d_r4(:,:,varLevIndex,stepIndex) = &
1131 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,:,:)
1132 end do
1133 end do
1134 end if
1135 end if
1136
1137 else
1138
1139 do varIndex = 1, size(varNamesInGsv)
1140 varName = varNamesInGsv(varIndex)
1141 if (.not. ens_varExist(ens,varName)) cycle
1142 nLev = gsv_getNumLev(statevector,vnl_varLevelFromVarname(varName),varName)
1143 if (ens%dataKind == 8) then
1144 call gsv_getField(statevector,ptr4d_r8,varName_opt=varName)
1145 do stepIndex = 1, numStep
1146 do gsvLevIndex = 1, nLev
1147 ensVarLevIndex = gsvLevIndex + ens_getOffsetFromVarName(ens,varName)
1148 ptr4d_r8(:,:,gsvLevIndex,stepIndex) = &
1149 ens%allLev_r8(ensVarLevIndex)%onelevel(memberIndex,stepIndex,:,:)
1150 end do
1151 end do
1152 else if (ens%dataKind == 4) then
1153 if (gsv_getDataKind(statevector) == 8) then
1154 call gsv_getField(statevector,ptr4d_r8,varName_opt=varName)
1155 do stepIndex = 1, numStep
1156 do gsvLevIndex = 1, nLev
1157 ensVarLevIndex = gsvLevIndex + ens_getOffsetFromVarName(ens,varName)
1158 ptr4d_r8(:,:,gsvLevIndex,stepIndex) = &
1159 real(ens%allLev_r4(ensVarLevIndex)%onelevel(memberIndex,stepIndex,:,:),8)
1160 end do
1161 end do
1162 else
1163 call gsv_getField(statevector,ptr4d_r4,varName_opt=varName)
1164 do stepIndex = 1, numStep
1165 do gsvLevIndex = 1, nLev
1166 ensVarLevIndex = gsvLevIndex + ens_getOffsetFromVarName(ens,varName)
1167 ptr4d_r4(:,:,gsvLevIndex,stepIndex) = &
1168 ens%allLev_r4(ensVarLevIndex)%onelevel(memberIndex,stepIndex,:,:)
1169 end do
1170 end do
1171 end if
1172 end if
1173 end do
1174
1175 end if
1176
1177 if (associated(varNamesInGsv)) deallocate(varNamesInGsv)
1178 if (associated(varNamesInEns)) deallocate(varNamesInEns)
1179
1180 end subroutine ens_copyMember
1181
1182 !--------------------------------------------------------------------------
1183 ! ens_insertMember
1184 !--------------------------------------------------------------------------
1185 subroutine ens_insertMember(ens, statevector, memberIndex)
1186 !
1187 !:Purpose: Copy the supplied stateVector in the selected ensemble member.
1188 !
1189 implicit none
1190
1191 ! Arguments:
1192 type(struct_ens), intent(inout) :: ens
1193 type(struct_gsv), intent(inout) :: statevector
1194 integer, intent(in) :: memberIndex
1195
1196 ! Locals:
1197 real(4), pointer :: ptr4d_r4(:,:,:,:)
1198 real(8), pointer :: ptr4d_r8(:,:,:,:)
1199 integer :: k1, k2, varLevIndex, stepIndex, numStep, varIndex
1200 integer :: gsvLevIndex, ensVarLevIndex, nLev
1201 character(len=4), pointer :: varNamesInEns(:)
1202 character(len=4), pointer :: varNamesInGsv(:)
1203 character(len=4) :: varName
1204 logical :: sameVariables
1205
1206 if (.not. gsv_isAllocated(ens%statevector_work)) then
1207 call utl_abort('ens_insertMember: ens not allocated')
1208 end if
1209
1210 numStep = ens%statevector_work%numStep
1211
1212 nullify(varNamesInEns)
1213 call gsv_varNamesList(varNamesInEns, ens%statevector_work)
1214 nullify(varNamesInGsv)
1215 call gsv_varNamesList(varNamesInGsv, statevector)
1216
1217 sameVariables = .false.
1218 if (size(ens%statevector_work%varExistlist) == size(statevector%varExistlist)) then
1219 if (all(ens%statevector_work%varExistlist == statevector%varExistlist)) then
1220 sameVariables = .true.
1221 end if
1222 end if
1223
1224 if (sameVariables) then
1225
1226 k1 = ens%statevector_work%mykBeg
1227 k2 = ens%statevector_work%mykEnd
1228
1229 if (ens%dataKind == 8) then
1230 call gsv_getField(statevector,ptr4d_r8)
1231 do stepIndex = 1, numStep
1232 do varLevIndex = k1, k2
1233 ens%allLev_r8(varLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
1234 ptr4d_r8(:,:,varLevIndex,stepIndex)
1235 end do
1236 end do
1237 else if (ens%dataKind == 4) then
1238 if (gsv_getDataKind(statevector) == 8) then
1239 call gsv_getField(statevector,ptr4d_r8)
1240 do stepIndex = 1, numStep
1241 do varLevIndex = k1, k2
1242 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
1243 real(ptr4d_r8(:,:,varLevIndex,stepIndex),4)
1244 end do
1245 end do
1246 else
1247 call gsv_getField(statevector,ptr4d_r4)
1248 do stepIndex = 1, numStep
1249 do varLevIndex = k1, k2
1250 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
1251 ptr4d_r4(:,:,varLevIndex,stepIndex)
1252 end do
1253 end do
1254 end if
1255 end if
1256
1257 else
1258
1259 do varIndex = 1, size(varNamesInGsv)
1260 varName = varNamesInGsv(varIndex)
1261 nLev = gsv_getNumLev(statevector,vnl_varLevelFromVarname(varName),varName)
1262 if (ens%dataKind == 8) then
1263 call gsv_getField(statevector,ptr4d_r8,varName_opt=varName)
1264 do stepIndex = 1, numStep
1265 do gsvLevIndex = 1, nLev
1266 ensVarLevIndex = gsvLevIndex + ens_getOffsetFromVarName(ens,varName)
1267 ens%allLev_r8(ensVarLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
1268 ptr4d_r8(:,:,gsvLevIndex,stepIndex)
1269 end do
1270 end do
1271 else if (ens%dataKind == 4) then
1272 if (gsv_getDataKind(statevector) == 8) then
1273 call gsv_getField(statevector,ptr4d_r8,varName_opt=varName)
1274 do stepIndex = 1, numStep
1275 do gsvLevIndex = 1, nLev
1276 ensVarLevIndex = gsvLevIndex + ens_getOffsetFromVarName(ens,varName)
1277 ens%allLev_r4(ensVarLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
1278 real(ptr4d_r8(:,:,gsvLevIndex,stepIndex),4)
1279 end do
1280 end do
1281 else
1282 call gsv_getField(statevector,ptr4d_r4,varName_opt=varName)
1283 do stepIndex = 1, numStep
1284 do gsvLevIndex = 1, nLev
1285 ensVarLevIndex = gsvLevIndex + ens_getOffsetFromVarName(ens,varName)
1286 ens%allLev_r4(ensVarLevIndex)%onelevel(memberIndex,stepIndex,:,:) = &
1287 ptr4d_r4(:,:,gsvLevIndex,stepIndex)
1288 end do
1289 end do
1290 end if
1291 end if
1292 end do ! varIndex
1293
1294 end if
1295
1296 if (associated(varNamesInGsv)) deallocate(varNamesInGsv)
1297 if (associated(varNamesInEns)) deallocate(varNamesInEns)
1298
1299 end subroutine ens_insertMember
1300
1301 !--------------------------------------------------------------------------
1302 ! ens_getMask
1303 !--------------------------------------------------------------------------
1304 subroutine ens_getMask(ens,oceanMask)
1305 !
1306 !:Purpose: Copy the instance of oceanMask from inside the ens object
1307 ! to the supplied instance of oceanMask.
1308 !
1309 implicit none
1310
1311 ! Arguments:
1312 type(struct_ens), intent(inout) :: ens
1313 type(struct_ocm), intent(inout) :: oceanMask
1314
1315 call ocm_copyMask(ens%statevector_work%oceanMask,oceanMask)
1316
1317 end subroutine ens_getMask
1318
1319 !--------------------------------------------------------------------------
1320 ! ens_copyMaskToGsv
1321 !--------------------------------------------------------------------------
1322 subroutine ens_copyMaskToGsv(ens,statevector)
1323 !
1324 !:Purpose: Copy the instance of oceanMask from inside the ens object
1325 ! to the instance inside the supplied stateVector object.
1326 !
1327 implicit none
1328
1329 ! Arguments:
1330 type(struct_ens), intent(inout) :: ens
1331 type(struct_gsv), intent(inout) :: statevector
1332
1333 call gsv_copyMask(ens%statevector_work,statevector)
1334
1335 end subroutine ens_copyMaskToGsv
1336
1337 !--------------------------------------------------------------------------
1338 ! ens_varExist
1339 !--------------------------------------------------------------------------
1340 function ens_varExist(ens,varName) result(varExist)
1341 !
1342 !:Purpose: Return true if the specified variable name exists in the ensemble.
1343 !
1344 implicit none
1345
1346 ! Arguments:
1347 type(struct_ens), intent(in) :: ens
1348 character(len=*), intent(in) :: varName
1349 ! Result:
1350 logical :: varExist
1351
1352 varExist = gsv_varExist(ens%statevector_work, varName)
1353
1354 end function ens_varExist
1355
1356 !--------------------------------------------------------------------------
1357 ! ens_varNamesList
1358 !--------------------------------------------------------------------------
1359 subroutine ens_varNamesList(varNames,ens_opt)
1360 !
1361 !:Purpose: Return a list of the variable names that exist in the ensemble.
1362 !
1363 implicit none
1364
1365 ! Arguments:
1366 type(struct_ens), optional, intent(in) :: ens_opt
1367 character(len=4), pointer, intent(inout) :: varNames(:)
1368
1369 if (associated(varNames)) then
1370 call utl_abort('ens_varNamesList: varNames must be NULL pointer on input')
1371 end if
1372
1373 if (present(ens_opt)) then
1374 call gsv_varNamesList(varNames, ens_opt%statevector_work)
1375 else
1376 call gsv_varNamesList(varNames)
1377 end if
1378
1379 end subroutine ens_varNamesList
1380
1381 !--------------------------------------------------------------------------
1382 ! ens_getNumLev
1383 !--------------------------------------------------------------------------
1384 function ens_getNumLev(ens,varLevel,varName_opt) result(nlev)
1385 !
1386 !:Purpose: Return the number of vertical levels of the ensemble.
1387 !
1388 implicit none
1389
1390 ! Arguments:
1391 type(struct_ens), intent(in) :: ens
1392 character(len=*), intent(in) :: varLevel
1393 character(len=*), optional, intent(in) :: varName_opt
1394 ! Result:
1395 integer :: nlev
1396
1397 nlev = vco_getNumLev(ens%statevector_work%vco,varLevel,varName_opt)
1398
1399 end function ens_getNumLev
1400
1401 !--------------------------------------------------------------------------
1402 ! ens_getNumMembers
1403 !--------------------------------------------------------------------------
1404 function ens_getNumMembers(ens) result(numMembers)
1405 !
1406 !:Purpose: Return the number of members in the ensemble.
1407 !
1408 implicit none
1409
1410 ! Arguments:
1411 type(struct_ens), intent(in) :: ens
1412 ! Result:
1413 integer :: numMembers
1414
1415 numMembers = ens%numMembers
1416
1417 end function ens_getNumMembers
1418
1419
1420 !--------------------------------------------------------------------------
1421 ! ens_getNumSubEns
1422 !--------------------------------------------------------------------------
1423 function ens_getNumSubEns(ens) result(numMembers)
1424 !
1425 !:Purpose: Return the number of sub-ensembles in the ensemble.
1426 !
1427 implicit none
1428
1429 ! Arguments:
1430 type(struct_ens), intent(in) :: ens
1431 ! Result:
1432 integer :: numMembers
1433
1434 numMembers = ens%numSubEns
1435
1436 end function ens_getNumSubEns
1437
1438 !--------------------------------------------------------------------------
1439 ! ens_getNumK
1440 !--------------------------------------------------------------------------
1441 function ens_getNumK(ens) result(numK)
1442 !
1443 !:Purpose: Return the number of kIndex (a.k.a. varLevs) values of the ensemble.
1444 !
1445 implicit none
1446
1447 ! Arguments:
1448 type(struct_ens), intent(in) :: ens
1449 ! Result:
1450 integer :: numK
1451
1452 numK = 1 + ens%statevector_work%mykEnd - ens%statevector_work%mykBeg
1453
1454 end function ens_getNumK
1455
1456 !--------------------------------------------------------------------------
1457 ! ens_getDataKind
1458 !--------------------------------------------------------------------------
1459 function ens_getDataKind(ens) result(dataKind)
1460 !
1461 !:Purpose: Return the floating point kind of the ensemble (4 or 8).
1462 !
1463 implicit none
1464
1465 ! Arguments:
1466 type(struct_ens), intent(in) :: ens
1467 ! Result:
1468 integer :: dataKind
1469
1470 dataKind = ens%dataKind
1471
1472 end function ens_getDataKind
1473
1474 !--------------------------------------------------------------------------
1475 ! ens_getPathName
1476 !--------------------------------------------------------------------------
1477 function ens_getPathName(ens) result(pathName)
1478 !
1479 !:Purpose: Return the path name for the ensemble files.
1480 !
1481 implicit none
1482
1483 ! Arguments:
1484 type(struct_ens), intent(in) :: ens
1485 ! Result:
1486 character(len=256) :: pathName
1487
1488 pathName = ens%ensPathName
1489
1490 end function ens_getPathName
1491
1492 !--------------------------------------------------------------------------
1493 ! ens_getOffsetFromVarName
1494 !--------------------------------------------------------------------------
1495 function ens_getOffsetFromVarName(ens,varName) result(offset)
1496 !
1497 !:Purpose: Return the offset of the kIndex for the specified variable name.
1498 !
1499 implicit none
1500
1501 ! Arguments:
1502 type(struct_ens), intent(in) :: ens
1503 character(len=*), intent(in) :: varName
1504 ! Result:
1505 integer :: offset
1506
1507 if (.not. ens_varExist(ens,varName)) then
1508 call utl_abort('ens_getOffsetFromVarName: this varName is not present in ens: '//trim(varName))
1509 end if
1510
1511 offset=gsv_getOffsetFromVarName(ens%statevector_work,varName)
1512
1513 end function ens_getOffsetFromVarName
1514
1515 !--------------------------------------------------------------------------
1516 ! ens_getLevFromK
1517 !--------------------------------------------------------------------------
1518 function ens_getLevFromK(ens,kIndex) result(levIndex)
1519 !
1520 !:Purpose: Return the level index from the kIndex value.
1521 !
1522 implicit none
1523
1524 ! Arguments:
1525 type(struct_ens), intent(in) :: ens
1526 integer, intent(in) :: kIndex
1527 ! Result:
1528 integer :: levIndex
1529
1530 levIndex = gsv_getLevFromK(ens%statevector_work,kIndex)
1531
1532 end function ens_getLevFromK
1533
1534 !--------------------------------------------------------------------------
1535 ! ens_getKFromLevVarName
1536 !--------------------------------------------------------------------------
1537 function ens_getKFromLevVarName(ens, levIndex, varName) result(kIndex)
1538 !
1539 !:Purpose: Return the kIndex value for the specified level index
1540 ! and variable name.
1541 !
1542 implicit none
1543
1544 ! Arguments:
1545 type(struct_ens), intent(in) :: ens
1546 integer, intent(in) :: levIndex
1547 character(len=*), intent(in) :: varName
1548 ! Result:
1549 integer :: kIndex
1550
1551 kIndex = levIndex + gsv_getOffsetFromVarName(ens%statevector_work,trim(varName))
1552
1553 end function ens_getKFromLevVarName
1554
1555 !--------------------------------------------------------------------------
1556 ! ens_getVarNameFromK
1557 !--------------------------------------------------------------------------
1558 function ens_getVarNameFromK(ens,kIndex) result(varName)
1559 !
1560 !:Purpose: Return the variable name from the specified kIndex value.
1561 !
1562 implicit none
1563
1564 ! Arguments:
1565 type(struct_ens), intent(in) :: ens
1566 integer, intent(in) :: kIndex
1567 ! Result:
1568 character(len=4) :: varName
1569
1570 varName = gsv_getVarNameFromK(ens%statevector_work,kIndex)
1571
1572 end function ens_getVarNameFromK
1573
1574 !--------------------------------------------------------------------------
1575 ! ens_getVco
1576 !--------------------------------------------------------------------------
1577 function ens_getVco(ens) result(vco_ptr)
1578 !
1579 !:Purpose: Return a pointer to the verticalCoord object associate with the ensemble.
1580 !
1581 implicit none
1582
1583 ! Arguments:
1584 type(struct_ens), intent(in) :: ens
1585 ! Result:
1586 type(struct_vco), pointer :: vco_ptr
1587
1588 vco_ptr => ens%statevector_work%vco
1589
1590 end function ens_getVco
1591
1592 !--------------------------------------------------------------------------
1593 ! ens_getHco
1594 !--------------------------------------------------------------------------
1595 function ens_getHco(ens) result(hco_ptr)
1596 !
1597 !:Purpose: Return a pointer to the horizontalCoord object associate with the ensemble.
1598 !
1599 implicit none
1600
1601 ! Arguments:
1602 type(struct_ens), intent(in) :: ens
1603 ! Result:
1604 type(struct_hco), pointer :: hco_ptr
1605
1606 hco_ptr => ens%statevector_work%hco
1607
1608 end function ens_getHco
1609
1610 !--------------------------------------------------------------------------
1611 ! ens_getLatLonBounds
1612 !--------------------------------------------------------------------------
1613 subroutine ens_getLatLonBounds(ens, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
1614 !
1615 !:Purpose: Return the longitude and latitude index bounds for this mpi task.
1616 !
1617 implicit none
1618
1619 ! Arguments:
1620 type(struct_ens), intent(in) :: ens
1621 integer, intent(out) :: myLonBeg
1622 integer, intent(out) :: myLonEnd
1623 integer, intent(out) :: myLatBeg
1624 integer, intent(out) :: myLatEnd
1625
1626 myLonBeg = ens%statevector_work%myLonBeg
1627 myLonEnd = ens%statevector_work%myLonEnd
1628 myLatBeg = ens%statevector_work%myLatBeg
1629 myLatEnd = ens%statevector_work%myLatEnd
1630
1631 end subroutine ens_getLatLonBounds
1632
1633 !--------------------------------------------------------------------------
1634 ! ens_getNumStep
1635 !--------------------------------------------------------------------------
1636 function ens_getNumStep(ens) result(numStep)
1637 !
1638 !:Purpose: Return the number of time steps stored in the ensemble.
1639 !
1640 implicit none
1641
1642 ! Arguments:
1643 type(struct_ens), intent(in) :: ens
1644 ! Result:
1645 integer :: numStep
1646
1647 numStep = ens%statevector_work%numStep
1648
1649 end function ens_getNumStep
1650
1651 !--------------------------------------------------------------------------
1652 ! ens_computeMean
1653 !--------------------------------------------------------------------------
1654 subroutine ens_computeMean(ens, computeSubEnsMeans_opt, numSubEns_opt)
1655 !
1656 !:Purpose: Internally compute the ensemble mean.
1657 !
1658 implicit none
1659
1660 ! Arguments:
1661 type(struct_ens), intent(inout) :: ens
1662 logical, optional, intent(in) :: computeSubEnsMeans_opt
1663 integer, optional, intent(out) :: numSubEns_opt
1664
1665 ! Locals:
1666 logical :: computeSubEnsMeans, lExists
1667 character(len=256), parameter :: subEnsIndexFileName = 'subEnsembleIndex.txt'
1668 integer :: kulin, ierr, memberIndex, memberIndex2, stepIndex, subEnsIndex
1669 integer :: k1, k2, varLevIndex, lon1, lon2, lat1, lat2, numStep, lonIndex, latIndex
1670 integer :: fnom, fclos
1671
1672 if (present(computeSubEnsMeans_opt)) then
1673 computeSubEnsMeans = computeSubEnsMeans_opt
1674 else
1675 computeSubEnsMeans = .false.
1676 end if
1677
1678 ! Read sub-ensemble index list from file, if it exists
1679 if (.not. allocated(ens%subEnsIndexList)) then
1680 allocate(ens%subEnsIndexList(ens%numMembers))
1681 end if
1682 if ( computeSubEnsMeans ) then
1683 write(*,*) 'ens_computeMean: checking in ensemble directory if file with sub-ensemble index list exists: ',subEnsIndexFileName
1684 inquire(file=trim(ens%enspathname) // trim(subEnsIndexFileName),exist=lExists)
1685 if ( lExists ) then
1686 kulin = 0
1687 ierr = fnom(kulin,trim(ens%enspathname) // trim(subEnsIndexFileName),'FMT+SEQ+R/O',0)
1688 do memberIndex = 1, ens%numMembers
1689 read(kulin,*) memberIndex2, ens%subEnsIndexList(memberIndex)
1690 write(*,*) 'read from sub-ensemble index list: ',memberIndex, memberIndex2, ens%subEnsIndexList(memberIndex)
1691 end do
1692 ierr = fclos(kulin)
1693 else
1694 call utl_abort('ens_computeMean: could not find file with sub-ensemble index list')
1695 end if
1696 else
1697 ens%subEnsIndexList(:) = 1
1698 end if
1699 ens%numSubEns = maxval(ens%subEnsIndexList(:))
1700 if (.not. allocated(ens%nEnsSubEns)) then
1701 allocate(ens%nEnsSubEns(ens%numSubEns))
1702 end if
1703 ens%nEnsSubEns(:) = 0
1704 do memberIndex = 1, ens%numMembers
1705 ens%nEnsSubEns(ens%subEnsIndexList(memberIndex)) = ens%nEnsSubEns(ens%subEnsIndexList(memberIndex)) + 1
1706 end do
1707 write(*,*) 'ens_computeMean: number of sub-ensembles = ', ens%numSubEns
1708 write(*,*) 'ens_computeMean: number of members in each sub-ensemble = ', ens%nensSubEns(:)
1709
1710 lon1 = ens%statevector_work%myLonBeg
1711 lon2 = ens%statevector_work%myLonEnd
1712 lat1 = ens%statevector_work%myLatBeg
1713 lat2 = ens%statevector_work%myLatEnd
1714 k1 = ens%statevector_work%mykBeg
1715 k2 = ens%statevector_work%mykEnd
1716 numStep = ens%statevector_work%numStep
1717
1718 if (.not. allocated(ens%allLev_ensMean_r8)) then
1719 call ens_allocateMean(ens)
1720 else
1721 do varLevIndex = k1, k2
1722 ens%allLev_ensMean_r8(varLevIndex)%onelevel(:,:,:,:) = 0.0d0
1723 end do
1724 end if
1725 ens%meanIsComputed = .true.
1726
1727 ! Compute ensemble mean(s)
1728 !$OMP PARALLEL DO PRIVATE (varLevIndex,latIndex,lonIndex,stepIndex,memberIndex,subEnsIndex)
1729 do varLevIndex = k1, k2
1730 do latIndex = lat1, lat2
1731 do lonIndex = lon1, lon2
1732 do stepIndex = 1, ens%statevector_work%numStep
1733 do memberIndex = 1, ens%numMembers
1734 ens%allLev_ensMean_r8(varLevIndex)%onelevel(ens%subEnsIndexList(memberIndex),stepIndex,lonIndex,latIndex) = &
1735 ens%allLev_ensMean_r8(varLevIndex)%onelevel(ens%subEnsIndexList(memberIndex),stepIndex,lonIndex,latIndex) + &
1736 dble(ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex))
1737 end do
1738 do subEnsIndex = 1, ens%numSubEns
1739 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,lonIndex,latIndex) = &
1740 ens%allLev_ensMean_r8(varLevIndex)%onelevel(subEnsIndex,stepIndex,lonIndex,latIndex) / &
1741 dble(ens%nEnsSubEns(subEnsIndex))
1742 end do
1743 end do
1744 end do
1745 end do
1746 end do
1747 !$OMP END PARALLEL DO
1748
1749 ! provide output argument value
1750 if ( present(numSubEns_opt) ) numSubEns_opt = ens%numSubEns
1751
1752 end subroutine ens_computeMean
1753
1754 !--------------------------------------------------------------------------
1755 ! ens_computeStdDev
1756 !--------------------------------------------------------------------------
1757 subroutine ens_computeStdDev(ens, containsScaledPerts_opt)
1758 !
1759 !:Purpose: Internally compute the ensemble stdDev.
1760 !
1761 implicit none
1762
1763 ! Arguments:
1764 type(struct_ens), intent(inout) :: ens
1765 logical, optional, intent(in) :: containsScaledPerts_opt
1766
1767 ! Locals:
1768 integer :: memberIndex, stepIndex, subEnsIndex
1769 integer :: k1, k2, varLevIndex, lon1, lon2, lat1, lat2, numStep, lonIndex, latIndex
1770 real(8), allocatable :: subEnsStdDev(:)
1771 logical :: containsScaledPerts
1772
1773 if ( present(containsScaledPerts_opt) ) then
1774 containsScaledPerts = containsScaledPerts_opt
1775 else
1776 containsScaledPerts = .false.
1777 if (.not.ens%meanIsRemoved .and. .not.ens%meanIsComputed) then
1778 if (mmpi_myid == 0) write(*,*) 'ens_computeStdDev: compute Mean since it was not already done'
1779 call ens_computeMean( ens )
1780 end if
1781 end if
1782
1783 ! Read sub-ensemble index list from file, if it exists
1784 ! The sub-ensembles should have been already read in routine 'ens_computeMean'
1785 write(*,*) 'ens_computeStdDev: number of sub-ensembles = ', ens%numSubEns
1786 write(*,*) 'ens_computeStdDev: number of members in each sub-ensemble = ', ens%nensSubEns(:)
1787
1788 lon1 = ens%statevector_work%myLonBeg
1789 lon2 = ens%statevector_work%myLonEnd
1790 lat1 = ens%statevector_work%myLatBeg
1791 lat2 = ens%statevector_work%myLatEnd
1792 k1 = ens%statevector_work%mykBeg
1793 k2 = ens%statevector_work%mykEnd
1794 numStep = ens%statevector_work%numStep
1795
1796 if (.not. allocated(ens%allLev_ensStdDev_r8)) then
1797 call ens_allocateStdDev(ens)
1798 else
1799 do varLevIndex = k1, k2
1800 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(:,:,:,:) = 0.0d0
1801 end do
1802 end if
1803 ens%StdDevIsComputed = .true.
1804
1805 if (containsScaledPerts) then
1806
1807 if (ens%numSubEns /= 1) then
1808 call utl_abort('ens_computeStdDev: sub-ensemble approach not compatible with scale perturbations')
1809 end if
1810
1811 ! Compute the ensemble StdDev from previously scale ensemble perturbations
1812 ! (i.e. pert = (fcst-mean)/(nEns-1) )
1813
1814 !$OMP PARALLEL DO PRIVATE (varLevIndex,latIndex,lonIndex,stepIndex,memberIndex)
1815 do varLevIndex = k1, k2
1816 do latIndex = lat1, lat2
1817 do lonIndex = lon1, lon2
1818 do stepIndex = 1, ens%statevector_work%numStep
1819
1820 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) = 0.d0
1821
1822 do memberIndex = 1, ens%numMembers
1823 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) = &
1824 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) + &
1825 dble(ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex))**2
1826 end do
1827
1828 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) = &
1829 sqrt(ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex))
1830
1831 end do
1832 end do
1833 end do
1834 end do
1835 !$OMP END PARALLEL DO
1836
1837 else
1838
1839 allocate(subEnsStdDev(ens%numSubEns))
1840
1841 ! Compute global ensemble StdDev as the sqrt of the mean of each suchensemble variance
1842 ! var_subens(i) = sum( ( ens(j) - mean_subens(i) )**2, j=1..numEns_subens(i) ) / ( numEns_subens(i) - 1 )
1843 ! var_allensensemble = sum( numEns_subens(i) * var_subens(i), i=1..numSubEns)
1844 ! stddev = sqrt( var_allensensemble / numEnsTotal )
1845
1846 !$OMP PARALLEL DO PRIVATE (varLevIndex,latIndex,lonIndex,stepIndex,memberIndex,subEnsIndex,subEnsStdDev)
1847 do varLevIndex = k1, k2
1848 do latIndex = lat1, lat2
1849 do lonIndex = lon1, lon2
1850 do stepIndex = 1, ens%statevector_work%numStep
1851
1852 subEnsStdDev(:) = 0.0d0
1853
1854 if (ens%meanIsRemoved) then
1855 do memberIndex = 1, ens%numMembers
1856 subEnsStdDev(ens%subEnsIndexList(memberIndex)) = &
1857 subEnsStdDev(ens%subEnsIndexList(memberIndex)) + &
1858 (dble(ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex)))**2
1859 end do
1860 else
1861 do memberIndex = 1, ens%numMembers
1862 subEnsStdDev(ens%subEnsIndexList(memberIndex)) = &
1863 subEnsStdDev(ens%subEnsIndexList(memberIndex)) + &
1864 (dble(ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex)) - &
1865 ens%allLev_ensMean_r8(varLevIndex)%onelevel(ens%subEnsIndexList(memberIndex),stepIndex,lonIndex,latIndex))**2
1866 end do
1867 end if
1868
1869 do subEnsIndex = 1, ens%numSubEns
1870 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) = &
1871 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) + &
1872 ens%nEnsSubEns(subEnsIndex)*subEnsStdDev(subEnsIndex)/(ens%nEnsSubEns(subEnsIndex)-1)
1873
1874 end do
1875
1876 ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) = &
1877 sqrt( ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex) / dble(ens%numMembers) )
1878
1879 end do
1880 end do
1881 end do
1882 end do
1883 !$OMP END PARALLEL DO
1884
1885 deallocate(subEnsStdDev)
1886
1887 end if
1888
1889 end subroutine ens_computeStdDev
1890
1891 !--------------------------------------------------------------------------
1892 ! ens_normalize
1893 !--------------------------------------------------------------------------
1894 subroutine ens_normalize(ens)
1895 !
1896 !:Purpose: Normalize the ensemble by the 3D ensemble stdDev.
1897 !
1898 implicit none
1899
1900 ! Arguments:
1901 type(struct_ens), intent(inout) :: ens
1902
1903 ! Locals:
1904 integer :: lon1, lon2, lat1, lat2, k1, k2, numStep
1905 integer :: varLevIndex, latIndex, lonIndex, stepIndex, memberIndex
1906 real(8) :: factor
1907
1908 if (.not. ens%StdDevIsComputed) then
1909 if (mmpi_myid == 0) write(*,*) 'ens_normalize: compute Std. Dev. since it was not already done'
1910 call ens_computeStdDev( ens )
1911 end if
1912
1913 lon1 = ens%statevector_work%myLonBeg
1914 lon2 = ens%statevector_work%myLonEnd
1915 lat1 = ens%statevector_work%myLatBeg
1916 lat2 = ens%statevector_work%myLatEnd
1917 k1 = ens%statevector_work%mykBeg
1918 k2 = ens%statevector_work%mykEnd
1919 numStep = ens%statevector_work%numStep
1920
1921 !$OMP PARALLEL DO PRIVATE (varLevIndex,latIndex,lonIndex,stepIndex,memberIndex,factor)
1922 do varLevIndex = k1, k2
1923 do latIndex = lat1, lat2
1924 do lonIndex = lon1, lon2
1925 do stepIndex = 1, numStep
1926 do memberIndex = 1, ens%numMembers
1927
1928 if (ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(ens%subEnsIndexList(memberIndex),stepIndex,lonIndex,latIndex) > 0.0d0 ) then
1929 factor = 1.0d0/ens%allLev_ensStdDev_r8(varLevIndex)%onelevel(ens%subEnsIndexList(memberIndex),stepIndex,lonIndex,latIndex)
1930 else
1931 factor = 0.0d0
1932 endif
1933
1934 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
1935 real( real(ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex),8) * factor, 4)
1936 end do
1937 end do
1938 end do
1939 end do
1940 end do
1941 !$OMP END PARALLEL DO
1942
1943 end subroutine ens_normalize
1944
1945 !--------------------------------------------------------------------------
1946 ! ens_removeMean
1947 !--------------------------------------------------------------------------
1948 subroutine ens_removeMean(ens)
1949 !
1950 !:Purpose: Subtract the ensemble mean from each member.
1951 !
1952 implicit none
1953
1954 ! Arguments:
1955 type(struct_ens), intent(inout) :: ens
1956
1957 ! Locals:
1958 integer :: lon1, lon2, lat1, lat2, k1, k2, numStep
1959 integer :: varLevIndex, latIndex, lonIndex, stepIndex, memberIndex
1960
1961 lon1 = ens%statevector_work%myLonBeg
1962 lon2 = ens%statevector_work%myLonEnd
1963 lat1 = ens%statevector_work%myLatBeg
1964 lat2 = ens%statevector_work%myLatEnd
1965 k1 = ens%statevector_work%mykBeg
1966 k2 = ens%statevector_work%mykEnd
1967 numStep = ens%statevector_work%numStep
1968
1969 !$OMP PARALLEL DO PRIVATE (varLevIndex,latIndex,lonIndex,stepIndex,memberIndex)
1970 do varLevIndex = k1, k2
1971 do latIndex = lat1, lat2
1972 do lonIndex = lon1, lon2
1973 do stepIndex = 1, numStep
1974 do memberIndex = 1, ens%numMembers
1975 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
1976 real( (real(ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex),8) - &
1977 ens%allLev_ensMean_r8(varLevIndex)%onelevel(ens%subEnsIndexList(memberIndex),stepIndex,lonIndex,latIndex)), 4 )
1978 end do
1979 end do
1980 end do
1981 end do
1982 end do
1983 !$OMP END PARALLEL DO
1984
1985 ens%meanIsRemoved = .true.
1986
1987 end subroutine ens_removeMean
1988
1989 !--------------------------------------------------------------------------
1990 ! ens_removeGlobalMean
1991 !--------------------------------------------------------------------------
1992 subroutine ens_removeGlobalMean(ens)
1993 !
1994 !:Purpose: Subtract the 2D global mean from each member.
1995 !
1996 implicit none
1997
1998 ! Arguments:
1999 type(struct_ens), intent(inout) :: ens
2000
2001 ! Locals:
2002 integer :: lon1, lon2, lat1, lat2, k1, k2, numStep, ierr
2003 integer :: kIndex, latIndex, lonIndex, stepIndex, memberIndex
2004 real(8) :: globalMean, globalMean_mpiglobal
2005
2006 if ( .not. ens%statevector_work%hco%global ) then
2007 call utl_abort('ens_removeGlobalMean: must never by applied to limited-area ensembles')
2008 end if
2009
2010 lon1 = ens%statevector_work%myLonBeg
2011 lon2 = ens%statevector_work%myLonEnd
2012 lat1 = ens%statevector_work%myLatBeg
2013 lat2 = ens%statevector_work%myLatEnd
2014 k1 = ens%statevector_work%mykBeg
2015 k2 = ens%statevector_work%mykEnd
2016 numStep = ens%statevector_work%numStep
2017
2018 do kIndex = k1, k2
2019 do memberIndex = 1, ens%numMembers
2020 do stepIndex = 1, numStep
2021
2022 ! Compute the domain mean
2023 globalMean = 0.d0
2024 do latIndex = lat1, lat2
2025 do lonIndex = lon1, lon2
2026 globalMean = globalMean + &
2027 real(ens%allLev_r4(kIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex),8)
2028 end do
2029 end do
2030
2031 call rpn_comm_allreduce(globalMean, globalMean_mpiglobal,1,&
2032 "mpi_double_precision","mpi_sum","GRID",ierr)
2033 globalMean_mpiglobal = globalMean_mpiglobal / &
2034 (real(ens%statevector_work%ni,8)*real(ens%statevector_work%nj,8))
2035
2036 ! Remove it
2037 do latIndex = lat1, lat2
2038 do lonIndex = lon1, lon2
2039 ens%allLev_r4(kIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
2040 ens%allLev_r4(kIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) - real(globalMean_mpiglobal,4)
2041 end do
2042 end do
2043
2044 end do
2045 end do
2046 end do
2047
2048 end subroutine ens_removeGlobalMean
2049
2050 !--------------------------------------------------------------------------
2051 ! ens_recenter
2052 !--------------------------------------------------------------------------
2053 subroutine ens_recenter(ens, recenteringMean, recenteringCoeff_opt, &
2054 recenteringCoeffLand_opt, recenteringCoeffScalar_opt, &
2055 alternativeEnsembleMean_opt, &
2056 ensembleControlMember_opt, scaleFactor_opt, &
2057 numMembersToRecenter_opt)
2058 !
2059 !:Purpose:
2060 ! To compute:
2061 ! ..math::
2062 ! x_recentered =
2063 ! scaleFactor*x_original
2064 ! + recenteringCoeff*( x_recenteringMean
2065 ! - scaleFactor*x_ensembleMean
2066 ! )
2067 implicit none
2068
2069 ! Arguments:
2070 type(struct_ens), intent(inout) :: ens
2071 type(struct_gsv), intent(in) :: recenteringMean
2072 type(struct_gsv), optional, intent(in) :: alternativeEnsembleMean_opt
2073 type(struct_gsv), optional, intent(in) :: ensembleControlMember_opt
2074 real(8), optional, intent(in) :: recenteringCoeff_opt(:,:)
2075 real(8), optional, intent(in) :: recenteringCoeffLand_opt(:)
2076 real(8), optional, intent(in) :: recenteringCoeffScalar_opt
2077 real(8), optional, intent(in) :: scaleFactor_opt(:)
2078 integer, optional, intent(in) :: numMembersToRecenter_opt
2079
2080 ! Locals:
2081 real(8), pointer :: ptr4d_r8(:,:,:,:)
2082 real(8), pointer :: alternativeEnsembleMean_r8(:,:,:,:)
2083 real(8), pointer :: ptr4d_ensembleControlmember_r8(:,:,:,:)
2084 real(8) :: scaleFactor(vco_maxNumLevels)
2085 real(8) :: recenteringCoeffArray(vco_maxNumLevels,ens%numMembers)
2086 real(8) :: recenteringCoeffArrayLand(ens%numMembers)
2087 real(8) :: recenteringCoeffArrayUsed(ens%numMembers)
2088 real(8) :: increment, thisScaleFactor
2089 real(4), pointer :: ptr4d_r4(:,:,:,:)
2090 integer :: lon1, lon2, lat1, lat2, k1, k2, numStep, numMembersToRecenter
2091 integer :: varLevIndex, latIndex, lonIndex, stepIndex, memberIndex, levIndex
2092 character(len=4) :: varLevel
2093 character(len=2) :: varKind
2094
2095 ! if an alternative mean is not provided, we need to ensure ens mean is present
2096 if ( .not. present(alternativeEnsembleMean_opt)) then
2097 if ( .not. ens%meanIsComputed ) then
2098 if (mmpi_myid == 0) write(*,*) 'ens_recenter: compute Mean since it was not already done'
2099 call ens_computeMean( ens )
2100 end if
2101 end if
2102
2103 if ( present(recenteringCoeff_opt) ) then
2104 recenteringCoeffArray(:,:) = recenteringCoeff_opt(:,:)
2105 else if ( present(recenteringCoeffScalar_opt) ) then
2106 recenteringCoeffArray(:,:) = recenteringCoeffScalar_opt
2107 else
2108 call utl_abort('ens_recenter: Must specify recenteringCoeff_opt or recenteringCoeffScalar_opt')
2109 end if
2110
2111 if ( present(scaleFactor_opt) ) then
2112 ! scaleFactor cannot be used at the same time as a recenteringCoeff different from 1.0
2113 if ( any (abs(recenteringCoeffArray(:,:) - 1.0D0) > 1.0D-5) ) then
2114 call utl_abort('ens_recenter: recenteringCoeff must be equal to 1.0 when using scaleFactor')
2115 end if
2116 scaleFactor = scaleFactor_opt
2117 else
2118 scaleFactor(:) = 1.0D0
2119 end if
2120
2121 if ( present(recenteringCoeffLand_opt) ) then
2122 if (any(recenteringCoeffLand_opt < 0.0D0)) then
2123 ! negative coeff specified for land, apply same coeff as other variables
2124 recenteringCoeffArrayLand(:) = recenteringCoeffArray(max(1,ens_getNumLev(ens,'MM')),:)
2125 else
2126 ! specified coeff for land variables used for all members
2127 write(*,*) 'ens_recenter: different recentering applied to land variables:', &
2128 recenteringCoeffLand_opt(:)
2129 recenteringCoeffArrayLand(:) = recenteringCoeffLand_opt(:)
2130 end if
2131 else
2132 ! coeff for land not specified, apply same coeff as other variables
2133 recenteringCoeffArrayLand(:) = recenteringCoeffArray(max(1,ens_getNumLev(ens,'MM')),:)
2134 end if
2135
2136 if (present(numMembersToRecenter_opt)) then
2137 numMembersToRecenter = numMembersToRecenter_opt
2138 else
2139 numMembersToRecenter = ens%numMembers
2140 end if
2141
2142 lon1 = ens%statevector_work%myLonBeg
2143 lon2 = ens%statevector_work%myLonEnd
2144 lat1 = ens%statevector_work%myLatBeg
2145 lat2 = ens%statevector_work%myLatEnd
2146 k1 = ens%statevector_work%mykBeg
2147 k2 = ens%statevector_work%mykEnd
2148 numStep = ens%statevector_work%numStep
2149
2150 nullify(ptr4d_r4, ptr4d_r8)
2151 if (gsv_getDataKind(recenteringMean) == 8) then
2152 call gsv_getField(recenteringMean,ptr4d_r8)
2153 else
2154 call gsv_getField(recenteringMean,ptr4d_r4)
2155 end if
2156 if(present(alternativeEnsembleMean_opt)) then
2157 call gsv_getField(alternativeEnsembleMean_opt,alternativeEnsembleMean_r8)
2158 else
2159 nullify(alternativeEnsembleMean_r8)
2160 end if
2161
2162 if (present(ensembleControlMember_opt)) then
2163 call gsv_getField(ensembleControlMember_opt,ptr4d_ensembleControlmember_r8)
2164 else
2165 nullify(ptr4d_ensembleControlmember_r8)
2166 end if
2167
2168 !$OMP PARALLEL DO PRIVATE(varLevIndex,varLevel,varKind,levIndex,thisScaleFactor), &
2169 !$OMP PRIVATE(latIndex,lonIndex,stepIndex,memberIndex,increment,recenteringCoeffArrayUsed)
2170 do varLevIndex = k1, k2
2171
2172 ! define scaling factor as a function of vertical level and variable type
2173 varLevel = vnl_varLevelFromVarname(ens_getVarNameFromK(ens, varLevIndex))
2174 if ( trim(varLevel) == 'SF' .or. trim(varLevel) == 'SFMM' .or. trim(varLevel) == 'SFTH' ) then
2175 ! use lowest momentum level for surface variables
2176 levIndex = ens_getNumLev(ens, 'MM')
2177 else if ( (trim(varLevel) == 'MM') .and. (ens%statevector_work%vco%Vcode == 5002) ) then
2178 levIndex = ens_getLevFromK(ens, varLevIndex) + 1
2179 else
2180 levIndex = ens_getLevFromK(ens, varLevIndex)
2181 end if
2182 thisScaleFactor = scaleFactor(levIndex)
2183
2184 ! determine which recentering coeff are used: general or land-specific
2185 varKind = vnl_varKindFromVarname(ens_getVarNameFromK(ens, varLevIndex))
2186 if ( varKind == 'LD' ) then
2187 recenteringCoeffArrayUsed(:) = recenteringCoeffArrayLand(:)
2188 else
2189 recenteringCoeffArrayUsed(:) = recenteringCoeffArray(levIndex,:)
2190 end if
2191
2192 do latIndex = lat1, lat2
2193 do lonIndex = lon1, lon2
2194 do stepIndex = 1, numStep
2195 if(present(alternativeEnsembleMean_opt)) then
2196 if (associated(ptr4d_r8)) then
2197 increment = ptr4d_r8(lonIndex,latIndex,varLevIndex,stepIndex) - &
2198 thisScaleFactor * &
2199 alternativeEnsembleMean_r8(lonIndex,latIndex,varLevIndex,stepIndex)
2200 else
2201 increment = real(ptr4d_r4(lonIndex,latIndex,varLevIndex,stepIndex),8) - &
2202 thisScaleFactor * &
2203 alternativeEnsembleMean_r8(lonIndex,latIndex,varLevIndex,stepIndex)
2204 end if
2205 else
2206 if (associated(ptr4d_r8)) then
2207 increment = ptr4d_r8(lonIndex,latIndex,varLevIndex,stepIndex) - &
2208 thisScaleFactor * &
2209 ens%allLev_ensMean_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex)
2210 else
2211 increment = real(ptr4d_r4(lonIndex,latIndex,varLevIndex,stepIndex),8) - &
2212 thisScaleFactor * &
2213 ens%allLev_ensMean_r8(varLevIndex)%onelevel(1,stepIndex,lonIndex,latIndex)
2214 end if
2215 end if
2216 if (present(ensembleControlMember_opt)) then
2217 ptr4d_ensembleControlMember_r8(lonIndex,latIndex,varLevIndex,stepIndex) = &
2218 thisScaleFactor * &
2219 ptr4d_ensembleControlMember_r8(lonIndex,latIndex,varLevIndex,stepIndex) + &
2220 recenteringCoeffArrayUsed(1)*increment
2221 else
2222 do memberIndex = 1, numMembersToRecenter
2223 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex) = &
2224 real( real(thisScaleFactor * &
2225 ens%allLev_r4(varLevIndex)%onelevel(memberIndex,stepIndex,lonIndex,latIndex),8) + &
2226 recenteringCoeffArrayUsed(memberIndex)*increment, 4)
2227 end do
2228 end if
2229 end do
2230 end do
2231 end do
2232 end do
2233 !$OMP END PARALLEL DO
2234
2235 end subroutine ens_recenter
2236
2237 !--------------------------------------------------------------------------
2238 ! ens_readEnsemble
2239 !--------------------------------------------------------------------------
2240 subroutine ens_readEnsemble(ens, ensPathName, biPeriodic, &
2241 vco_file_opt, varNames_opt, checkModelTop_opt, &
2242 containsFullField_opt, ignoreDate_opt)
2243 !
2244 !:Purpose: Read the ensemble from disk in parallel and do mpi communication
2245 ! so that all members for a given lat-lon tile are present on each
2246 ! mpi task.
2247 !
2248 implicit none
2249
2250 ! Arguments:
2251 type(struct_ens), intent(inout) :: ens
2252 character(len=*), intent(in) :: ensPathName
2253 logical, intent(in) :: biPeriodic
2254 character(len=*), optional, intent(in) :: varNames_opt(:)
2255 type(struct_vco), pointer, optional, intent(in) :: vco_file_opt
2256 logical, optional, intent(in) :: checkModelTop_opt
2257 logical, optional, intent(in) :: containsFullField_opt
2258 logical, optional, intent(in) :: ignoreDate_opt
2259
2260 ! Locals:
2261 type(struct_gsv) :: statevector_file_r4, statevector_hint_r4, statevector_member_r4
2262 type(struct_hco), pointer :: hco_file, hco_ens, hco_coregrid
2263 type(struct_vco), pointer :: vco_file, vco_ens
2264 real(4), allocatable :: gd_send_r4(:,:,:,:)
2265 real(4), allocatable :: gd_recv_r4(:,:,:,:)
2266 real(4), pointer :: ptr3d_r4(:,:,:)
2267 integer,pointer :: dateStampList(:)
2268 integer :: batchIndex, nsize, ierr
2269 integer :: yourid, youridx, youridy
2270 integer, allocatable :: readFilePE(:), memberIndexFromMemberStep(:), stepIndexFromMemberStep(:)
2271 integer, allocatable :: batchIndexFromMemberStep(:)
2272 integer :: sendsizes(mmpi_nprocs), recvsizes(mmpi_nprocs), senddispls(mmpi_nprocs), recvdispls(mmpi_nprocs)
2273 integer :: lonPerPEmax, latPerPEmax, ni, nj, numK, numStep, numMembers, numLevelsToSend2
2274 integer :: memberIndex, memberIndex2, stepIndex, stepIndex2, procIndex, memberStepIndex, memberStepIndex2
2275 integer :: kIndexBeg, kIndexEnd, kCount, memberStepIndexStart, lastReadFilePE
2276 character(len=256) :: ensFileName
2277 character(len=2) :: typvar
2278 character(len=12) :: etiket
2279 character(len=4), pointer :: anlVar(:)
2280 logical :: thisProcIsAsender(mmpi_nprocs), doMpiCommunication
2281 logical :: verticalInterpNeeded, horizontalInterpNeeded, horizontalPaddingNeeded
2282 logical :: checkModelTop, containsFullField, ignoreDate
2283 character(len=4), pointer :: varNames(:)
2284 integer, parameter :: numLevelsToSend = 10
2285
2286 write(*,*) 'ens_readEnsemble: starting'
2287 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2288
2289 if ( .not. ens%allocated ) then
2290 call utl_abort('ens_readEnsemble: ensemble object not allocated!')
2291 end if
2292
2293 !
2294 !- 1. Initial setup
2295 !
2296 lonPerPEmax = ens%statevector_work%lonPerPEmax
2297 latPerPEmax = ens%statevector_work%latPerPEmax
2298 ni = ens%statevector_work%ni
2299 nj = ens%statevector_work%nj
2300 numK = ens%statevector_work%nk
2301 numStep = ens%statevector_work%numStep
2302 numMembers = ens%numMembers
2303
2304 dateStampList => ens%statevector_work%dateStampList
2305
2306 ens%ensPathName = trim(ensPathName)
2307
2308 ! Determine which MPI tasks read which members/steps to minimize file copies to ram disk
2309 allocate(batchIndexFromMemberStep(numMembers*numStep))
2310 allocate(readFilePE(numMembers*numStep))
2311 allocate(stepIndexFromMemberStep(numMembers*numStep))
2312 allocate(memberIndexFromMemberStep(numMembers*numStep))
2313 do stepIndex = 1, numStep
2314 do memberIndex = 1, numMembers
2315 memberStepIndex = ((stepIndex-1)*numMembers) + memberIndex
2316 stepIndexFromMemberStep(memberStepIndex) = stepIndex
2317 memberIndexFromMemberStep(memberStepIndex) = memberIndex
2318
2319 if (memberStepIndex == 1) then
2320 ! Very first member/step
2321 readFilePE(memberStepIndex) = 0
2322 batchIndexFromMemberStep(memberStepIndex) = 1
2323 else
2324 ! Increment MPI task ID and keep same batch
2325 readFilePE(memberStepIndex) = readFilePE(memberStepIndex-1) + 1
2326 batchIndexFromMemberStep(memberStepIndex) = batchIndexFromMemberStep(memberStepIndex-1)
2327 end if
2328
2329 ! Decide if we need to move to the next batch
2330 if (memberIndex == 1) then
2331 if (readFilePE(memberStepIndex) == 0) then
2332 ! First MPI task reading member 1, try to fit all others in this batch
2333 lastReadFilePE = numMembers - 1
2334 else
2335 ! Check if we can fit this full time step in this batch
2336 if (lastReadFilePE + numMembers < mmpi_nprocs) then
2337 lastReadFilePE = lastReadFilePE + numMembers
2338 end if
2339 ! If numMembers > nprocs, move to next batch
2340 if (numMembers > mmpi_nprocs) then
2341 readFilePE(memberStepIndex) = 0
2342 batchIndexFromMemberStep(memberStepIndex) = batchIndexFromMemberStep(memberStepIndex-1)+ 1
2343 lastReadFilePE = numMembers - 1
2344 end if
2345 end if
2346 ! Ensure we limit ourselves to the total number of MPI tasks
2347 lastReadFilePE = min(lastreadFilePE, mmpi_nprocs - 1)
2348 end if
2349 ! Move to next batch if we reached lastReadFilePE
2350 if (readFilePE(memberStepIndex) == lastReadFilePE + 1) then
2351 readFilePE(memberStepIndex) = 0
2352 batchIndexFromMemberStep(memberStepIndex) = batchIndexFromMemberStep(memberStepIndex-1)+ 1
2353 lastReadFilePE = min(numMembers - memberIndex, mmpi_nprocs - 1)
2354 end if
2355
2356 if (mmpi_myid == 0) then
2357 write(*,*) 'ens_readEnsemble: batchIndex, memberIndex, stepIndex, memberStepIndex, readFilePE = ', &
2358 batchIndexFromMemberStep(memberStepIndex), memberIndex, stepIndex, memberStepIndex, &
2359 readFilePE(memberStepIndex)
2360 end if
2361
2362 end do
2363 end do
2364
2365 ! the default is to check whether output grid has a higher top than input grid during vertical interpolation
2366 if ( present(checkModelTop_opt) ) then
2367 checkModelTop = checkModelTop_opt
2368 else
2369 checkModelTop = .true.
2370 end if
2371
2372 ! the default is to NOT ignore the date - can only ignore if numStep == 1
2373 if ( present(ignoreDate_opt) ) then
2374 ignoreDate = ignoreDate_opt
2375 else
2376 ignoreDate = .false.
2377 end if
2378 if ( ignoreDate .and. (numStep > 1) ) then
2379 call utl_abort('ens_readEnsemble: cannot ignore date if numStep > 1')
2380 end if
2381
2382 ! Set up hco and vco for ensemble files
2383 call fln_ensFileName(ensFileName, ensPathName, memberIndex_opt=1, copyToRamDisk_opt=.false., &
2384 fileMemberIndex1_opt=ens%fileMemberIndex1)
2385
2386 nullify(anlVar)
2387 call gsv_varNamesList(anlVar)
2388 nullify(hco_file)
2389 call hco_SetupFromFile(hco_file, ensFileName, ' ', 'ENSFILEGRID', varName_opt=anlVar(1))
2390 if ( present(vco_file_opt) ) then
2391 ! use the input vertical grid provided
2392 vco_file => vco_file_opt
2393 else
2394 ! find the info from the ensemble files
2395 nullify(vco_file)
2396 if ( mmpi_myid == 0 ) then
2397 call vco_SetupFromFile(vco_file, ensFileName)
2398 end if
2399 call vco_mpiBcast(vco_file)
2400 end if
2401 hco_ens => gsv_getHco(ens%statevector_work)
2402 hco_coregrid => ens%hco_core
2403 vco_ens => gsv_getVco(ens%statevector_work)
2404 horizontalInterpNeeded = (.not. hco_equal(hco_ens, hco_file))
2405 verticalInterpNeeded = (.not. vco_equal(vco_ens, vco_file))
2406
2407 ! More efficient handling of common case where input is on Z grid, analysis in on G grid
2408 if ( hco_file%grtyp == 'Z' .and. hco_ens%grtyp == 'G' ) then
2409 if ( hco_file%ni == (hco_ens%ni+1) ) then
2410 write(*,*) 'ens_readEnsemble: no interpolation done for equivalent Gaussian grid stored as a Z grid'
2411 horizontalInterpNeeded = .false.
2412 end if
2413 end if
2414 ! In limited-area mode, avoid horizontal interpolation when the ensemble is on the coregrid
2415 horizontalPaddingNeeded = .false.
2416 if ( .not. hco_file%global ) then
2417 if ( hco_file%ni == hco_coregrid%ni .and. hco_file%nj == hco_coregrid%nj ) then
2418 if (mmpi_myid == 0) then
2419 write(*,*) 'ens_readEnsemble: no interpolation needed for ensemble on the limited-area coregrid'
2420 end if
2421 horizontalInterpNeeded = .false.
2422 if ( hco_file%ni /= hco_ens%ni .or. hco_file%nj /= hco_ens%nj ) then
2423 if (mmpi_myid == 0) then
2424 write(*,*) 'ens_readEnsemble: horizontal padding needed for limited-area ensemble'
2425 end if
2426 horizontalPaddingNeeded = .true.
2427 end if
2428 end if
2429 end if
2430
2431 if (mmpi_myid == 0) then
2432 write(*,*)
2433 write(*,*) 'ens_readEnsemble: dateStampList=',dateStampList(1:numStep)
2434 write(*,*)
2435 if (horizontalInterpNeeded ) write(*,*) 'ens_readEnsemble: HORIZONTAL interpolation is needed'
2436 if (verticalInterpNeeded ) write(*,*) 'ens_readEnsemble: VERTICAL interpolation is needed'
2437 if (horizontalPaddingNeeded) write(*,*) 'ens_readEnsemble: HORIZONTAL padding is needed'
2438 end if
2439
2440 ! Input type
2441 if (present(containsFullField_opt)) then
2442 containsFullField = containsFullField_opt
2443 else
2444 containsFullField = .true.
2445 end if
2446 if (mmpi_myid == 0) then
2447 write(*,*)
2448 write(*,*) 'ens_readEnsemble: containsFullField = ', containsFullField
2449 end if
2450
2451 !
2452 !- 2. Ensemble forecasts reading loop
2453 !
2454
2455 nullify(varNames)
2456 if (present(varNames_opt)) then
2457 allocate(varNames(size(varNames_opt)))
2458 varNames(:) = varNames_opt(:)
2459 else
2460 call gsv_varNamesList(varNames)
2461 end if
2462
2463 !- 2.1 Loop on time, ensemble member, variable, level
2464 memberStepIndexStart = 1
2465 stepLoop: do stepIndex = 1, numStep
2466 write(*,*) ' '
2467 write(*,*) 'ens_readEnsemble: starting to read time level ', stepIndex
2468
2469 memberLoop: do memberIndex = 1, numMembers
2470
2471 memberStepIndex = ((stepIndex-1)*numMembers) + memberIndex
2472 batchIndex = batchIndexFromMemberStep(memberStepIndex)
2473 if (mmpi_myid == readFilePE(memberStepIndex)) then
2474
2475 ! allocate the needed statevector objects
2476 call gsv_allocate(statevector_member_r4, 1, hco_ens, vco_ens, &
2477 datestamp_opt = dateStampList(stepIndex), mpi_local_opt = .false., &
2478 varNames_opt = varNames, dataKind_opt = 4, &
2479 hInterpolateDegree_opt = ens%hInterpolateDegree)
2480 if (horizontalInterpNeeded .or. verticalInterpNeeded .or. horizontalPaddingNeeded) then
2481 call gsv_allocate(statevector_file_r4, 1, hco_file, vco_file, &
2482 datestamp_opt = dateStampList(stepIndex), mpi_local_opt = .false., &
2483 varNames_opt = varNames, dataKind_opt = 4, &
2484 hInterpolateDegree_opt = ens%hInterpolateDegree)
2485 end if
2486 if (verticalInterpNeeded) then
2487 call gsv_allocate(statevector_hint_r4, 1, hco_ens, vco_file, &
2488 datestamp_opt = dateStampList(stepIndex), mpi_local_opt = .false., &
2489 varNames_opt = varNames, dataKind_opt = 4, &
2490 hInterpolateDegree_opt = ens%hInterpolateDegree)
2491 end if
2492 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2493
2494 ! Read the file
2495 call fln_ensFileName(ensFileName, ensPathName, memberIndex_opt=memberIndex, &
2496 fileMemberIndex1_opt=ens%fileMemberIndex1)
2497 typvar = ' '
2498 etiket = ' '
2499 if (.not. horizontalInterpNeeded .and. &
2500 .not. verticalInterpNeeded .and. &
2501 .not. horizontalPaddingNeeded ) then
2502 call gio_readFile(statevector_member_r4, ensFileName, etiket, typvar, &
2503 containsFullField, ignoreDate_opt=ignoreDate)
2504 else
2505 call gio_readFile(statevector_file_r4, ensFileName, etiket, typvar, &
2506 containsFullField, ignoreDate_opt=ignoreDate)
2507 end if
2508
2509 ! Remove file from ram disk if no longer needed
2510 if ( all(readFilePE(memberStepIndex+1:numMembers*numStep) /= mmpi_myid) .or. &
2511 (stepIndex == numStep) .or. &
2512 (batchIndex == maxval(batchIndexFromMemberStep(:))) ) then
2513 ierr = ram_remove(ensFileName)
2514 end if
2515
2516 ! do any required interpolation
2517 if (horizontalInterpNeeded .and. verticalInterpNeeded) then
2518 call int_hInterp_gsv(statevector_file_r4, statevector_hint_r4)
2519 call int_vInterp_gsv( statevector_hint_r4, statevector_member_r4, &
2520 Ps_in_hPa_opt=.true.,checkModelTop_opt=checkModelTop)
2521
2522 else if (horizontalInterpNeeded .and. .not. verticalInterpNeeded) then
2523 call int_hInterp_gsv(statevector_file_r4, statevector_member_r4)
2524
2525 else if (.not. horizontalInterpNeeded .and. verticalInterpNeeded) then
2526 if (horizontalPaddingNeeded) then
2527 call gsv_hPad(statevector_file_r4, statevector_hint_r4)
2528 else
2529 call gsv_copy(statevector_file_r4, statevector_hint_r4)
2530 end if
2531 call int_vInterp_gsv( statevector_hint_r4, statevector_member_r4, &
2532 Ps_in_hPa_opt=.true.,checkModelTop_opt=checkModelTop)
2533
2534 else if (horizontalPaddingNeeded) then
2535 call gsv_hPad(statevector_file_r4, statevector_member_r4)
2536 end if
2537
2538 ! unit conversion
2539 call gio_fileUnitsToStateUnits(statevector_member_r4, containsFullField)
2540
2541 ! Create bi-periodic forecasts when using scale-dependent localization in LAM mode
2542 if ( .not. hco_ens%global .and. biperiodic ) then
2543 call gsv_getField(statevector_member_r4,ptr3d_r4)
2544 call lgt_mach_r4(ptr3d_r4, & ! INOUT
2545 ni, nj, statevector_member_r4%nk) ! IN
2546 end if
2547
2548 ! copy over some time related and other parameters
2549 ens%statevector_work%deet = statevector_member_r4%deet
2550 ens%statevector_work%dateOriginList(stepIndex) = statevector_member_r4%dateOriginList(1)
2551 ens%statevector_work%npasList(stepIndex) = statevector_member_r4%npasList(1)
2552 ens%statevector_work%ip2List(stepIndex) = statevector_member_r4%ip2List(1)
2553 ens%statevector_work%etiket = statevector_member_r4%etiket
2554 ens%statevector_work%onPhysicsGrid(:) = statevector_member_r4%onPhysicsGrid(:)
2555 ens%statevector_work%hco_physics => statevector_member_r4%hco_physics
2556 ! if it exists, copy over mask from member read on task 0, which should always read
2557 if(mmpi_myid == 0) then
2558 call gsv_copyMask(stateVector_member_r4, ens%stateVector_work)
2559 end if
2560
2561 end if ! locally read one member
2562
2563 ! MPI communication: from 1 ensemble member per process to 1 lat-lon tile per process
2564 if (memberStepIndex == numStep*numMembers) then
2565 ! last member/step was read, do last communication
2566 doMpiCommunication = .true.
2567 else if (batchIndex < batchIndexFromMemberStep(memberStepIndex+1)) then
2568 ! next member/step is in next batch, do communication
2569 doMpiCommunication = .true.
2570 else
2571 ! do not do communication, still reading members/steps
2572 doMpiCommunication = .false.
2573 end if
2574
2575 if (doMpiCommunication) then
2576 write(*,*) 'ens_readEnsemble: Do communication for batchIndex = ', batchIndex
2577 write(*,*) ' for the memberStepIndex range = ', memberStepIndexStart, memberStepIndex
2578
2579 ! determine which tasks have something to send and let everyone know
2580 do procIndex = 1, mmpi_nprocs
2581 thisProcIsAsender(procIndex) = .false.
2582 if ( mmpi_myid == (procIndex-1) .and. gsv_isAllocated(stateVector_member_r4) ) then
2583 thisProcIsAsender(procIndex) = .true.
2584 end if
2585 call rpn_comm_bcast(thisProcIsAsender(procIndex), 1, &
2586 'MPI_LOGICAL', procIndex-1, 'GRID', ierr)
2587 end do
2588
2589 do kIndexBeg = 1, numK, numLevelsToSend
2590 kIndexEnd = min(numK,kIndexBeg+numLevelsToSend-1)
2591 numLevelsToSend2 = kIndexEnd - kIndexBeg + 1
2592
2593 ! prepare for alltoallv
2594 nsize = lonPerPEmax * latPerPEmax * numLevelsToSend2
2595
2596 ! only send the data from tasks with data, same amount to all
2597 sendsizes(:) = 0
2598 if ( gsv_isAllocated(stateVector_member_r4) ) then
2599 do procIndex = 1, mmpi_nprocs
2600 sendsizes(procIndex) = nsize
2601 end do
2602 end if
2603 senddispls(1) = 0
2604 do procIndex = 2, mmpi_nprocs
2605 senddispls(procIndex) = senddispls(procIndex-1) + sendsizes(procIndex-1)
2606 end do
2607
2608 ! all tasks recv only from those with data
2609 recvsizes(:) = 0
2610 do procIndex = 1, mmpi_nprocs
2611 if ( thisProcIsAsender(procIndex) ) then
2612 recvsizes(procIndex) = nsize
2613 end if
2614 end do
2615 recvdispls(1) = 0
2616 do procIndex = 2, mmpi_nprocs
2617 recvdispls(procIndex) = recvdispls(procIndex-1) + recvsizes(procIndex-1)
2618 end do
2619
2620 if (gsv_isAllocated(statevector_member_r4)) then
2621 allocate(gd_send_r4(lonPerPEmax,latPerPEmax,numLevelsToSend2,mmpi_nprocs))
2622 gd_send_r4(:,:,:,:) = 0.0
2623 call gsv_getField(statevector_member_r4,ptr3d_r4)
2624 !$OMP PARALLEL DO PRIVATE(youridy,youridx,yourid)
2625 do youridy = 0, (mmpi_npey-1)
2626 do youridx = 0, (mmpi_npex-1)
2627 yourid = youridx + youridy*mmpi_npex
2628 gd_send_r4(1:ens%statevector_work%allLonPerPE(youridx+1), &
2629 1:ens%statevector_work%allLatPerPE(youridy+1), :, yourid+1) = &
2630 ptr3d_r4(ens%statevector_work%allLonBeg(youridx+1):ens%statevector_work%allLonEnd(youridx+1), &
2631 ens%statevector_work%allLatBeg(youridy+1):ens%statevector_work%allLatEnd(youridy+1), kIndexBeg:kIndexEnd)
2632 end do
2633 end do
2634 !$OMP END PARALLEL DO
2635 else
2636 allocate(gd_send_r4(1,1,1,1))
2637 gd_send_r4(:,:,:,:) = 0.0
2638 end if
2639 allocate(gd_recv_r4(lonPerPEmax,latPerPEmax,numLevelsToSend2,max(mmpi_nprocs,numMembers)))
2640 gd_recv_r4(:,:,:,:) = 0.0
2641
2642 if (mmpi_nprocs.gt.1) then
2643 call mpi_alltoallv(gd_send_r4, sendsizes, senddispls, mmpi_datyp_real4, &
2644 gd_recv_r4, recvsizes, recvdispls, mmpi_datyp_real4, &
2645 mmpi_comm_grid, ierr)
2646 else
2647 gd_recv_r4(:,:,:,1) = gd_send_r4(:,:,:,1)
2648 end if
2649
2650 !$OMP PARALLEL DO PRIVATE(kCount,memberStepIndex2,memberIndex2,stepIndex2,yourid)
2651 do kCount = 1, numLevelsToSend2
2652 do memberStepIndex2 = memberStepIndexStart, memberStepIndex
2653 memberIndex2 = memberIndexFromMemberStep(memberStepIndex2)
2654 stepIndex2 = stepIndexFromMemberStep(memberStepIndex2)
2655 yourid = readFilePE(memberStepIndex2)
2656 ens%allLev_r4(kCount+kIndexBeg-1)%onelevel(memberIndex2,stepIndex2, :, :) = &
2657 gd_recv_r4(1:ens%statevector_work%lonPerPE, 1:ens%statevector_work%latPerPE, &
2658 kCount, yourid+1)
2659 end do
2660 end do
2661 !$OMP END PARALLEL DO
2662
2663 deallocate(gd_send_r4)
2664 deallocate(gd_recv_r4)
2665
2666 end do ! kIndexBeg
2667
2668 ! deallocate the needed statevector objects
2669 if (gsv_isAllocated(statevector_member_r4)) then
2670 call gsv_deallocate(statevector_member_r4)
2671 end if
2672 if (gsv_isAllocated(statevector_file_r4)) then
2673 call gsv_deallocate(statevector_file_r4)
2674 end if
2675 if (gsv_isAllocated(statevector_hint_r4)) then
2676 call gsv_deallocate(statevector_hint_r4)
2677 end if
2678
2679 memberStepIndexStart = memberStepIndex + 1
2680
2681 end if ! MPI communication
2682
2683 end do memberLoop
2684
2685 end do stepLoop
2686
2687 call gsv_communicateTimeParams(ens%statevector_work)
2688 call ocm_communicateMask(ens%statevector_work%oceanMask)
2689
2690 deallocate(datestamplist)
2691 call hco_deallocate(hco_file)
2692 if ( .not. present(vco_file_opt) ) then
2693 call vco_deallocate(vco_file)
2694 end if
2695 deallocate(readFilePE)
2696 deallocate(batchIndexFromMemberStep)
2697 deallocate(stepIndexFromMemberStep)
2698 deallocate(memberIndexFromMemberStep)
2699
2700 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2701 write(*,*) 'ens_readEnsemble: finished reading and communicating ensemble members...'
2702
2703 end subroutine ens_readEnsemble
2704
2705 !--------------------------------------------------------------------------
2706 ! ens_writeEnsemble
2707 !--------------------------------------------------------------------------
2708 subroutine ens_writeEnsemble(ens, ensPathName, ensFileNamePrefix, &
2709 etiket, typvar, &
2710 etiketAppendMemberNumber_opt, varNames_opt, &
2711 ip3_opt, containsFullField_opt, numBits_opt, &
2712 resetTimeParams_opt)
2713 !
2714 !:Purpose: Write the ensemble to disk by doing mpi transpose so that
2715 ! each mpi task can write a single member in parallel.
2716 !
2717 implicit none
2718
2719 ! Arguments:
2720 type(struct_ens), intent(inout) :: ens
2721 character(len=*), intent(in) :: ensPathName
2722 character(len=*), intent(in) :: ensFileNamePrefix
2723 character(len=*), intent(in) :: etiket
2724 character(len=*), intent(in) :: typvar
2725 character(len=*), optional, intent(in) :: varNames_opt(:)
2726 integer, optional, intent(in) :: ip3_opt
2727 integer, optional, intent(in) :: numBits_opt
2728 logical, optional, intent(in) :: etiketAppendMemberNumber_opt
2729 logical, optional, intent(in) :: containsFullField_opt
2730 logical, optional, intent(in) :: resetTimeParams_opt
2731
2732 ! Locals:
2733 type(struct_gsv) :: statevector_member_r4
2734 type(struct_hco), pointer :: hco_ens
2735 type(struct_vco), pointer :: vco_ens
2736 real(4), allocatable :: gd_send_r4(:,:,:,:)
2737 real(4), allocatable :: gd_recv_r4(:,:,:,:)
2738 real(4), pointer :: ptr3d_r4(:,:,:)
2739 integer, allocatable :: dateStampList(:)
2740 integer :: batchIndex, nsize, ierr
2741 integer :: yourid, youridx, youridy
2742 integer :: writeFilePE(1000)
2743 integer :: lonPerPE, lonPerPEmax, latPerPE, latPerPEmax, ni, nj
2744 integer :: numK, numStep, numlevelstosend, numlevelstosend2
2745 integer :: memberIndex, memberIndex2, stepIndex, kIndexBeg, kIndexEnd, kCount
2746 integer :: ip3, ensFileExtLength, maximumBaseEtiketLength
2747 character(len=256) :: ensFileName
2748 character(len=12) :: etiketStr ! this is the etiket that will be used to write files
2749 !! The two next declarations are sufficient until we reach 10^10 members
2750 character(len=10) :: memberIndexStr ! this is the member number in a character string
2751 character(len=10) :: ensFileExtLengthStr ! this is a string containing the same number as 'ensFileExtLength'
2752 character(len=4), pointer :: varNamesInEns(:)
2753 logical :: containsFullField
2754
2755 write(*,*) 'ens_writeEnsemble: starting'
2756 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2757
2758 if ( .not. ens%allocated ) then
2759 call utl_abort('ens_writeEnsemble: ensemble object not allocated!')
2760 end if
2761
2762 !- 1. Initial setup
2763
2764 nullify(varNamesInEns)
2765 if (present(varNames_opt)) then
2766 allocate(varNamesInEns(size(varNames_opt)))
2767 varNamesInEns(:) = varNames_opt(:)
2768 else
2769 call gsv_varNamesList(varNamesInEns, ens%statevector_work)
2770 end if
2771
2772 if (present(ip3_opt)) then
2773 ip3 = ip3_opt
2774 else
2775 ip3 = 0
2776 end if
2777
2778 if (present(resetTimeParams_opt)) then
2779 if (resetTimeParams_opt) then
2780 call gsv_resetTimeParams(ens%statevector_work)
2781 end if
2782 end if
2783
2784 lonPerPE = ens%statevector_work%lonPerPE
2785 latPerPE = ens%statevector_work%latPerPE
2786 lonPerPEmax = ens%statevector_work%lonPerPEmax
2787 latPerPEmax = ens%statevector_work%latPerPEmax
2788 ni = ens%statevector_work%ni
2789 nj = ens%statevector_work%nj
2790 numK = ens%statevector_work%nk
2791 numStep = ens%statevector_work%numStep
2792
2793 ens%ensPathName = trim(ensPathName)
2794
2795 ! Memory allocation
2796 numLevelsToSend = 10
2797 allocate(gd_send_r4(lonPerPEmax,latPerPEmax,numLevelsToSend,mmpi_nprocs))
2798 allocate(gd_recv_r4(lonPerPEmax,latPerPEmax,numLevelsToSend,mmpi_nprocs))
2799 gd_send_r4(:,:,:,:) = 0.0
2800 gd_recv_r4(:,:,:,:) = 0.0
2801
2802 allocate(dateStampList(numStep))
2803 call tim_getstamplist(dateStampList,numStep,tim_getDatestamp())
2804
2805 do memberIndex = 1, ens%numMembers
2806 writeFilePE(memberIndex) = mod(memberIndex-1,mmpi_nprocs)
2807 end do
2808
2809 hco_ens => gsv_getHco(ens%statevector_work)
2810 vco_ens => gsv_getVco(ens%statevector_work)
2811
2812 if (mmpi_myid == 0) then
2813 write(*,*)
2814 write(*,*) 'ens_writeEnsemble: dateStampList=',dateStampList(1:numStep)
2815 write(*,*)
2816 end if
2817
2818 !
2819 !- 2. Ensemble forecasts writing loop
2820 !
2821
2822 !- 2.1 Loop on time, ensemble member, variable, level
2823 do stepIndex = 1, numStep
2824 write(*,*) ' '
2825 write(*,*) 'ens_writeEnsemble: starting to write time level ', stepIndex
2826
2827 ! allocate the needed statevector objects
2828 call gsv_allocate(statevector_member_r4, 1, hco_ens, vco_ens, &
2829 datestamp_opt=dateStampList(stepIndex), mpi_local_opt=.false., &
2830 varNames_opt=varNamesInEns, dataKind_opt=4)
2831
2832 ! copy over some time related parameters
2833 statevector_member_r4%deet = ens%statevector_work%deet
2834 statevector_member_r4%dateOriginList(1) = ens%statevector_work%dateOriginList(stepIndex)
2835 statevector_member_r4%npasList(1) = ens%statevector_work%npasList(stepIndex)
2836 statevector_member_r4%ip2List(1) = ens%statevector_work%ip2List(stepIndex)
2837 ! if it exists, copy over mask from work statevector to member being written
2838 call gsv_copyMask(ens%stateVector_work, stateVector_member_r4)
2839
2840 do memberIndex = 1, ens%numMembers
2841
2842 ! MPI communication: from 1 lat-lon tile per process to 1 ensemble member per process
2843 if (writeFilePE(memberIndex) == 0) then
2844
2845 batchIndex = ceiling(dble(memberIndex + mmpi_nprocs - 1)/dble(mmpi_nprocs))
2846
2847 do kIndexBeg = 1, numK, numLevelsToSend
2848 kIndexEnd = min(numK,kIndexBeg+numLevelsToSend-1)
2849 numLevelsToSend2 = kIndexEnd - kIndexBeg + 1
2850
2851 if ( ens%dataKind == 8 ) then
2852 !$OMP PARALLEL DO PRIVATE(kCount,memberIndex2,yourid)
2853 do kCount = 1, numLevelsToSend2
2854 do memberIndex2 = 1+(batchIndex-1)*mmpi_nprocs, min(ens%numMembers, batchIndex*mmpi_nprocs)
2855 yourid = writeFilePE(memberIndex2)
2856 gd_send_r4(1:lonPerPE,1:latPerPE,kCount,yourid+1) = &
2857 real(ens%allLev_r8(kCount+kIndexBeg-1)%onelevel(memberIndex2,stepIndex,:,:),4)
2858 end do
2859 end do
2860 !$OMP END PARALLEL DO
2861 else
2862 !$OMP PARALLEL DO PRIVATE(kCount,memberIndex2,yourid)
2863 do kCount = 1, numLevelsToSend2
2864 do memberIndex2 = 1+(batchIndex-1)*mmpi_nprocs, min(ens%numMembers, batchIndex*mmpi_nprocs)
2865 yourid = writeFilePE(memberIndex2)
2866 gd_send_r4(1:lonPerPE,1:latPerPE,kCount,yourid+1) = &
2867 ens%allLev_r4(kCount+kIndexBeg-1)%onelevel(memberIndex2,stepIndex,:,:)
2868 end do
2869 end do
2870 !$OMP END PARALLEL DO
2871 end if
2872
2873 nsize = lonPerPEmax * latPerPEmax * numLevelsToSend2
2874 if (mmpi_nprocs > 1) then
2875 call rpn_comm_alltoall(gd_send_r4(:,:,1:numLevelsToSend2,:),nsize,"mpi_real4", &
2876 gd_recv_r4(:,:,1:numLevelsToSend2,:),nsize,"mpi_real4","GRID",ierr)
2877 else
2878 gd_recv_r4(:,:,1:numLevelsToSend2,1) = gd_send_r4(:,:,1:numLevelsToSend2,1)
2879 end if
2880
2881 call gsv_getField(statevector_member_r4,ptr3d_r4)
2882 !$OMP PARALLEL DO PRIVATE(youridy,youridx,yourid)
2883 do youridy = 0, (mmpi_npey-1)
2884 do youridx = 0, (mmpi_npex-1)
2885 yourid = youridx + youridy*mmpi_npex
2886 ptr3d_r4(ens%statevector_work%allLonBeg(youridx+1):ens%statevector_work%allLonEnd(youridx+1), &
2887 ens%statevector_work%allLatBeg(youridy+1):ens%statevector_work%allLatEnd(youridy+1), kIndexBeg:kIndexEnd) = &
2888 gd_recv_r4(1:ens%statevector_work%allLonPerPE(youridx+1), &
2889 1:ens%statevector_work%allLatPerPE(youridy+1), 1:numLevelsToSend2, yourid+1)
2890
2891 end do
2892 end do
2893 !$OMP END PARALLEL DO
2894
2895 end do ! kIndexBeg
2896
2897 end if ! MPI communication
2898
2899
2900 ! Write statevector to file
2901 if (mmpi_myid == writeFilePE(memberIndex)) then
2902
2903 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2904
2905 if ( typvar == 'A' .or. typvar == 'R' ) then
2906 if ( typvar == 'R' ) then
2907 call fln_ensAnlFileName( ensFileName, ensPathName, tim_getDateStamp(), &
2908 memberIndex_opt=memberIndex, &
2909 ensFileNamePrefix_opt=ensFileNamePrefix, &
2910 ensFileNameSuffix_opt='inc' )
2911 else
2912 call fln_ensAnlFileName( ensFileName, ensPathName, tim_getDateStamp(), &
2913 memberIndex_opt=memberIndex, &
2914 ensFileNamePrefix_opt=ensFileNamePrefix )
2915 end if
2916 ensFileExtLength = 4
2917 else
2918 call fln_ensFileName( ensFileName, ensPathName, memberIndex_opt=memberIndex, &
2919 ensFileNamePrefix_opt=ensFileNamePrefix, &
2920 shouldExist_opt=.false., ensembleFileExtLength_opt=ensFileExtLength, &
2921 fileMemberIndex1_opt=ens%fileMemberIndex1 )
2922 end if
2923
2924 ! Determine if ensemble is full fields (if yes, will be converted from K to C)
2925 if (present(containsFullField_opt)) then
2926 containsFullField = containsFullField_opt
2927 else
2928 containsFullField = (.not. ens%meanIsRemoved)
2929 end if
2930
2931 etiketStr = etiket
2932
2933 if (present(etiketAppendMemberNumber_opt)) then
2934 if (etiketAppendMemberNumber_opt .and. etiketStr /= 'UNDEFINED') then
2935 write(ensFileExtLengthStr,"(I1)") ensFileExtLength
2936 write(memberIndexStr,'(I0.' // trim(ensFileExtLengthStr) // ')') memberIndex
2937 ! 12 is the maximum length of an etiket for RPN fstd files
2938 maximumBaseEtiketLength = 12 - ensFileExtLength
2939 if ( len(trim(etiketStr)) >= maximumBaseEtiketLength ) then
2940 etiketStr = etiketStr(1:maximumBaseEtiketLength) // trim(memberIndexStr)
2941 else
2942 etiketStr = trim(etiketStr) // trim(memberIndexStr)
2943 end if
2944 end if
2945 end if
2946
2947 ! The routine 'gio_writeToFile' ignores the supplied
2948 ! argument for the etiket, here 'etiketStr', if
2949 ! 'statevector_member_r4%etiket' is different from
2950 ! 'UNDEFINED'. So we must define it explicitely in the
2951 ! 'statevector_member_r4'.
2952 statevector_member_r4%etiket = etiketStr
2953
2954 call gio_writeToFile( statevector_member_r4, ensFileName, etiketStr, ip3_opt = ip3, &
2955 typvar_opt = typvar, numBits_opt = numBits_opt, &
2956 containsFullField_opt = containsFullField )
2957
2958 end if ! locally written one member
2959
2960 end do ! memberIndex
2961
2962 ! deallocate the needed statevector objects
2963 call gsv_deallocate(statevector_member_r4)
2964
2965 end do ! time
2966
2967 deallocate(varNamesInEns)
2968 deallocate(gd_send_r4)
2969 deallocate(gd_recv_r4)
2970 deallocate(datestamplist)
2971
2972 write(*,*) 'Memory Used: ',get_max_rss()/1024,'Mb'
2973 write(*,*) 'ens_writeEnsemble: finished communicating and writing ensemble members...'
2974
2975 end subroutine ens_writeEnsemble
2976
2977 !--------------------------------------------------------------------------
2978 ! ens_applyMaskLAM
2979 !--------------------------------------------------------------------------
2980 subroutine ens_applyMaskLAM(ensIncrement, stateVectorAnalIncMask)
2981 !:Purpose: To apply a mask to an ensemble state vector for LAM grid
2982 !
2983 implicit none
2984
2985 ! Arguments
2986 type(struct_ens), intent(inout) :: ensIncrement
2987 type(struct_gsv), intent(in) :: stateVectorAnalIncMask
2988
2989 ! Locals
2990 real(4), pointer :: increment_ptr(:,:,:,:)
2991 real(pre_incrReal), pointer :: analIncMask_ptr(:,:,:)
2992 integer :: nEns, numVarLev, myLonBeg, myLonEnd, myLatBeg, myLatEnd
2993 integer :: varLevIndex, lonIndex, latIndex, stepIndex, memberIndex
2994
2995 write(*,*) 'ens_applyMaskLAM: starting'
2996
2997 if (.not.(ens_isAllocated(ensIncrement).and.(gsv_isAllocated(stateVectorAnalIncMask)))) then
2998 call utl_abort('epp_applyMaskLAM: increment and mask must be avaliable.')
2999 end if
3000
3001 call gsv_getField(stateVectorAnalIncMask, analIncMask_ptr)
3002
3003 nEns = ens_getNumMembers(ensIncrement)
3004 numVarLev = ens_getNumK(ensIncrement)
3005 call ens_getLatLonBounds(ensIncrement, myLonBeg, myLonEnd, myLatBeg, myLatEnd)
3006 do varLevIndex = 1, numVarLev
3007 increment_ptr => ens_getOneLev_r4(ensIncrement,varLevIndex)
3008 !$OMP PARALLEL DO PRIVATE (stepIndex,latIndex,lonIndex,memberIndex)
3009 do latIndex = myLatBeg, myLatEnd
3010 do lonIndex = myLonBeg, myLonEnd
3011 do stepIndex = 1, tim_nstepobsinc
3012 do memberIndex = 1, nEns
3013 increment_ptr(memberIndex,stepIndex,lonIndex,latIndex) = &
3014 increment_ptr(memberIndex,stepIndex,lonIndex,latIndex) * &
3015 analIncMask_ptr(lonIndex,latIndex,1)
3016 end do
3017 end do
3018 end do
3019 end do
3020 !$OMP END PARALLEL DO
3021 end do
3022 write(*,*) 'ens_applyMaskLAM: finished to mask each member of increments'
3023
3024 end subroutine ens_applyMaskLAM
3025
3026end module ensembleStateVector_mod
