-
Notifications
You must be signed in to change notification settings - Fork 150
/
Copy pathassim_tools_mod.f90
2849 lines (2327 loc) · 114 KB
/
assim_tools_mod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! DART software - Copyright UCAR. This open source software is provided
! by UCAR, "as is", without charge, subject to all terms of use at
! http://www.image.ucar.edu/DAReS/DART/DART_download
!> A variety of operations required by assimilation.
module assim_tools_mod
!> \defgroup assim_tools assim_tools_mod
!>
!> @{
use types_mod, only : r8, i8, digits12, PI, missing_r8
use options_mod, only : get_missing_ok_status
use utilities_mod, only : file_exist, get_unit, check_namelist_read, do_output, &
find_namelist_in_file, error_handler, &
E_ERR, E_MSG, nmlfileunit, do_nml_file, do_nml_term, &
open_file, close_file, timestamp
use sort_mod, only : index_sort
use random_seq_mod, only : random_seq_type, random_gaussian, init_random_seq, &
random_uniform
use obs_sequence_mod, only : obs_sequence_type, obs_type, get_num_copies, get_num_qc, &
init_obs, get_obs_from_key, get_obs_def, get_obs_values, &
destroy_obs
use obs_def_mod, only : obs_def_type, get_obs_def_location, get_obs_def_time, &
get_obs_def_error_variance, get_obs_def_type_of_obs
use obs_kind_mod, only : get_num_types_of_obs, get_index_for_type_of_obs, &
get_quantity_for_type_of_obs, assimilate_this_type_of_obs
use cov_cutoff_mod, only : comp_cov_factor
use reg_factor_mod, only : comp_reg_factor
use obs_impact_mod, only : allocate_impact_table, read_impact_table, free_impact_table
use sampling_error_correction_mod, only : get_sampling_error_table_size, &
read_sampling_error_correction
use location_mod, only : location_type, get_close_type, query_location, &
operator(==), set_location_missing, write_location, &
LocationDims, is_vertical, vertical_localization_on, &
set_vertical, has_vertical_choice, get_close_init, &
get_vertical_localization_coord, get_close_destroy, &
set_vertical_localization_coord
use ensemble_manager_mod, only : ensemble_type, get_my_num_vars, get_my_vars, &
compute_copy_mean_var, get_var_owner_index, &
map_pe_to_task
use mpi_utilities_mod, only : my_task_id, broadcast_send, broadcast_recv, &
sum_across_tasks, task_count, start_mpi_timer, &
read_mpi_timer, task_sync
use adaptive_inflate_mod, only : do_obs_inflate, do_single_ss_inflate, do_ss_inflate, &
do_varying_ss_inflate, &
update_inflation, update_single_state_space_inflation, &
update_varying_state_space_inflation, &
inflate_ens, adaptive_inflate_type, &
deterministic_inflate, solve_quadratic
use time_manager_mod, only : time_type, get_time
use assim_model_mod, only : get_state_meta_data, &
get_close_obs, get_close_state, &
convert_vertical_obs, convert_vertical_state
use distributed_state_mod, only : create_mean_window, free_mean_window
use quality_control_mod, only : good_dart_qc, DARTQC_FAILED_VERT_CONVERT
implicit none
private
public :: filter_assim, &
set_assim_tools_trace, &
test_state_copies, &
update_ens_from_weights
! Indicates if module initialization subroutine has been called yet
logical :: module_initialized = .false.
integer :: print_timestamps = 0
integer :: print_trace_details = 0
! True if random sequence needs to be initialized
logical :: first_inc_ran_call = .true.
type (random_seq_type) :: inc_ran_seq
integer :: num_types = 0
real(r8), allocatable :: cutoff_list(:)
logical :: has_special_cutoffs
logical :: close_obs_caching = .true.
! true if we have multiple vert choices and we're doing vertical localization
! (make it a local variable so we don't keep making subroutine calls)
logical :: is_doing_vertical_conversion = .false.
character(len=512) :: msgstring, msgstring2, msgstring3
! Need to read in table for off-line based sampling correction and store it
integer :: sec_table_size
real(r8), allocatable :: exp_true_correl(:), alpha(:)
! if adjust_obs_impact is true, read in triplets from the ascii file
! and fill this 2d impact table.
real(r8), allocatable :: obs_impact_table(:,:)
character(len=*), parameter :: source = 'assim_tools_mod.f90'
!============================================================================
!---- namelist with default values
! Filter kind selects type of observation space filter
! 1 = EAKF filter
! 2 = ENKF
! 3 = Kernel filter
! 4 = particle filter
! 5 = random draw from posterior
! 6 = deterministic draw from posterior with fixed kurtosis
! 8 = Rank Histogram Filter (see Anderson 2011)
!
! special_localization_obs_types -> Special treatment for the specified observation types
! special_localization_cutoffs -> Different cutoff value for each specified obs type
!
integer :: filter_kind = 1
real(r8) :: cutoff = 0.2_r8
logical :: sort_obs_inc = .false.
logical :: spread_restoration = .false.
logical :: sampling_error_correction = .false.
integer :: adaptive_localization_threshold = -1
real(r8) :: adaptive_cutoff_floor = 0.0_r8
integer :: print_every_nth_obs = 0
! since this is in the namelist, it has to have a fixed size.
integer, parameter :: MAX_ITEMS = 300
character(len = 129) :: special_localization_obs_types(MAX_ITEMS)
real(r8) :: special_localization_cutoffs(MAX_ITEMS)
logical :: output_localization_diagnostics = .false.
character(len = 129) :: localization_diagnostics_file = "localization_diagnostics"
! Following only relevant for filter_kind = 8
logical :: rectangular_quadrature = .true.
logical :: gaussian_likelihood_tails = .false.
! False by default; if true, expect to read in an ascii table
! to adjust the impact of obs on other state vector and obs values.
logical :: adjust_obs_impact = .false.
character(len=256) :: obs_impact_filename = ''
logical :: allow_any_impact_values = .false.
! These next two only affect models with multiple options
! for vertical localization:
!
! "convert_state" is false by default; it depends on the model
! what is faster - do the entire state up front and possibly
! do unneeded work, or do the conversion during the assimilation
! loop. we think this depends heavily on how much of the state
! is going to be adjusted by the obs. for a global model
! we think false may be better; for a regional model with
! a lot of obs and full coverage true may be better.
!
! "convert_obs" is true by default; in general it seems to
! be better for each task to convert the obs vertical before
! going into the loop but again this depends on how many
! obs per task and whether the mean is distributed or
! replicated on each task.
logical :: convert_all_state_verticals_first = .false.
logical :: convert_all_obs_verticals_first = .true.
! Not in the namelist; this var disables the experimental
! linear and spherical case code in the adaptive localization
! sections. to try out the alternatives, set this to .false.
logical :: only_area_adapt = .true.
! Option to distribute the mean. If 'false' each task will have a full
! copy of the ensemble mean, which speeds models doing vertical conversion.
! If 'true' the mean will be spread across all tasks which reduces the
! memory needed per task but requires communication if the mean is used
! for vertical conversion. We have changed the default to be .false.
! compared to previous versions of this namelist item.
logical :: distribute_mean = .false.
namelist / assim_tools_nml / filter_kind, cutoff, sort_obs_inc, &
spread_restoration, sampling_error_correction, &
adaptive_localization_threshold, adaptive_cutoff_floor, &
print_every_nth_obs, rectangular_quadrature, gaussian_likelihood_tails, &
output_localization_diagnostics, localization_diagnostics_file, &
special_localization_obs_types, special_localization_cutoffs, &
distribute_mean, close_obs_caching, &
adjust_obs_impact, obs_impact_filename, allow_any_impact_values, &
convert_all_state_verticals_first, convert_all_obs_verticals_first
!============================================================================
contains
!-------------------------------------------------------------
subroutine assim_tools_init()
integer :: iunit, io, i, j
integer :: num_special_cutoff, type_index
logical :: cache_override = .false.
! do this up front
module_initialized = .true.
! give these guys initial values at runtime *before* we read
! in the namelist. this is to help detect how many items are
! actually given in the namelist.
special_localization_obs_types(:) = 'null'
special_localization_cutoffs(:) = missing_r8
! Read the namelist entry
call find_namelist_in_file("input.nml", "assim_tools_nml", iunit)
read(iunit, nml = assim_tools_nml, iostat = io)
call check_namelist_read(iunit, io, "assim_tools_nml")
! Write the namelist values to the log file
if (do_nml_file()) write(nmlfileunit, nml=assim_tools_nml)
if (do_nml_term()) write( * , nml=assim_tools_nml)
! Forcing distributed_mean for single processor.
! Note null_win_mod.f90 ignores distibute_mean.
if (task_count() == 1) distribute_mean = .true.
! FOR NOW, can only do spread restoration with filter option 1 (need to extend this)
if(spread_restoration .and. .not. filter_kind == 1) then
write(msgstring, *) 'cannot combine spread_restoration and filter_kind ', filter_kind
call error_handler(E_ERR,'assim_tools_init:', msgstring, source)
endif
! allocate a list in all cases - even the ones where there is only
! a single cutoff value. note that in spite of the name these
! are specific types (e.g. RADIOSONDE_TEMPERATURE, AIRCRAFT_TEMPERATURE)
! because that's what get_close() is passed. and because i've confused
! myself several times -- we define generic kinds starting at 0, but
! the specific types are autogenerated and always start at 1. so the
! cutoff list is never (0:num_types); it is always (num_types).
num_types = get_num_types_of_obs()
allocate(cutoff_list(num_types))
cutoff_list(:) = cutoff
has_special_cutoffs = .false.
! Go through special-treatment observation kinds, if any.
num_special_cutoff = 0
j = 0
do i = 1, MAX_ITEMS
if(special_localization_obs_types(i) == 'null') exit
if(special_localization_cutoffs(i) == MISSING_R8) then
write(msgstring, *) 'cutoff value', i, ' is uninitialized.'
call error_handler(E_ERR,'assim_tools_init:', &
'special cutoff namelist for types and distances do not match', &
source, &
text2='kind = '//trim(special_localization_obs_types(i)), &
text3=trim(msgstring))
endif
j = j + 1
enddo
num_special_cutoff = j
if (num_special_cutoff > 0) has_special_cutoffs = .true.
do i = 1, num_special_cutoff
type_index = get_index_for_type_of_obs(special_localization_obs_types(i))
if (type_index < 0) then
write(msgstring, *) 'unrecognized TYPE_ in the special localization namelist:'
call error_handler(E_ERR,'assim_tools_init:', msgstring, source, &
text2=trim(special_localization_obs_types(i)))
endif
cutoff_list(type_index) = special_localization_cutoffs(i)
end do
! cannot cache previous obs location if different obs types have different
! localization radii. change it to false, and warn user why.
if (has_special_cutoffs .and. close_obs_caching) then
cache_override = .true.
close_obs_caching = .false.
endif
if(sampling_error_correction) then
sec_table_size = get_sampling_error_table_size()
allocate(exp_true_correl(sec_table_size), alpha(sec_table_size))
! we can't read the table here because we don't have access to the ens_size
endif
is_doing_vertical_conversion = (has_vertical_choice() .and. vertical_localization_on())
call log_namelist_selections(num_special_cutoff, cache_override)
end subroutine assim_tools_init
!-------------------------------------------------------------
subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys, &
ens_size, num_groups, obs_val_index, inflate, ENS_MEAN_COPY, ENS_SD_COPY, &
ENS_INF_COPY, ENS_INF_SD_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END, OBS_PRIOR_VAR_START, &
OBS_PRIOR_VAR_END, inflate_only)
type(ensemble_type), intent(inout) :: ens_handle, obs_ens_handle
type(obs_sequence_type), intent(in) :: obs_seq
integer, intent(in) :: keys(:)
integer, intent(in) :: ens_size, num_groups, obs_val_index
! JLA: At present, this only needs to be inout because of the possible use of
! non-determinstic obs_space adaptive inflation that is not currently supported.
! Implementing that would require communication of the info about the inflation
! values as each observation updated them.
type(adaptive_inflate_type), intent(inout) :: inflate
integer, intent(in) :: ENS_MEAN_COPY, ENS_SD_COPY, ENS_INF_COPY
integer, intent(in) :: ENS_INF_SD_COPY
integer, intent(in) :: OBS_KEY_COPY, OBS_GLOBAL_QC_COPY
integer, intent(in) :: OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END
integer, intent(in) :: OBS_PRIOR_VAR_START, OBS_PRIOR_VAR_END
logical, intent(in) :: inflate_only
! changed the ensemble sized things here to allocatable
real(r8) :: obs_prior(ens_size), obs_inc(ens_size), updated_ens(ens_size)
real(r8) :: final_factor
real(r8) :: net_a(num_groups), correl(num_groups)
real(r8) :: obs(1), obs_err_var, my_inflate, my_inflate_sd
real(r8) :: obs_qc, cutoff_rev, cutoff_orig
real(r8) :: orig_obs_prior_mean(num_groups), orig_obs_prior_var(num_groups)
real(r8) :: obs_prior_mean(num_groups), obs_prior_var(num_groups)
real(r8) :: vertvalue_obs_in_localization_coord, whichvert_real
real(r8), allocatable :: close_obs_dist(:)
real(r8), allocatable :: close_state_dist(:)
integer(i8) :: state_index
integer(i8), allocatable :: my_state_indx(:)
integer(i8), allocatable :: my_obs_indx(:)
integer :: my_num_obs, i, j, owner, owners_index, my_num_state
integer :: obs_mean_index, obs_var_index
integer :: grp_beg(num_groups), grp_end(num_groups), grp_size, grp_bot, grp_top, group
integer :: num_close_obs, obs_index, num_close_states
integer :: last_num_close_obs, last_num_close_states
integer :: base_obs_kind, base_obs_type, nth_obs
integer :: num_close_obs_cached, num_close_states_cached
integer :: num_close_obs_calls_made, num_close_states_calls_made
integer :: whichvert_obs_in_localization_coord
integer :: istatus, localization_unit
integer, allocatable :: close_obs_ind(:)
integer, allocatable :: close_state_ind(:)
integer, allocatable :: my_obs_kind(:)
integer, allocatable :: my_obs_type(:)
integer, allocatable :: my_state_kind(:)
integer, allocatable :: vstatus(:)
type(location_type) :: base_obs_loc, last_base_obs_loc, last_base_states_loc
type(location_type) :: dummyloc
type(location_type), allocatable :: my_obs_loc(:)
type(location_type), allocatable :: my_state_loc(:)
type(get_close_type) :: gc_obs, gc_state
type(obs_type) :: observation
type(obs_def_type) :: obs_def
type(time_type) :: obs_time
logical :: allow_missing_in_state
logical :: local_single_ss_inflate
logical :: local_varying_ss_inflate
logical :: local_ss_inflate
logical :: local_obs_inflate
! allocate rather than dump all this on the stack
allocate(close_obs_dist( obs_ens_handle%my_num_vars), &
close_obs_ind( obs_ens_handle%my_num_vars), &
vstatus( obs_ens_handle%my_num_vars), &
my_obs_indx( obs_ens_handle%my_num_vars), &
my_obs_kind( obs_ens_handle%my_num_vars), &
my_obs_type( obs_ens_handle%my_num_vars), &
my_obs_loc( obs_ens_handle%my_num_vars))
allocate(close_state_dist( ens_handle%my_num_vars), &
close_state_ind( ens_handle%my_num_vars), &
my_state_indx( ens_handle%my_num_vars), &
my_state_kind( ens_handle%my_num_vars), &
my_state_loc( ens_handle%my_num_vars))
! end alloc
! Initialize assim_tools_module if needed
if (.not. module_initialized) call assim_tools_init()
!HK make window for mpi one-sided communication
! used for vertical conversion in get_close_obs
! Need to give create_mean_window the mean copy
call create_mean_window(ens_handle, ENS_MEAN_COPY, distribute_mean)
! filter kinds 1 and 8 return sorted increments, however non-deterministic
! inflation can scramble these. the sort is expensive, so help users get better
! performance by rejecting namelist combinations that do unneeded work.
if (sort_obs_inc) then
if(deterministic_inflate(inflate) .and. ((filter_kind == 1) .or. (filter_kind == 8))) then
write(msgstring, *) 'With a deterministic filter [assim_tools_nml:filter_kind = ',filter_kind,']'
write(msgstring2, *) 'and deterministic inflation [filter_nml:inf_deterministic = .TRUE.]'
write(msgstring3, *) 'assim_tools_nml:sort_obs_inc = .TRUE. is not needed and is expensive.'
call error_handler(E_MSG,'', '') ! whitespace
call error_handler(E_MSG,'WARNING filter_assim:', msgstring, source, &
text2=msgstring2,text3=msgstring3)
call error_handler(E_MSG,'', '') ! whitespace
sort_obs_inc = .FALSE.
endif
endif
! Open the localization diagnostics file
if(output_localization_diagnostics .and. my_task_id() == 0) &
localization_unit = open_file(localization_diagnostics_file, action = 'append')
! For performance, make local copies of these settings which
! are really in the inflate derived type.
local_single_ss_inflate = do_single_ss_inflate(inflate)
local_varying_ss_inflate = do_varying_ss_inflate(inflate)
local_ss_inflate = do_ss_inflate(inflate)
local_obs_inflate = do_obs_inflate(inflate)
! Default to printing nothing
nth_obs = -1
! Divide ensemble into num_groups groups.
! make sure the number of groups and ensemble size result in
! at least 2 members in each group (to avoid divide by 0) and
! that the groups all have the same number of members.
grp_size = ens_size / num_groups
if ((grp_size * num_groups) /= ens_size) then
write(msgstring, *) 'The number of ensemble members must divide into the number of groups evenly.'
write(msgstring2, *) 'Ensemble size = ', ens_size, ' Number of groups = ', num_groups
write(msgstring3, *) 'Change number of groups or ensemble size to avoid remainders.'
call error_handler(E_ERR,'filter_assim:', msgstring, source, &
text2=msgstring2,text3=msgstring3)
endif
if (grp_size < 2) then
write(msgstring, *) 'There must be at least 2 ensemble members in each group.'
write(msgstring2, *) 'Ensemble size = ', ens_size, ' Number of groups = ', num_groups
write(msgstring3, *) 'results in < 2 members/group. Decrease number of groups or increase ensemble size'
call error_handler(E_ERR,'filter_assim:', msgstring, source, &
text2=msgstring2,text3=msgstring3)
endif
do group = 1, num_groups
grp_beg(group) = (group - 1) * grp_size + 1
grp_end(group) = grp_beg(group) + grp_size - 1
enddo
! Put initial value of state space inflation in copy normally used for SD
! This is to avoid weird storage footprint in filter
ens_handle%copies(ENS_SD_COPY, :) = ens_handle%copies(ENS_INF_COPY, :)
! For single state or obs space inflation, the inflation is like a token
! Gets passed from the processor with a given obs on to the next
if(local_single_ss_inflate) then
my_inflate = ens_handle%copies(ENS_INF_COPY, 1)
my_inflate_sd = ens_handle%copies(ENS_INF_SD_COPY, 1)
end if
! Get info on my number and indices for obs
my_num_obs = get_my_num_vars(obs_ens_handle)
call get_my_vars(obs_ens_handle, my_obs_indx)
! Construct an observation temporary
call init_obs(observation, get_num_copies(obs_seq), get_num_qc(obs_seq))
! Get the locations for all of my observations
! HK I would like to move this to before the calculation of the forward operator so you could
! overwrite the vertical location with the required localization vertical coordinate when you
! do the forward operator calculation
call get_my_obs_loc(obs_ens_handle, obs_seq, keys, my_obs_loc, my_obs_kind, my_obs_type, obs_time)
if (convert_all_obs_verticals_first .and. is_doing_vertical_conversion) then
! convert the vertical of all my observations to the localization coordinate
if (obs_ens_handle%my_num_vars > 0) then
call convert_vertical_obs(ens_handle, obs_ens_handle%my_num_vars, my_obs_loc, &
my_obs_kind, my_obs_type, get_vertical_localization_coord(), vstatus)
do i = 1, obs_ens_handle%my_num_vars
if (good_dart_qc(nint(obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, i)))) then
!> @todo Can I just use the OBS_GLOBAL_QC_COPY? Is it ok to skip the loop?
if (vstatus(i) /= 0) obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, i) = DARTQC_FAILED_VERT_CONVERT
endif
enddo
endif
endif
! Get info on my number and indices for state
my_num_state = get_my_num_vars(ens_handle)
call get_my_vars(ens_handle, my_state_indx)
! Get the location and kind of all my state variables
do i = 1, ens_handle%my_num_vars
call get_state_meta_data(my_state_indx(i), my_state_loc(i), my_state_kind(i))
end do
!> optionally convert all state location verticals
if (convert_all_state_verticals_first .and. is_doing_vertical_conversion) then
if (ens_handle%my_num_vars > 0) then
call convert_vertical_state(ens_handle, ens_handle%my_num_vars, my_state_loc, my_state_kind, &
my_state_indx, get_vertical_localization_coord(), istatus)
endif
endif
! Get mean and variance of each group's observation priors for adaptive inflation
! Important that these be from before any observations have been used
if(local_ss_inflate) then
do group = 1, num_groups
obs_mean_index = OBS_PRIOR_MEAN_START + group - 1
obs_var_index = OBS_PRIOR_VAR_START + group - 1
call compute_copy_mean_var(obs_ens_handle, grp_beg(group), grp_end(group), &
obs_mean_index, obs_var_index)
end do
endif
! Initialize the method for getting state variables close to a given ob on my process
if (has_special_cutoffs) then
call get_close_init(gc_state, my_num_state, 2.0_r8*cutoff, my_state_loc, 2.0_r8*cutoff_list)
else
call get_close_init(gc_state, my_num_state, 2.0_r8*cutoff, my_state_loc)
endif
! Initialize the method for getting obs close to a given ob on my process
if (has_special_cutoffs) then
call get_close_init(gc_obs, my_num_obs, 2.0_r8*cutoff, my_obs_loc, 2.0_r8*cutoff_list)
else
call get_close_init(gc_obs, my_num_obs, 2.0_r8*cutoff, my_obs_loc)
endif
if (close_obs_caching) then
! Initialize last obs and state get_close lookups, to take advantage below
! of sequential observations at the same location (e.g. U,V, possibly T,Q)
! (this is getting long enough it probably should go into a subroutine. nsc.)
last_base_obs_loc = set_location_missing()
last_base_states_loc = set_location_missing()
last_num_close_obs = -1
last_num_close_states = -1
num_close_obs_cached = 0
num_close_states_cached = 0
num_close_obs_calls_made = 0
num_close_states_calls_made = 0
endif
allow_missing_in_state = get_missing_ok_status()
! Loop through all the (global) observations sequentially
SEQUENTIAL_OBS: do i = 1, obs_ens_handle%num_vars
! Some compilers do not like mod by 0, so test first.
if (print_every_nth_obs > 0) nth_obs = mod(i, print_every_nth_obs)
! If requested, print out a message every Nth observation
! to indicate progress is being made and to allow estimates
! of how long the assim will take.
if (nth_obs == 0) then
write(msgstring, '(2(A,I8))') 'Processing observation ', i, &
' of ', obs_ens_handle%num_vars
if (print_timestamps == 0) then
call error_handler(E_MSG,'filter_assim',msgstring)
else
call timestamp(trim(msgstring), pos="brief")
endif
endif
! Every pe has information about the global obs sequence
call get_obs_from_key(obs_seq, keys(i), observation)
call get_obs_def(observation, obs_def)
base_obs_loc = get_obs_def_location(obs_def)
obs_err_var = get_obs_def_error_variance(obs_def)
base_obs_type = get_obs_def_type_of_obs(obs_def)
if (base_obs_type > 0) then
base_obs_kind = get_quantity_for_type_of_obs(base_obs_type)
else
call get_state_meta_data(-1 * int(base_obs_type,i8), dummyloc, base_obs_kind) ! identity obs
endif
! Get the value of the observation
call get_obs_values(observation, obs, obs_val_index)
! Find out who has this observation and where it is
call get_var_owner_index(ens_handle, int(i,i8), owner, owners_index)
! Following block is done only by the owner of this observation
!-----------------------------------------------------------------------
if(ens_handle%my_pe == owner) then
! each task has its own subset of all obs. if they were converted in the
! vertical up above, then we need to broadcast the new values to all the other
! tasks so they're computing the right distances when applying the increments.
if (is_doing_vertical_conversion) then
vertvalue_obs_in_localization_coord = query_location(my_obs_loc(owners_index), "VLOC")
whichvert_obs_in_localization_coord = query_location(my_obs_loc(owners_index), "WHICH_VERT")
else
vertvalue_obs_in_localization_coord = 0.0_r8
whichvert_obs_in_localization_coord = 0
endif
obs_qc = obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, owners_index)
! Only value of 0 for DART QC field should be assimilated
IF_QC_IS_OKAY: if(nint(obs_qc) ==0) then
obs_prior = obs_ens_handle%copies(1:ens_size, owners_index)
! Note that these are before DA starts, so can be different from current obs_prior
orig_obs_prior_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START: &
OBS_PRIOR_MEAN_END, owners_index)
orig_obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START: &
OBS_PRIOR_VAR_END, owners_index)
endif IF_QC_IS_OKAY
!Broadcast the info from this obs to all other processes
! orig_obs_prior_mean and orig_obs_prior_var only used with adaptive inflation
! my_inflate and my_inflate_sd only used with single state space inflation
! vertvalue_obs_in_localization_coord and whichvert_real only used for vertical
! coordinate transformation
whichvert_real = real(whichvert_obs_in_localization_coord, r8)
call broadcast_send(map_pe_to_task(ens_handle, owner), obs_prior, &
orig_obs_prior_mean, orig_obs_prior_var, &
scalar1=obs_qc, scalar2=vertvalue_obs_in_localization_coord, &
scalar3=whichvert_real, scalar4=my_inflate, scalar5=my_inflate_sd)
! Next block is done by processes that do NOT own this observation
!-----------------------------------------------------------------------
else
call broadcast_recv(map_pe_to_task(ens_handle, owner), obs_prior, &
orig_obs_prior_mean, orig_obs_prior_var, &
scalar1=obs_qc, scalar2=vertvalue_obs_in_localization_coord, &
scalar3=whichvert_real, scalar4=my_inflate, scalar5=my_inflate_sd)
whichvert_obs_in_localization_coord = nint(whichvert_real)
endif
!-----------------------------------------------------------------------
! Everybody is doing this section, cycle if qc is bad
if(nint(obs_qc) /= 0) cycle SEQUENTIAL_OBS
!> all tasks must set the converted vertical values into the 'base' version of this loc
!> because that's what we pass into the get_close_xxx() routines below.
if (is_doing_vertical_conversion) &
call set_vertical(base_obs_loc, vertvalue_obs_in_localization_coord, whichvert_obs_in_localization_coord)
! Compute observation space increments for each group
do group = 1, num_groups
grp_bot = grp_beg(group); grp_top = grp_end(group)
call obs_increment(obs_prior(grp_bot:grp_top), grp_size, obs(1), &
obs_err_var, obs_inc(grp_bot:grp_top), inflate, my_inflate, &
my_inflate_sd, net_a(group))
! Also compute prior mean and variance of obs for efficiency here
obs_prior_mean(group) = sum(obs_prior(grp_bot:grp_top)) / grp_size
obs_prior_var(group) = sum((obs_prior(grp_bot:grp_top) - obs_prior_mean(group))**2) / &
(grp_size - 1)
if (obs_prior_var(group) < 0.0_r8) obs_prior_var(group) = 0.0_r8
end do
! Compute updated values for single state space inflation
if(local_single_ss_inflate) then
! Update for each group separately
do group = 1, num_groups
call update_single_state_space_inflation(inflate, my_inflate, my_inflate_sd, &
ens_handle%copies(ENS_SD_COPY, 1), orig_obs_prior_mean(group), &
orig_obs_prior_var(group), obs(1), obs_err_var, grp_size, inflate_only)
end do
endif
! Adaptive localization needs number of other observations within localization radius.
! Do get_close_obs first, even though state space increments are computed before obs increments.
call get_close_obs_cached(gc_obs, base_obs_loc, base_obs_type, &
my_obs_loc, my_obs_kind, my_obs_type, num_close_obs, close_obs_ind, close_obs_dist, &
ens_handle, last_base_obs_loc, last_num_close_obs, num_close_obs_cached, &
num_close_obs_calls_made)
! set the cutoff default, keep a copy of the original value, and avoid
! looking up the cutoff in a list if the incoming obs is an identity ob
! (and therefore has a negative kind). specific types can never be 0;
! generic kinds (not used here) start their numbering at 0 instead of 1.
if (base_obs_type > 0) then
cutoff_orig = cutoff_list(base_obs_type)
else
cutoff_orig = cutoff
endif
! JLA, could also cache for adaptive_localization which may be expensive?
call adaptive_localization_and_diags(cutoff_orig, cutoff_rev, adaptive_localization_threshold, &
adaptive_cutoff_floor, num_close_obs, close_obs_ind, close_obs_dist, my_obs_type, &
i, base_obs_loc, obs_def, localization_unit)
! Find state variables on my process that are close to observation being assimilated
call get_close_state_cached(gc_state, base_obs_loc, base_obs_type, &
my_state_loc, my_state_kind, my_state_indx, num_close_states, close_state_ind, close_state_dist, &
ens_handle, last_base_states_loc, last_num_close_states, num_close_states_cached, &
num_close_states_calls_made)
!call test_close_obs_dist(close_state_dist, num_close_states, i)
! Loop through to update each of my state variables that is potentially close
STATE_UPDATE: do j = 1, num_close_states
state_index = close_state_ind(j)
if ( allow_missing_in_state ) then
! Don't allow update of state ensemble with any missing values
if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE
endif
! Compute the covariance localization and adjust_obs_impact factors (module storage)
final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_state_loc(state_index), &
my_state_kind(state_index), close_state_dist(j), cutoff_rev)
if(final_factor <= 0.0_r8) cycle STATE_UPDATE
call obs_updates_ens(ens_size, num_groups, ens_handle%copies(1:ens_size, state_index), &
my_state_loc(state_index), my_state_kind(state_index), obs_prior, obs_inc, &
obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, &
net_a, grp_size, grp_beg, grp_end, i, &
my_state_indx(state_index), final_factor, correl, local_varying_ss_inflate, inflate_only)
! Compute spatially-varying state space inflation
if(local_varying_ss_inflate) then
do group = 1, num_groups
call update_varying_state_space_inflation(inflate, &
ens_handle%copies(ENS_INF_COPY, state_index), &
ens_handle%copies(ENS_INF_SD_COPY, state_index), &
ens_handle%copies(ENS_SD_COPY, state_index), &
orig_obs_prior_mean(group), orig_obs_prior_var(group), obs(1), &
obs_err_var, grp_size, final_factor, correl(group), inflate_only)
end do
endif
end do STATE_UPDATE
if(.not. inflate_only) then
! Now everybody updates their obs priors (only ones after this one)
OBS_UPDATE: do j = 1, num_close_obs
obs_index = close_obs_ind(j)
! Only have to update obs that have not yet been used
if(my_obs_indx(obs_index) > i) then
! If forward observation operator failed, no need to update unassimilated observations
if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE
! Compute the covariance localization and adjust_obs_impact factors (module storage)
final_factor = cov_and_impact_factors(base_obs_loc, base_obs_type, my_obs_loc(obs_index), &
my_obs_kind(obs_index), close_obs_dist(j), cutoff_rev)
if(final_factor <= 0.0_r8) cycle OBS_UPDATE
call obs_updates_ens(ens_size, num_groups, obs_ens_handle%copies(1:ens_size, obs_index), &
my_obs_loc(obs_index), my_obs_kind(obs_index), obs_prior, obs_inc, &
obs_prior_mean, obs_prior_var, base_obs_loc, base_obs_type, obs_time, &
net_a, grp_size, grp_beg, grp_end, i, &
-1*my_obs_indx(obs_index), final_factor, correl, .false., inflate_only)
endif
end do OBS_UPDATE
endif
end do SEQUENTIAL_OBS
! Every pe needs to get the current my_inflate and my_inflate_sd back
if(local_single_ss_inflate) then
ens_handle%copies(ENS_INF_COPY, :) = my_inflate
ens_handle%copies(ENS_INF_SD_COPY, :) = my_inflate_sd
end if
! Free up the storage
call destroy_obs(observation)
call get_close_destroy(gc_state)
call get_close_destroy(gc_obs)
! do some stats - being aware that unless we do a reduce() operation
! this is going to be per-task. so only print if something interesting
! shows up in the stats? maybe it would be worth a reduce() call here?
! Assure user we have done something
if (print_trace_details >= 0) then
write(msgstring, '(A,I8,A)') 'Processed', obs_ens_handle%num_vars, ' total observations'
call error_handler(E_MSG,'filter_assim:',msgstring)
endif
! diagnostics for stats on saving calls by remembering obs at the same location.
! change .true. to .false. in the line below to remove the output completely.
if (close_obs_caching) then
if (num_close_obs_cached > 0 .and. do_output()) then
print *, "Total number of calls made to get_close_obs for obs/states: ", &
num_close_obs_calls_made + num_close_states_calls_made
print *, "Total number of calls avoided to get_close_obs for obs/states: ", &
num_close_obs_cached + num_close_states_cached
if (num_close_obs_cached+num_close_obs_calls_made+ &
num_close_states_cached+num_close_states_calls_made > 0) then
print *, "Percent saved: ", 100.0_r8 * &
(real(num_close_obs_cached+num_close_states_cached, r8) / &
(num_close_obs_calls_made+num_close_obs_cached + &
num_close_states_calls_made+num_close_states_cached))
endif
endif
endif
!call test_state_copies(ens_handle, 'end')
! Close the localization diagnostics file
if(output_localization_diagnostics .and. my_task_id() == 0) call close_file(localization_unit)
! get rid of mpi window
call free_mean_window()
! deallocate space
deallocate(close_obs_dist, &
my_obs_indx, &
my_obs_kind, &
my_obs_type, &
close_obs_ind, &
vstatus, &
my_obs_loc)
deallocate(close_state_dist, &
my_state_indx, &
close_state_ind, &
my_state_kind, &
my_state_loc)
! end dealloc
end subroutine filter_assim
!-------------------------------------------------------------
subroutine obs_increment(ens_in, ens_size, obs, obs_var, obs_inc, &
inflate, my_cov_inflate, my_cov_inflate_sd, net_a)
! Given the ensemble prior for an observation, the observation, and
! the observation error variance, computes increments and adjusts
! observation space inflation values
integer, intent(in) :: ens_size
real(r8), intent(in) :: ens_in(ens_size), obs, obs_var
real(r8), intent(out) :: obs_inc(ens_size)
type(adaptive_inflate_type), intent(inout) :: inflate
real(r8), intent(inout) :: my_cov_inflate, my_cov_inflate_sd
real(r8), intent(out) :: net_a
real(r8) :: ens(ens_size), inflate_inc(ens_size)
real(r8) :: prior_mean, prior_var, new_val(ens_size)
integer :: i, ens_index(ens_size), new_index(ens_size)
real(r8) :: rel_weights(ens_size)
! Copy the input ensemble to something that can be modified
ens = ens_in
! Null value of net spread change factor is 1.0
net_a = 0.0_r8
! Compute prior variance and mean from sample
prior_mean = sum(ens) / ens_size
prior_var = sum((ens - prior_mean)**2) / (ens_size - 1)
! If observation space inflation is being done, compute the initial
! increments and update the inflation factor and its standard deviation
! as needed. my_cov_inflate < 0 means don't do any of this.
if(do_obs_inflate(inflate)) then
! If my_cov_inflate_sd is <= 0, just retain current my_cov_inflate setting
if(my_cov_inflate_sd > 0.0_r8) &
! Gamma set to 1.0 because no distance for observation space
call update_inflation(inflate, my_cov_inflate, my_cov_inflate_sd, prior_mean, &
prior_var, ens_size, obs, obs_var, gamma_corr = 1.0_r8)
! Now inflate the ensemble and compute a preliminary inflation increment
call inflate_ens(inflate, ens, prior_mean, my_cov_inflate, prior_var)
! Keep the increment due to inflation alone
inflate_inc = ens - ens_in
! Need to recompute variance if non-deterministic inflation (mean is unchanged)
if(.not. deterministic_inflate(inflate)) &
prior_var = sum((ens - prior_mean)**2) / (ens_size - 1)
endif
! If obs_var == 0, delta function. The mean becomes obs value with no spread.
! If prior_var == 0, obs has no effect. The increments are 0.
! If both obs_var and prior_var == 0 there is no right thing to do, so Stop.
if ((obs_var == 0.0_r8) .and. (prior_var == 0.0_r8)) then
! fail if both obs variance and prior spreads are 0.
write(msgstring, *) 'Observation value is ', obs, ' ensemble mean value is ', prior_mean
write(msgstring2, *) 'The observation has 0.0 error variance, and the ensemble members have 0.0 spread.'
write(msgstring3, *) 'These require inconsistent actions and the algorithm cannot continue.'
call error_handler(E_ERR, 'obs_increment', msgstring, &
source, text2=msgstring2, text3=msgstring3)
else if (obs_var == 0.0_r8) then
! new mean is obs value, so increments are differences between obs
! value and current value. after applying obs, all state will equal obs.
obs_inc(:) = obs - ens
else if (prior_var == 0.0_r8) then
! if all state values are the same, nothing changes.
obs_inc(:) = 0.0_r8
else
! Call the appropriate filter option to compute increments for ensemble
! note that at this point we've taken care of the cases where either the
! obs_var or the prior_var is 0, so the individual routines no longer need
! to have code to test for those cases.
if(filter_kind == 1) then
call obs_increment_eakf(ens, ens_size, prior_mean, prior_var, &
obs, obs_var, obs_inc, net_a)
else if(filter_kind == 2) then
call obs_increment_enkf(ens, ens_size, prior_var, obs, obs_var, obs_inc)
else if(filter_kind == 3) then
call obs_increment_kernel(ens, ens_size, obs, obs_var, obs_inc)
else if(filter_kind == 4) then
call obs_increment_particle(ens, ens_size, obs, obs_var, obs_inc)
else if(filter_kind == 5) then
call obs_increment_ran_kf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
else if(filter_kind == 6) then
call obs_increment_det_kf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
else if(filter_kind == 7) then
call obs_increment_boxcar(ens, ens_size, obs, obs_var, obs_inc, rel_weights)
else if(filter_kind == 8) then
call obs_increment_rank_histogram(ens, ens_size, prior_var, obs, obs_var, obs_inc)
else
call error_handler(E_ERR,'obs_increment', &
'Illegal value of filter_kind in assim_tools namelist [1-8 OK]', source)
endif
endif
! Add in the extra increments if doing observation space covariance inflation
if(do_obs_inflate(inflate)) obs_inc = obs_inc + inflate_inc
! To minimize regression errors, may want to sort to minimize increments
! This makes sense for any of the non-deterministic algorithms
! By doing it here, can take care of both standard non-deterministic updates
! plus non-deterministic obs space covariance inflation. This is expensive, so
! don't use it if it's not needed.
if (sort_obs_inc) then
new_val = ens_in + obs_inc
! Sorting to make increments as small as possible
call index_sort(ens_in, ens_index, ens_size)
call index_sort(new_val, new_index, ens_size)
do i = 1, ens_size
obs_inc(ens_index(i)) = new_val(new_index(i)) - ens_in(ens_index(i))
end do
endif
! Get the net change in spread if obs space inflation was used
if(do_obs_inflate(inflate)) net_a = net_a * sqrt(my_cov_inflate)
end subroutine obs_increment
subroutine obs_increment_eakf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc, a)
!========================================================================
!
! EAKF version of obs increment
integer, intent(in) :: ens_size
real(r8), intent(in) :: ens(ens_size), prior_mean, prior_var, obs, obs_var
real(r8), intent(out) :: obs_inc(ens_size)
real(r8), intent(out) :: a
real(r8) :: new_mean, var_ratio
! Compute the new mean
var_ratio = obs_var / (prior_var + obs_var)
new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
! Compute sd ratio and shift ensemble
a = sqrt(var_ratio)
obs_inc = a * (ens - prior_mean) + new_mean - ens
end subroutine obs_increment_eakf
subroutine obs_increment_ran_kf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
!========================================================================
!
! Forms a random sample of the Gaussian from the update equations.
! This is very close to what a true 'ENSEMBLE' Kalman Filter would
! look like. Note that outliers, multimodality, etc., get tossed.
integer, intent(in) :: ens_size
real(r8), intent(in) :: prior_mean, prior_var
real(r8), intent(in) :: ens(ens_size), obs, obs_var
real(r8), intent(out) :: obs_inc(ens_size)
real(r8) :: new_mean, var_ratio
real(r8) :: temp_mean, temp_var, new_ens(ens_size), new_var
integer :: i
var_ratio = obs_var / (prior_var + obs_var)
new_var = var_ratio * prior_var
new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
! This will reproduce exactly for multiple runs with the same task count,
! but WILL NOT reproduce for a different number of MPI tasks.
! To make it independent of the number of MPI tasks, it would need to
! use the global ensemble number or something else that remains constant
! as the processor count changes. this is not currently an argument to
! this function and so we are not trying to make it task-count invariant.
! Form a random sample from the updated distribution
! Then adjust the mean (what about adjusting the variance?)!
! Definitely need to sort with this; sort is done in main obs_increment