1program midas_prepcma
2 !
3 !:Purpose: Read the observation files (usually after output by the background check)
4 ! and apply further quality control and thinning for use by the ``LETKF`` program
5 !
6 ! ---
7 !
8 !:Algorithm: After reading the observation files that have been processed by the background check,
9 ! the ``prepcma`` program rejects more observations based on the
10 ! various flags and conditions. It also performs further thinnings
11 ! on the data types of aircraft (AI), scatterometer (SC),
12 ! satellite winds (SW) and some radiance (TO). The rejection and thinning
13 ! are controlled by the options in the namelist of ``NAMPREPCMA``.
14 !
15 ! ---
16 !
17 !:File I/O: The required input files and produced output files are listed as follows.
18 !
19 ! --
20 !
21 !============================================== ==============================================================
22 ! Input and Output Files Description of file
23 !============================================== ==============================================================
24 ! ``flnml`` In - Main namelist file with parameters user may modify
25 ! ``flnml_static`` In - The "static" namelist that should not be modified
26 ! ``obserr`` In - Observation error statistics
27 ! ``obsfiles_$FAM/obs$FAM_$NNNN_$NNNN`` In - Observation file for each "family"
28 ! ``stats_tovs`` In - Satellite radiance observation errors
29 ! ``rtcoef_$PLATFORM_$SENSOR.dat`` In - RTTOV coefficient files
30 ! ``obsfiles_$FAM.updated/obs$FAM_$NNNN_$NNNN`` Out - final observation file for each family
31 !============================================== ==============================================================
32 !
33 ! ---
34 !
35 !:Synopsis: Below is a summary of the ``prepcma`` program calling sequence:
36 !
37 ! - **Initial setups:**
38 !
39 ! - Read the NAMPREPCMA namelist and check/modify some values.
40 !
41 ! - ``filt_setup``: set up list of elements to be assimilated and flags for rejection
42 !
43 ! - ``obsf_setup``: get observation file names and datestamp
44 !
45 ! - **Computation:**
46 !
47 ! - ``obsf_readFiles``: get the observations
48 !
49 ! - ``filt_suprep``: select the elements to assimilate and apply rejection flags
50 !
51 ! - ``oer_setObsErrors``: initialize obs error covariances and set flag
52 !
53 ! - ``oti_setup``: reject any observations outside the data assimilation window
54 !
55 ! - ``enkf_rejectHighLatIR``: reject all IR radiance observation in arctic and antarctic
56 !
57 ! - ``enkf_modifyAmsubObsError``: modify the obs error stddev for AMSUB in the tropics
58 !
59 ! - ``thinning_fam``: perform thinning for aircraft (AI), scatterometer (SC),
60 ! satellite winds (SW) and some radiance (TO)
61 !
62 ! - **Final steps:**
63 !
64 ! - ``obsf_writeFiles``: write to burp/sqlite files
65 !
66 ! - ``obsf_printFiles``: print to ascci file and to unformatted files
67 !
68 ! ---
69 !
70 !:Options: `List of namelist blocks <../namelists_in_each_program.html#prepcma>`_
71 ! that can affect the ``prepcma`` program.
72 !
73 ! * Some of the relevant namelist blocks used to configure the
74 ! prepcma are listed in the following table:
75 !
76 ! --
77 !
78 !=================== ====================== ===========================================
79 ! Program/Module Namelist Description of what is controlled
80 !=================== ====================== ===========================================
81 ! ``midas_prepcma`` ``NAMPREPCMA`` parameters for CMA format and others
82 ! to modify, reject and thinning some
83 ! observation data
84 ! ``timeCoord_mod`` ``NAMTIME`` assimilation time window length, temporal
85 ! resolution of the background state and the
86 ! analysis
87 ! ``tovs_nl_mod`` ``NAMTOV`` The list of satellite and instrument
88 !=================== ====================== ===========================================
89 !
90 !
91 use version_mod
92 use obsSpaceData_mod
93 use obsFiles_mod
94 use obsFilter_mod
95 use obsTimeInterp_mod
96 use obsErrors_mod
97 use tovsNL_mod
98 use timeCoord_mod
99 use enkf_mod
100 use utilities_mod
101 use midasMpi_mod
102 use ramDisk_mod
103 use regions_mod
104 use burpRead_mod
105
106 implicit none
107
108 integer :: fnom, fclos, nulnam, ierr, dateStampFromObs
109 type(struct_obs), target :: obsSpaceData
110 type(struct_oti), pointer :: oti => null()
111 real(kind=8) :: hx_dummy(1,1)
112 integer :: ncmahdr, ncmahx, ncmabdy, ncmadim, nobsout, nbrpform
113 logical :: qcvar, numHeader, numBody
114 character(len=7) :: resumeType
115
116 ! number of pressure ranges used for the thinning of aircraft (and other) data:
117 integer, parameter :: npres_ai = 5
118 integer, parameter :: npres_sw = 2
119 integer, parameter :: nai_target = 10
120 integer, parameter :: nsc_target = 10
121 integer, parameter :: nsw_target = 6
122 integer, parameter :: nto_target = 6
123 real(8) :: nai_pmax(npres_ai) = (/ 25000.0, 40000.0, 60000.0, 80000.0, 110000.0/)
124 real(8) :: nsw_pmax(npres_sw) = (/ 60000.0, 110000.0/)
125 ! For a scalar array, no layer selection will be done
126 real(8) :: nsc_pmax(1) = (/ 0.0 /)
127 real(8) :: nto_pmax(1) = (/ 0.0 /)
128
129 ! Namelist variables:
130 character(len=256) :: cmahdr ! should not be used anymore
131 character(len=256) :: cmabdy ! should not be used anymore
132 character(len=256) :: cmadim ! should not be used anymore
133 character(len=256) :: obsout ! file name for ascii output
134 character(len=256) :: brpform ! should not be used anymore
135 logical :: suprep ! choose to execute 'suprep' obs filtering
136 logical :: rejectOutsideTimeWindow ! choose to reject obs outside time window
137 logical :: thinning ! choose to apply 'extra' thinning of some obs types
138 logical :: applySatUtil ! choose to reject satellite obs based on 'util' column of stats_tovs
139 logical :: modifyAmsubObsError ! choose to modify the obs error stddev for AMSUB/MHS in the tropics
140 logical :: rejectHighLatIR ! choose to reject IR data in high latitudes
141 logical :: obsClean ! choose to remove rejected observations from files
142 logical :: writeObsFiles ! choose to update the (burp or sqlite) observation files
143 logical :: writeAsciiCmaFiles ! choose to write ascii output
144
145 NAMELIST /NAMPREPCMA/ cmahdr, cmabdy, cmadim, obsout, brpform, &
146 suprep, rejectOutsideTimeWindow, thinning, &
147 applySatUtil, modifyAmsubObsError, rejectHighLatIR, &
148 obsClean, writeObsFiles, writeAsciiCmaFiles
149
150 call ver_printNameAndVersion('prepcma','Prepare observations for LETKF')
151
152 !- 1.0 mpi
153 call mmpi_initialize
154
155 !- 1.1 timings
156 call tmg_init(mmpi_myid, 'TMG_INFO')
157 call utl_tmg_start(0,'Main')
158
159 if ( mmpi_myid == 0 ) call utl_writeStatus('PREPCMA_BEG')
160
161 !- Specify default values for namelist variables
162 cmahdr = 'NOT_DEFINED'
163 cmabdy = 'NOT_DEFINED'
164 cmadim = 'NOT_DEFINED'
165 obsout = 'NOT_DEFINED'
166 brpform = 'brpform'
167 suprep = .true.
168 rejectOutsideTimeWindow = .true.
169 thinning = .true.
170 applySatUtil = .true.
171 modifyAmsubObsError = .true.
172 rejectHighLatIR = .true.
173 obsClean = .true.
174 writeObsFiles = .false.
175 writeAsciiCmaFiles = .false.
176
177 nulnam = 0
178 ierr = fnom(nulnam,'./flnml','FTN+SEQ+R/O',0)
179 read(nulnam,nml=namprepcma,iostat=ierr)
180 if (ierr /= 0) call utl_abort('midas-prepcma: Error reading namelist')
181 if (mmpi_myid == 0) write(*,nml=namprepcma)
182 ierr = fclos(nulnam)
183
184 !- RAM disk usage
185 call ram_setup
186
187 call utl_tmg_start(10,'--Observations')
188
189 !- Set up list of elements to be assimilated and flags for rejection (from namelist)
190 call filt_setup('prepcma')
191
192 !- Observation file names and get datestamp
193 call obsf_setup(dateStampFromObs, 'prepcma')
194
195 !- Allocate obsSpaceData
196 call obs_class_initialize('ENKFMIDAS')
197 call obs_initialize( obsSpaceData, mpi_local_opt=obsf_filesSplit() )
198
199 !- Read observations
200 call utl_tmg_start(11,'----ReadObsFiles')
201 call obsf_readFiles( obsSpaceData )
202 call utl_tmg_stop(11)
203
204 numHeader = obs_numheader(obsSpaceData)
205 numBody = obs_numbody(obsSpaceData)
206 write(*,*) 'midas-prepcma: obs_numheader =', numheader
207 write(*,*) 'midas-prepcma: obs_numbody =', numbody
208
209 !- Determine if qcvar flag is expected to be present
210 resumeType = brpr_getTypeResume()
211 write(*,*) 'midas_prepcma: RESUME type =', resumeType
212 qcvar = (resumeType == 'POSTALT')
213 if (qcvar) then
214 write(*,*) 'midas_prepcma: The input file is a postalt file'
215 else
216 write(*,*) 'midas_prepcma: The input file is NOT a postalt file'
217 end if
218
219 !- Initialize TOVS processing
220 if (obs_famExist(obsSpaceData,'TO')) call tvs_setup
221
222 !- Select the elements to assimilate and apply rejection flags
223 if (suprep) call filt_suprep(obsSpaceData)
224
225 !- Allocation for TOVS
226 if (obs_famExist(obsSpaceData,'TO')) call tvs_setupAlloc(obsSpaceData)
227
228 !- Initialize obs error covariances and set flag using 'util' column of stats_tovs
229 call oer_setObsErrors(obsSpaceData, 'analysis', useTovsUtil_opt=applySatUtil) ! IN
230
231 !- Call suprep again to 'black list' channels according to 'util' column of stats_tovs
232 if (applySatUtil) call filt_suprep(obsSpaceData)
233
234 call utl_tmg_stop(10)
235
236 !- Setup timeCoord module and set dateStamp from env variable
237 call tim_setup()
238 if (tim_getDateStamp() == 0) then
239 if (dateStampFromObs > 0) then
240 ! use dateStamp from obs if not already set by env variable
241 call tim_setDateStamp(dateStampFromObs)
242 else
243 call utl_abort('midas-prepcma: DateStamp was not set')
244 end if
245 end if
246
247 !- Reject any observation outside the data assimilation window
248 if (rejectOutsideTimeWindow) then
249 call oti_setup( oti, obsSpaceData, numStep=1, &
250 headerIndexBeg=1, headerIndexEnd=obs_numheader(obsSpaceData), &
251 flagObsOutside_opt=.true. )
252 end if
253
254 !- Reject all IR radiance observation in arctic and antarctic (.i.e |lat|>60. )
255 if (rejectHighLatIR) call enkf_rejectHighLatIR(obsSpaceData)
256
257 !- Modify the obs error stddev for AMSUB in the tropics
258 if (modifyAmsubObsError) call enkf_modifyAmsubObsError(obsSpaceData)
259
260 !- Perform thinning for several observation types
261 if (thinning) then
262 ! perform thinning for aircraft observations
263 call thinning_fam(obsSpaceData, nai_pmax, nai_target, 'AI')
264 ! perform thinning for scatterometer observations
265 call thinning_fam(obsSpaceData, nsc_pmax, nsc_target, 'SC')
266 ! perform thinning for radiance observations
267 call thinning_fam(obsSpaceData, nto_pmax, nto_target, 'TO')
268 ! perform thinning for satwind observations
269 call thinning_fam(obsSpaceData, nsw_pmax, nsw_target, 'SW')
270 end if
271
272 !- Write the results
273 write(*,*)
274 write(*,*) '> midas-prepcma: writing to files'
275
276 !- Write to burp/sqlite files if requested
277 if (writeObsFiles) then
278 call obsf_writeFiles(obsSpaceData)
279 if ( obsClean ) call obsf_cleanObsFiles()
280 end if
281
282 if (writeAsciiCmaFiles) then
283
284 !- Remove all observations from obsSpaceData that will not be assimilated
285 !- But, unlike the EnKF program, do not check value of OBS_ZHA
286 if (obsClean) then
287 call obs_clean(obsSpaceData, hx_dummy, 0, -1, qcvar, checkZha_opt=.false.)
288 end if
289
290 if (mmpi_nprocs > 1) then
291 call obs_expandToMpiGlobal(obsSpaceData)
292 end if
293
294 if (mmpi_myid == 0) then
295 !- Open file for ascii output
296 nobsout = 0
297 ierr = fnom(nobsout, obsout, 'FMT+SEQ+R/W', 0)
298 call obs_print(obsSpaceData,nobsout)
299 close(nobsout)
300
301 !- Write the results in CMA format
302 ncmahdr = 0
303 ierr = fnom(ncmahdr, cmahdr, 'FTN+SEQ+UNF+R/W', 0)
304 ncmabdy = 0
305 ierr = fnom(ncmabdy, cmabdy, 'FTN+SEQ+UNF+R/W', 0)
306 ncmadim = 0
307 ierr = fnom(ncmadim, cmadim, 'FTN+SEQ+R/W', 0)
308 ncmahx = -1
309 call obs_write(obsSpaceData, hx_dummy, 0, ncmahdr, ncmabdy, ncmahx, ncmadim)
310 close(ncmahdr)
311 close(ncmabdy)
312 close(ncmadim)
313
314 !! This used to contain a .true. or .false. value indicating if observations passed the QCVar
315 !! Since, this is not the case, we can write .false.
316 nbrpform = 0
317 ierr = fnom(nbrpform, brpform, 'FTN+SEQ+R/W', 0)
318 write(nbrpform,*) .false.
319 close(nbrpform)
320
321 end if
322
323 end if
324
325 !
326 !- 3. Ending
327 !
328 write(*,*)
329 write(*,*) '> midas-prepcma: Ending'
330 call obs_finalize(obsSpaceData) ! deallocate obsSpaceData
331
332 call utl_tmg_stop(0)
333 call tmg_terminate(mmpi_myid, 'TMG_INFO')
334
335 call rpn_comm_finalize(ierr)
336
337 if ( mmpi_myid == 0 ) then
338 call utl_writeStatus('PREPCMA_END')
339 end if
340
341contains
342
343 subroutine thinning_fam(obsSpaceData, n_pmax, n_target, cfam)
344 !
345 ! :Purpose: thin the observations of the selected family
346 !
347 implicit none
348
349 ! Arguments:
350 type (struct_obs), intent(inout) :: obsSpaceData ! the data in observation space
351 real(8), intent(in) :: n_pmax(:) ! pressure levels that separate vertical layers for the thinning
352 integer, intent(in) :: n_target ! maximum desired amount of data per 3-D box
353 character(len=2), intent(in) :: cfam ! family type
354
355 ! Locals:
356 type(struct_reg) :: lsc
357 integer :: idist, n_count_thin, iai, iseed(4)
358 integer :: nobs_count, nobs_count_mpiGlobal, nobs_count_thin, nobs_count_thin_mpiGlobal
359 integer :: nrep_count, nrep_count_mpiGlobal, nrep_count_thin, nrep_count_thin_mpiGlobal
360 integer :: headerIndex, bodyIndex, bodyIndexBeg, bodyIndexEnd
361 integer :: iblock, codeType, ilat, incr, ipres, nblocksum, npres, nsize, num_stn
362 logical :: count_obs, allRejected
363 real(4) :: lat_r4, lon_r4
364 real(8) :: pressure, rannum
365 real(8), allocatable :: latcenter(:), latmin(:), latmax(:), ranvals(:)
366 integer, allocatable :: nblockoffset(:), nlonblock(:)
367 integer, allocatable :: ai_indices(:,:), nstation(:,:), nstationMpiGlobal(:,:)
368 real(8), allocatable :: keep_ai(:,:)
369
370 ! box size that is used for observation thinning
371 ! (the numerator is an approximate distance in km)
372 real(8), parameter :: r0_count_km = 200.0/(2**0.5)
373 ! next two parameters are not used in this program
374 real(8), parameter :: r1_dum = 1.0
375 real(8), parameter :: rz_dum = 1.0
376
377 num_stn = obs_numheader(obsSpaceData)
378 npres = size(n_pmax,1)
379 write(*,*) 'Start thinning for ', cfam, ' data'
380 ! at this stage we still have many radiance channels
381 ! that will be rejected at a later stage.
382 if (cfam == 'XX') then
383 ! never getting here
384 write(*,*) 'count individual observations'
385 count_obs = .true.
386 else
387 write(*,*) 'count the number of reports'
388 count_obs = .false.
389 end if
390
391 call reg_init_struct(lsc, r0_count_km, r1_dum, rz_dum)
392 if (mmpi_myid == 0) write(*,*) 'number of latitude bands: ', lsc%nlatband
393 nsize = lsc%nlatband
394 allocate(latmin(nsize))
395 allocate(latmax(nsize))
396 allocate(latcenter(nsize))
397 allocate(nlonblock(nsize))
398 allocate(nblockoffset(nsize))
399
400 call reg_getlatitude(lsc%r0_rad, lsc%nlatband, latmin, latcenter, latmax)
401 if (mmpi_myid == 0) write(*,*) 'number of latitude bands: ',lsc%nlatband
402 do ilat = 1, lsc%nlatband
403 if (mmpi_myid == 0) write(*,*) ' band: ', ilat, ' latitude between ', latmin(ilat), latmax(ilat)
404 end do
405 call reg_getblock(lsc%nlatband, lsc%r0_rad, latmin, latmax, nlonblock)
406 nblocksum = 0
407 do ilat = 1, lsc%nlatband
408 nblockoffset(ilat) = nblocksum
409 nblocksum = nblocksum + nlonblock(ilat)
410 if (mmpi_myid == 0) write(*,*) 'latband: ', ilat, ' no of blocks: ', nlonblock(ilat)
411 end do
412 if (mmpi_myid == 0) write(*,*) 'total number of blocks: ', nblocksum
413
414 nrep_count = 0
415 nobs_count = 0
416
417 allocate(nstation(nblocksum, npres))
418 allocate(keep_ai(nblocksum, npres))
419
420 nstation=0
421 ! keep_ai = 1 corresponds to no thinning
422 keep_ai(:,:) = 1.0
423
424 allocate(ai_indices(num_stn,3))
425
426 header_loop: do headerIndex = 1, num_stn
427 codeType= obs_headElem_i(obsSpaceData, obs_ity, headerIndex)
428 if ( (cfam=='AI' .and. (codeType==42 .or. codeType==128 .or. &
429 codeType==157 .or. codeType==177)) .or. &
430 (cfam=='SC' .and. codeType==254) .or. &
431 (cfam=='SW' .and. (codeType==88 .or. codeType==188)) .or. &
432 (cfam=='TO' .and. tvs_isIdBurpTovs(codeType)) ) then
433
434 ! skip this header if all observations already rejected
435 bodyIndexBeg = obs_headElem_i(obsSpaceData, obs_rln, headerIndex)
436 bodyIndexEnd = obs_headElem_i(obsSpaceData, obs_nlv, headerIndex) + bodyIndexBeg - 1
437 allRejected = .true.
438 body_loop: do bodyIndex = bodyIndexBeg, bodyIndexEnd
439 if (obs_bodyElem_i(obsSpaceData, obs_ass, bodyIndex) == obs_assimilated) then
440 allRejected = .false.
441 exit body_loop
442 end if
443 end do body_loop
444 if (allRejected) cycle header_loop
445
446 lat_r4 = obs_headElem_r(obsSpaceData, obs_lat, headerIndex)
447 lon_r4 = obs_headElem_r(obsSpaceData, obs_lon, headerIndex)
448 call reg_locatestn(lsc%r0_rad, lat_r4, lon_r4, &
449 lsc%nlatband, nlonblock, &
450 nblockoffset, iblock)
451 ! note that all data from the aircraft are at the same pressure
452 if (npres == 1) then
453 ipres = 1
454 else
455 pressure= obs_bodyElem_r(obsSpaceData, obs_ppp, bodyIndexBeg)
456 ipres = 1
457 do
458 if ((n_pmax(ipres) > pressure) .or. (ipres > npres)) exit
459 ipres = ipres + 1
460 end do
461 end if
462 if (ipres <= npres) then
463 if (count_obs) then
464 incr = obs_headElem_i(obsSpaceData, obs_nlv, headerIndex)
465 else
466 incr = 1
467 end if
468 nstation(iblock,ipres) = nstation(iblock,ipres) + incr
469 nobs_count = nobs_count + obs_headElem_i(obsSpaceData, obs_nlv, headerIndex)
470 nrep_count = nrep_count + 1
471 ai_indices(nrep_count, 1) = headerIndex
472 ai_indices(nrep_count, 2) = iblock
473 ai_indices(nrep_count, 3) = ipres
474 end if
475 end if
476 end do header_loop
477
478 ! do mpi communication of the accumulators
479 nsize = nblocksum * npres
480 allocate(nstationMpiGlobal(nblocksum, npres))
481 call rpn_comm_allreduce(nstation, nstationMpiGlobal, nsize, &
482 'mpi_integer','mpi_sum', 'GRID', ierr)
483 call rpn_comm_allreduce(nrep_count, nrep_count_mpiGlobal, 1, &
484 'mpi_integer','mpi_sum', 'GRID', ierr)
485 call rpn_comm_allreduce(nobs_count, nobs_count_mpiGlobal, 1, &
486 'mpi_integer','mpi_sum', 'GRID', ierr)
487
488 write(*,*) 'total number of ', cfam, ' reports (local and mpiglobal): ', &
489 nrep_count, nrep_count_mpiGlobal
490 allocate(ranvals(nrep_count))
491 write(*,*) 'total number of ', cfam, ' observations (local and mpiglobal): ', &
492 nobs_count, nobs_count_mpiGlobal
493
494 n_count_thin = 0
495 do iblock = 1, nblocksum
496 do ipres = 1, npres
497 if (nstationMpiGlobal(iblock,ipres) .ge. 1) then
498 if (mmpi_myid == 0) write(*,*) 'block ipres and count: ',iblock,ipres, &
499 nstationMpiGlobal(iblock,ipres)
500 if (nstationMpiGlobal(iblock,ipres) > n_target) then
501 keep_ai(iblock,ipres) = dble(n_target) / dble(nstationMpiGlobal(iblock,ipres))
502 n_count_thin = n_count_thin + n_target
503 else
504 n_count_thin = n_count_thin + nstationMpiGlobal(iblock,ipres)
505 end if
506 end if
507 end do
508 end do
509
510 if (count_obs) then
511 write(*,*) 'Estimated remaining number of ', cfam, ' observations (mpiGlobal): ', n_count_thin
512 else
513 write(*,*) 'Estimated remaining number of ', cfam, ' reports (mpiGlobal): ', n_count_thin
514 end if
515
516 nrep_count_thin = 0
517 nobs_count_thin = 0
518
519 idist = 1
520 iseed(1) = 1
521 iseed(2) = 5
522 iseed(3) = 9
523 iseed(4) = 11
524 call dlarnv(idist,iseed,nrep_count,ranvals)
525 do iai = 1, nrep_count
526 headerIndex = ai_indices(iai,1)
527 iblock = ai_indices(iai,2)
528 ipres = ai_indices(iai,3)
529 rannum = ranvals(iai)
530 if (rannum <= keep_ai(iblock,ipres)) then
531 nrep_count_thin = nrep_count_thin + 1
532 nobs_count_thin = nobs_count_thin + obs_headElem_i(obsSpaceData, obs_nlv, headerIndex)
533 else
534 ! reject the profile
535 bodyIndexBeg = obs_headElem_i(obsSpaceData, obs_rln, headerIndex)
536 bodyIndexEnd = obs_headElem_i(obsSpaceData, obs_nlv, headerIndex) + bodyIndexBeg - 1
537 do bodyIndex = bodyIndexBeg, bodyIndexEnd
538 call obs_bodySet_i(obsSpaceData, obs_ass, bodyIndex, obs_notAssimilated)
539 ! also set the 'rejected by selection process' flag (bit 11)
540 call obs_bodySet_i( obsSpaceData, obs_flg, bodyIndex, &
541 ibset( obs_bodyElem_i( obsSpaceData, obs_flg, bodyIndex ), 11) )
542 end do
543 end if
544 end do
545
546 ! mpi communication of accumulators
547 call rpn_comm_allreduce(nrep_count_thin, nrep_count_thin_mpiGlobal, 1, &
548 'mpi_integer','mpi_sum', 'GRID', ierr)
549 call rpn_comm_allreduce(nobs_count_thin, nobs_count_thin_mpiGlobal, 1, &
550 'mpi_integer','mpi_sum', 'GRID', ierr)
551
552 write(*,*) 'True remaining number of ', cfam, ' reports (local, mpiGlobal): ', &
553 nrep_count_thin, nrep_count_thin_mpiGlobal
554 write(*,*) 'True remaining number of ', cfam, ' observations (local, mpiGlobal): ', &
555 nobs_count_thin, nobs_count_thin_mpiGlobal
556
557 deallocate(ranvals)
558 deallocate(latmin)
559 deallocate(latmax)
560 deallocate(latcenter)
561 deallocate(nlonblock)
562 deallocate(nblockoffset)
563 deallocate(nstation)
564 deallocate(nstationMpiGlobal)
565 deallocate(keep_ai)
566 deallocate(ai_indices)
567
568 end subroutine thinning_fam
569
570end program
