47 integer,
public,
parameter ::
nv3d = 11
48 integer,
public,
parameter ::
nv2d = 0
49 integer,
public,
parameter ::
nid_obs = 16
50 integer,
public,
parameter ::
nobtype = 24
59 integer,
allocatable :: elm(:)
60 real(
rp),
allocatable :: lon(:)
61 real(
rp),
allocatable :: lat(:)
62 real(
rp),
allocatable :: lev(:)
63 real(
rp),
allocatable :: dat(:)
64 real(
rp),
allocatable :: err(:)
65 integer,
allocatable :: typ(:)
66 real(
rp),
allocatable :: dif(:)
68 real(
rp),
allocatable :: ri(:)
69 real(
rp),
allocatable :: rj(:)
70 integer,
allocatable :: rank(:)
75 integer :: nobs_in_key = 0
76 integer,
allocatable :: set(:)
77 integer,
allocatable :: idx(:)
78 integer,
allocatable :: key(:)
79 real(
rp),
allocatable :: val(:)
80 real(
rp),
allocatable :: ensval(:,:)
81 real(
rp),
allocatable :: eqv(:,:)
82 real(
rp),
allocatable :: qv(:)
83 real(
rp),
allocatable :: tm(:)
84 real(
rp),
allocatable :: pm(:)
85 integer,
allocatable :: qc(:)
97 integer,
allocatable :: n(:,:,:)
98 integer,
allocatable :: ac(:,:,:)
99 integer,
allocatable :: tot(:)
100 integer,
allocatable :: n_ext(:,:)
101 integer,
allocatable :: ac_ext(:,:)
105 integer,
allocatable :: next(:,:)
112 private :: letkf_core
114 private :: get_nobs_radar
116 private :: read_obs_radar
117 private :: read_obs_radar_toshiba_pawr
118 private :: read_obs_radar_toshiba_mp_pawr
119 private :: str_replace
121 private :: radar_georeference
122 private :: define_grid
123 private :: calc_ref_vr
125 private :: obs_local_range
126 private :: obs_local_cal
127 private :: relax_beta
128 private :: weight_rtpp
129 private :: weight_rtps
130 private :: com_distll_1
131 private :: ensmean_grd
132 private :: trans_xtoy
133 private :: trans_xtoy_radar
139 private :: itpl_2d_column
141 private :: merge_sort_parallel
142 private :: merge_sort_2threads
143 private :: merge_2threads
145 private :: merge_sort_mpi
147 private :: merge_mpi_no_nest
148 private :: scatter_grd_mpi
149 private :: scatter_grd_mpi_all2all
150 private :: gather_grd_mpi_all2all
151 private :: grd_to_buf
152 private :: buf_to_grd
153 private :: calc_z_grd
154 private :: qc_indexing_and_packing
156 private :: uid_obs_varlocal
157 private :: binary_search_i8
158 private :: rank_1d_2d
159 private :: rank_2d_1d
163 private :: ij_obsgrd_ext
164 private :: obs_choose
165 private :: obs_choose_ext
166 private :: obs_info_allocate
167 private :: obs_info_deallocate
168 private :: obs_da_value_allocate
169 private :: obs_da_value_deallocate
170 private :: obs_da_value_allreduce
171 private :: obs_da_value_partial_reduce_iter
172 private :: read_ens_mpi_addiinfl
174 private :: monit_obs_mpi
176 private :: monit_print
177 private :: state_to_history
184 integer,
parameter :: nobsfilemax = 10
185 integer,
parameter :: membermax = 10000
187 integer,
parameter :: iv3d_rho = 1
188 integer,
parameter :: iv3d_rhou = 2
189 integer,
parameter :: iv3d_rhov = 3
190 integer,
parameter :: iv3d_rhow = 4
191 integer,
parameter :: iv3d_rhot = 5
192 integer,
parameter :: iv3d_u = 1
193 integer,
parameter :: iv3d_v = 2
194 integer,
parameter :: iv3d_w = 3
195 integer,
parameter :: iv3d_t = 4
196 integer,
parameter :: iv3d_p = 5
197 integer,
parameter :: iv3d_q = 6
198 integer,
parameter :: iv3d_qc = 7
199 integer,
parameter :: iv3d_qr = 8
200 integer,
parameter :: iv3d_qi = 9
201 integer,
parameter :: iv3d_qs = 10
202 integer,
parameter :: iv3d_qg = 11
205 integer,
parameter :: nv3dd = 13
206 integer,
parameter :: nv2dd = 7
207 integer,
parameter :: iv3dd_u = 1
208 integer,
parameter :: iv3dd_v = 2
209 integer,
parameter :: iv3dd_w = 3
210 integer,
parameter :: iv3dd_t = 4
211 integer,
parameter :: iv3dd_p = 5
212 integer,
parameter :: iv3dd_q = 6
213 integer,
parameter :: iv3dd_qc = 7
214 integer,
parameter :: iv3dd_qr = 8
215 integer,
parameter :: iv3dd_qi = 9
216 integer,
parameter :: iv3dd_qs = 10
217 integer,
parameter :: iv3dd_qg = 11
218 integer,
parameter :: iv3dd_rh = 12
219 integer,
parameter :: iv3dd_hgt = 13
220 integer,
parameter :: iv2dd_topo = 1
221 integer,
parameter :: iv2dd_ps = 2
222 integer,
parameter :: iv2dd_rain = 3
223 integer,
parameter :: iv2dd_u10m = 4
224 integer,
parameter :: iv2dd_v10m = 5
225 integer,
parameter :: iv2dd_t2m = 6
226 integer,
parameter :: iv2dd_q2m = 7
228 integer,
parameter :: iqc_good = 0
229 integer,
parameter :: iqc_gross_err = 5
230 integer,
parameter :: iqc_ps_ter = 10
231 integer,
parameter :: iqc_ref_low = 11
232 integer,
parameter :: iqc_ref_mem = 12
233 integer,
parameter :: iqc_radar_vhi = 19
234 integer,
parameter :: iqc_out_vhi = 20
235 integer,
parameter :: iqc_out_vlo = 21
236 integer,
parameter :: iqc_obs_bad = 50
237 integer,
parameter :: iqc_otype = 90
238 integer,
parameter :: iqc_time = 97
239 integer,
parameter :: iqc_out_h = 98
240 integer,
parameter :: iqc_undef = 99
242 integer,
parameter :: nid_obs_varlocal = 8
246 integer,
parameter :: id_u_obs = 2819
247 integer,
parameter :: id_v_obs = 2820
248 integer,
parameter :: id_t_obs = 3073
249 integer,
parameter :: id_tv_obs = 3074
250 integer,
parameter :: id_q_obs = 3330
251 integer,
parameter :: id_rh_obs = 3331
255 integer,
parameter :: id_ps_obs = 14593
256 integer,
parameter :: id_rain_obs = 19999
257 integer,
parameter :: id_tclon_obs = 99991
258 integer,
parameter :: id_tclat_obs = 99992
259 integer,
parameter :: id_tcmip_obs = 99993
263 integer,
parameter :: id_radar_ref_obs = 4001
264 integer,
parameter :: id_radar_ref_zero_obs = 4004
265 integer,
parameter :: id_radar_vr_obs = 4002
266 integer,
parameter :: id_radar_prh_obs = 4003
270 integer,
parameter :: id_h08ir_obs = 8800
272 real(
rp),
parameter :: dist_zero_fac = 3.651483717
273 real(
rp),
parameter :: dist_zero_fac_square = 13.33333333
274 real(
rp),
parameter :: minz = 0.01_rp
275 real(
rp),
parameter :: vr_min_dist = 8000.0_rp
277 integer,
parameter :: elem_uid(
nid_obs) = &
278 (/ id_u_obs, id_v_obs, id_t_obs, id_tv_obs, id_q_obs, &
279 id_rh_obs, id_ps_obs, id_rain_obs, id_radar_ref_obs, id_radar_ref_zero_obs, &
280 id_radar_vr_obs, id_radar_prh_obs, id_h08ir_obs, id_tclon_obs, id_tclat_obs, &
283 character(3),
parameter :: obelmlist(
nid_obs) = &
284 (/
' U',
' V',
' T',
' Tv',
' Q', &
285 ' RH',
' PS',
'PRC',
'REF',
'RE0', &
286 ' Vr',
'PRH',
'H08',
'TCX',
'TCY', &
289 character(3),
parameter :: obelmlist_varlocal(nid_obs_varlocal) = &
290 (/
'WND',
' T',
'MOI',
' PS',
'PRC', &
291 'TCV',
'REF',
' Vr' /)
293 character(6),
parameter :: obtypelist(
nobtype) = &
294 (/
'ADPUPA',
'AIRCAR',
'AIRCFT',
'SATWND',
'PROFLR', &
295 'VADWND',
'SATEMP',
'ADPSFC',
'SFCSHP',
'SFCBOG', &
296 'SPSSMI',
'SYNDAT',
'ERS1DA',
'GOESND',
'QKSWND', &
297 'MSONET',
'GPSIPW',
'RASSDA',
'WDSATR',
'ASCATW', &
298 'TMPAPR',
'PHARAD',
'H08IRB',
'TCVITL' /)
300 logical :: ens_with_mdet = .false.
302 logical :: use_obs(
nobtype) = .true.
303 logical :: obs_postfix_timelabel = .false.
304 logical :: obsda_run(nobsfilemax) = .true.
305 logical :: obsda_out = .false.
306 character(len=H_LONG) :: obs_in_name(nobsfilemax) =
''
307 character(len=H_LONG) :: obsda_out_basename =
''
308 character(len=H_LONG) :: obsda_mean_out_basename =
''
309 character(len=H_LONG) :: obsda_mdet_out_basename =
''
311 integer :: slot_start = 1
312 integer :: slot_end = 1
313 integer :: slot_base = 1
314 real(
rp) :: slot_tinterval = 3600.0_rp
316 logical :: departure_stat_radar = .false.
317 real(
rp) :: departure_stat_t_range = 0.0
325 real(
rp) :: var_local_rain(
nv3d+
nv2d) = 1.0_rp
327 real(
rp) :: var_local_radar_ref(
nv3d+
nv2d) = 1.0_rp
328 real(
rp) :: var_local_radar_vr(
nv3d+
nv2d) = 1.0_rp
330 character(len=H_LONG) :: infl_mul_in_basename =
''
331 character(len=H_LONG) :: infl_mul_out_basename =
''
332 real(
rp) :: infl_mul = 1.0_rp
334 real(
rp) :: infl_mul_min = -1.0_rp
335 logical :: infl_mul_adaptive = .false.
337 character(len=H_LONG) :: infl_add_in_basename =
''
338 real(
rp) :: infl_add = 0.0_rp
339 logical :: infl_add_shuffle = .false.
340 logical :: infl_add_q_ratio = .false.
341 logical :: infl_add_ref_only = .false.
344 real(
rp) :: obserr_u = 1.0_rp
345 real(
rp) :: obserr_v = 1.0_rp
346 real(
rp) :: obserr_t = 1.0_rp
347 real(
rp) :: obserr_q = 0.001_rp
348 real(
rp) :: obserr_rh = 10.0_rp
349 real(
rp) :: obserr_ps = 100.0_rp
350 real(
rp) :: obserr_radar_ref = 5.0_rp
351 real(
rp) :: obserr_radar_vr = 3.0_rp
352 real(
rp) :: obserr_tcx = 50.e+3_rp
353 real(
rp) :: obserr_tcy = 50.e+3_rp
354 real(
rp) :: obserr_tcp = 5.e+2_rp
355 real(
rp) :: obserr_pq = 0.001_rp
358 integer :: nradartype = 1
359 real(
rp) :: radar_so_size_hori = 1000.0_rp
360 real(
rp) :: radar_so_size_vert = 1000.0_rp
361 real(
rp) :: radar_max_abs_vr = 100.0_rp
362 integer :: radar_thin_letkf_method = 0
368 integer :: radar_thin_letkf_hgrid = 1
369 integer :: radar_thin_letkf_vgrid = 1
370 integer :: radar_thin_letkf_hnear = 1
371 integer :: radar_thin_letkf_vnear = 1
372 integer :: radar_thin_hori = 1
373 integer :: radar_thin_vert = 1
374 logical :: radar_use_vr_std = .true.
375 logical :: radar_bias_cor_rain = .false.
376 logical :: radar_bias_cor_clr = .false.
377 real(
rp) :: radar_bias_rain_const_dbz = 0.0_rp
378 real(
rp) :: radar_bias_clr_const_dbz = 0.0_rp
379 real(
rp) :: vr_std_threshold = 2.5_rp
380 real(
rp) :: attenuation_threshold = 0.25_rp
382 logical :: use_radar_ref = .true.
383 logical :: use_radar_vr = .true.
384 logical :: use_radar_pseudo_rh = .false.
386 logical :: use_obserr_radar_ref = .false.
387 logical :: use_obserr_radar_vr = .false.
388 logical :: radar_obs_4d = .false.
389 logical :: radar_pqv = .false.
391 real(
rp) :: radar_pqv_omb = 25.0_rp
392 real(
rp) :: radar_ref_thres_dbz = 15.0_rp
393 integer :: min_radar_ref_member_obsrain = 1
394 integer :: min_radar_ref_member_obsnorain = 1
396 real(
rp) :: min_radar_ref_dbz = 0.0_rp
397 real(
rp) :: min_radar_ref_dbz_vr = 5.0_rp
398 real(
rp) :: min_radar_ref_vr = 0.0_rp
399 real(
rp) :: low_ref_shift = 0.0_rp
401 real(
rp) :: radar_zmax = 99.e+3_rp
402 real(
rp) :: radar_zmin = -99.e+3_rp
403 real(
rp) :: radar_prh_error = 0.1_rp
405 integer :: interpolation_technique = 1
406 integer :: method_ref_calc = 3
407 logical :: use_method3_ref_melt = .false.
408 logical :: use_t08_rs2014 = .false.
409 logical :: use_terminal_velocity = .false.
410 logical :: use_attenuation = .true.
411 logical :: use_qcflag = .true.
413 real(
rp) :: relax_alpha = 0.0_rp
414 real(
rp) :: relax_alpha_spread = 0.0_rp
415 logical :: relax_to_inflated_prior = .false.
417 logical :: relax_spread_out = .false.
418 character(len=H_LONG) :: relax_spread_out_basename =
''
420 real(
rp) :: gross_error = 5.0_rp
421 real(
rp) :: gross_error_rain = -1.0_rp
422 real(
rp) :: gross_error_radar_ref = -1.0_rp
423 real(
rp) :: gross_error_radar_vr = -1.0_rp
424 real(
rp) :: gross_error_radar_prh = -1.0_rp
425 real(
rp) :: gross_error_tcx = -1.0_rp
426 real(
rp) :: gross_error_tcy = -1.0_rp
427 real(
rp) :: gross_error_tcp = -1.0_rp
429 real(
rp) :: q_update_top = 0.0_rp
430 real(
rp) :: q_sprd_max = -1.0_rp
432 real(
rp) :: boundary_buffer_width = 0.0_rp
433 real(
rp) :: ps_adjust_thres = 100.0_rp
435 logical :: nobs_out = .false.
436 character(len=H_LONG) :: nobs_out_basename =
''
439 (/ 0.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
440 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
441 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
442 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
443 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp /)
445 (/ 0.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
446 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
447 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
448 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
449 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp /)
451 (/ 0.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
452 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
453 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
454 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
455 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp /)
457 real(
rp) :: hori_local_radar_obsnoref = -1.0_rp
458 real(
rp) :: hori_local_radar_vr = -1.0_rp
459 real(
rp) :: vert_local_radar_vr = -1.0_rp
460 real(
rp) :: vert_local_rain_base = 85000.0_rp
462 integer :: max_nobs_per_grid(
nobtype) = &
463 (/ 0, -1, -1, -1, -1, -1, -1, -1, -1, -1, &
464 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, &
467 integer :: max_nobs_per_grid_criterion = 1
472 (/ 0.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
473 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
474 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
475 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
476 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp /)
478 real(
rp) :: obs_sort_grid_spacing(
nobtype) = &
479 (/ 0.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
480 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
481 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
482 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp, &
483 -1.0_rp, -1.0_rp, -1.0_rp, -1.0_rp /)
485 integer :: utime_obs(6) = &
486 (/ -1, -1, -1, -1, -1, -1 /)
488 integer :: mmean = -99
489 integer :: mmdet = -99
490 integer :: mmdetin = -99
491 integer :: mmdetobs = -99
493 logical :: letkf_entire_grid_search_x = .false.
494 logical :: letkf_entire_grid_search_y = .false.
496 logical :: letkf_debug_log = .false.
498 type(
obs_info),
allocatable :: obs(:)
505 real(
rp) :: min_radar_ref
506 real(
rp) :: radar_ref_thres
507 real(
rp) :: radar_bias_rain_const
508 real(
rp) :: radar_bias_clr_const
510 real(
rp),
allocatable :: var_local(:,:)
512 integer,
allocatable :: n_merge(:)
513 integer,
allocatable :: ic_merge(:,:)
515 integer :: nobstotalg
517 integer :: maxnobs_per_ctype
522 integer,
allocatable :: elm_ctype(:)
523 integer,
allocatable :: elm_u_ctype(:)
524 integer,
allocatable :: typ_ctype(:)
525 real(
rp),
allocatable :: hori_loc_ctype(:)
526 real(
rp),
allocatable :: vert_loc_ctype(:)
529 integer :: obsdep_nobs
530 integer,
allocatable :: obsdep_set(:)
531 integer,
allocatable :: obsdep_idx(:)
532 integer,
allocatable :: obsdep_qc (:)
533 real(
rp),
allocatable :: obsdep_omb(:)
534 real(
rp),
allocatable :: obsdep_oma(:)
537 integer :: n_merge_max
538 logical :: radar_only
577 integer,
allocatable :: nij1node(:)
579 real(
rp),
allocatable :: rig1(:)
580 real(
rp),
allocatable :: rjg1(:)
581 real(
rp),
allocatable :: topo1(:)
582 real(
rp),
allocatable :: hgt1(:,:)
583 real(
rp),
allocatable :: v3dg(:,:,:,:)
584 real(
rp),
allocatable :: v2dg(:,:,:)
585 real(
rp),
allocatable :: v3d(:,:,:)
586 real(
rp),
allocatable :: v2d(:,:)
587 real(
rp),
allocatable :: topo2d(:,:)
622 integer,
intent(in) :: obs_in_num
623 integer,
intent(in) :: ensemble_comm
624 integer,
intent(in) :: ensemble_nprocs
625 integer,
intent(in) :: ensemble_myrank
626 integer,
intent(in) :: local_comm
627 integer,
intent(in) :: local_nprocs
628 integer,
intent(in) :: local_myrank
629 integer,
intent(in) :: prc_num_x
630 integer,
intent(in) :: prc_num_y
631 integer,
intent(in) :: ka, ks, ke
632 integer,
intent(in) :: ia, is, ie
633 integer,
intent(in) :: ja, js, je
634 integer,
intent(in) :: kmax
635 integer,
intent(in) :: imax
636 integer,
intent(in) :: jmax
637 integer,
intent(in) :: khalo
638 integer,
intent(in) :: ihalo
639 integer,
intent(in) :: jhalo
641 real(
rp),
intent(in) :: delta_x
642 real(
rp),
intent(in) :: delta_y
643 real(
rp),
intent(in) :: zsfc(:,:)
645 logical :: letkf_deterministic_run
647 integer :: letkf_mem_nodes = 0
649 namelist / param_letkf / &
651 letkf_deterministic_run, &
657 departure_stat_radar, &
666 relax_alpha_spread, &
667 relax_to_inflated_prior, &
670 gross_error_radar_ref, &
671 gross_error_radar_vr, &
672 gross_error_radar_prh, &
678 boundary_buffer_width, &
682 hori_local_radar_obsnoref, &
683 hori_local_radar_vr, &
684 vert_local_radar_vr, &
686 max_nobs_per_grid_criterion, &
688 obs_sort_grid_spacing, &
692 use_terminal_velocity, &
693 use_obserr_radar_ref, &
694 use_obserr_radar_vr, &
695 radar_ref_thres_dbz, &
696 min_radar_ref_member_obsrain, &
697 min_radar_ref_member_obsnorain, &
699 min_radar_ref_dbz_vr, &
703 radar_so_size_hori, &
704 radar_so_size_vert, &
706 use_method3_ref_melt, &
707 radar_bias_cor_rain, &
708 radar_bias_cor_clr, &
709 radar_bias_rain_const_dbz, &
710 radar_bias_clr_const_dbz, &
711 radar_thin_letkf_method, &
712 radar_thin_letkf_hgrid, &
713 radar_thin_letkf_vgrid, &
714 radar_thin_letkf_hnear, &
715 radar_thin_letkf_vnear, &
727 letkf_entire_grid_search_x, &
728 letkf_entire_grid_search_y
733 integer :: i, j, k, n
738 log_info(
"LETKF_setup",*)
'Setup'
742 else if(
rp ==
dp )
then
743 datatype = mpi_double_precision
745 log_error(
"obsope_tool",*)
'The precision has not been implemented yet:',
rp
749 letkf_deterministic_run = .false.
755 log_info(
"LETKF_setup",*)
'Not found namelist. Default used.'
756 elseif( ierr > 0 )
then
757 log_error(
"LETKF_setup",*)
'Not appropriate names in namelist PARAM_LETKF. Check!'
764 nens = ensemble_nprocs
765 nensobs = ensemble_nprocs
766 nmem = ensemble_nprocs
773 nlong = imax * prc_num_x
774 nlatg = jmax * prc_num_y
787 comm_ens = ensemble_comm
788 nprc_ens = ensemble_nprocs
789 rank_ens = ensemble_myrank
790 comm_lcl = local_comm
791 nprc_lcl = local_nprocs
792 rank_lcl = local_myrank
797 if( letkf_deterministic_run )
then
798 ens_with_mdet = .true.
806 if( letkf_mem_nodes == 0 )
then
809 if( letkf_mem_nodes > 1 )
then
810 n_mem = nprc_ens / letkf_mem_nodes
816 nitmax = ( nprc_ens - 1 ) / ( n_mem * n_mempn ) + 1
818 i = mod( nlon*nlat, nprc_ens )
819 nij1max = ( nlon*nlat - i ) / nprc_ens + 1
820 if( rank_ens < i )
then
826 allocate( obs( obs_in_num ) )
827 allocate( nij1node( nprc_ens ) )
831 nij1node(n) = nij1max
833 nij1node(n) = nij1max - 1
837 allocate( rig1( nij1 ) )
838 allocate( rjg1( nij1 ) )
839 allocate( topo1( nij1 ) )
840 allocate( hgt1( nij1, nlev ) )
842 allocate( v3dg( nlev, nlon, nlat,
nv3d ) )
843 allocate( v2dg( nlon, nlat,
nv2d ) )
845 allocate( v3d( nij1, nlev,
nv3d ) )
846 allocate( v2d( nij1,
nv2d ) )
850 v3dg(1,i,j,1) = real( i +
prc_2drank(local_myrank,1) * imax + ihalo, kind=
rp )
851 v3dg(1,i,j,2) = real( j +
prc_2drank(local_myrank,2) * jmax + jhalo, kind=
rp )
852 v3dg(1,i,j,3) = zsfc(i+ihalo,j+jhalo)
856 call scatter_grd_mpi( mod( nens, n_mem*n_mempn ), v3dg, v2dg, v3d, v2d )
860 topo1(:) = v3d(:,1,3)
862 call calc_z_grd( nij1, topo1, hgt1 )
864 allocate( var_local(
nv3d+
nv2d, nid_obs_varlocal ) )
866 if( radar_ref_thres_dbz < min_radar_ref_dbz )
then
867 radar_ref_thres_dbz = min_radar_ref_dbz
870 min_radar_ref = 10.0_rp ** ( min_radar_ref_dbz / 10.0_rp )
871 min_radar_ref_vr = 10.0_rp ** ( min_radar_ref_dbz_vr / 10.0_rp )
872 radar_ref_thres = 10.0_rp ** ( radar_ref_thres_dbz / 10.0_rp )
873 radar_bias_rain_const = 10.0_rp ** ( radar_bias_rain_const_dbz / 10.0_rp )
874 radar_bias_clr_const = 10.0_rp ** ( radar_bias_clr_const_dbz / 10.0_rp )
892 deallocate( var_local )
895 deallocate( nij1node )
912 integer,
intent(in) :: obs_in_num
914 character(len=H_LONG),
intent(in) :: obs_in_format(:)
915 character(len=H_LONG),
intent(in) :: obs_in_basename(:)
916 character(len=H_LONG),
intent(in) :: obs_in_maskfile
918 character(len=H_LONG) :: obsfile
919 character(len=H_LONG) :: timelabel_obsfile
927 log_progress(*)
'data-assimilation / LETKF / obs / readfile'
929 if( rank_ens == 0 .and. rank_lcl == 0 )
then
930 timelabel_obsfile =
'_???????????????????.dat'
934 obsfile = trim(obs_in_basename(n))//trim(timelabel_obsfile)
936 if (obs_in_format(n) /=
'PAWR_TOSHIBA' .and. &
937 obs_in_format(n) /=
'MP_PAWR_TOSHIBA' )
then
938 inquire( file=obsfile, exist=err )
940 log_info(
"LETKF_obs_readfile",*)
'Warning: File (',trim(obsfile),
') is not found. Skip.'
943 call obs_info_allocate(obs(n), extended=.true.)
948 select case( obs_in_format(n) )
950 call get_nobs( obsfile, 8, obs(n)%nobs )
951 call obs_info_allocate( obs(n), extended=.true. )
952 call read_obs( obsfile, obs(n) )
954 call get_nobs_radar( obsfile, obs(n)%nobs, obs(n)%meta(1), obs(n)%meta(2), obs(n)%meta(3) )
955 call obs_info_allocate( obs(n), extended=.true. )
956 call read_obs_radar( obsfile, obs(n) )
958 log_error(
"LETKF_obs_readfile",*)
'Error: This system has not been implemented yet. (OBS_IN_FORMAT(:) = PAWR_JRC)'
961 case (
'PAWR_TOSHIBA' )
962 call read_obs_radar_toshiba_pawr( obs(n), obsfile )
963 case (
'MP_PAWR_TOSHIBA' )
964 call read_obs_radar_toshiba_mp_pawr( obs(n), obsfile, obs_in_maskfile )
966 log_error(
"LETKF_obs_readfile",*)
'Error: Unsupported observation file format.'
974 if( rank_lcl == 0 )
then
975 call mpi_bcast( obs(n)%nobs, 1, mpi_integer, 0, comm_ens, ierr )
977 if( rank_ens /= 0 )
then
978 call obs_info_allocate( obs(n), extended=.true. )
981 call mpi_bcast( obs(n)%elm, obs(n)%nobs, mpi_integer, 0, comm_ens, ierr )
982 call mpi_bcast( obs(n)%lon, obs(n)%nobs, datatype, 0, comm_ens, ierr )
983 call mpi_bcast( obs(n)%lat, obs(n)%nobs, datatype, 0, comm_ens, ierr )
984 call mpi_bcast( obs(n)%lev, obs(n)%nobs, datatype, 0, comm_ens, ierr )
985 call mpi_bcast( obs(n)%dat, obs(n)%nobs, datatype, 0, comm_ens, ierr )
986 call mpi_bcast( obs(n)%err, obs(n)%nobs, datatype, 0, comm_ens, ierr )
987 call mpi_bcast( obs(n)%typ, obs(n)%nobs, mpi_integer, 0, comm_ens, ierr )
988 call mpi_bcast( obs(n)%dif, obs(n)%nobs, datatype, 0, comm_ens, ierr )
993 call mpi_bcast( obs(n)%nobs, 1, mpi_integer, 0, comm_lcl, ierr )
995 if( rank_lcl /= 0 )
then
996 call obs_info_allocate( obs(n), extended=.true. )
999 call mpi_bcast( obs(n)%elm, obs(n)%nobs, mpi_integer, 0, comm_lcl, ierr )
1000 call mpi_bcast( obs(n)%lon, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1001 call mpi_bcast( obs(n)%lat, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1002 call mpi_bcast( obs(n)%lev, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1003 call mpi_bcast( obs(n)%dat, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1004 call mpi_bcast( obs(n)%err, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1005 call mpi_bcast( obs(n)%typ, obs(n)%nobs, mpi_integer, 0, comm_lcl, ierr )
1006 call mpi_bcast( obs(n)%dif, obs(n)%nobs, datatype, 0, comm_lcl, ierr )
1009 log_info(
'LETKF_obs_readfile',*)
'observations input:', obs(n)%nobs
1018 integer,
intent(in) :: obs_in_num
1023 log_progress(*)
'data-assimilation / LETKF / cleaning'
1025 do n = 1, obs_in_num
1026 call obs_info_deallocate( obs(n) )
1029 call obs_da_value_deallocate( obsda_sort )
1038 OBS_IN_NUM, & ! [IN]
1039 OBS_IN_FORMAT, & ! [IN]
1063 integer,
intent(in) :: obs_in_num
1064 character(len=H_LONG),
intent(in) :: obs_in_format(:)
1066 real(
rp),
intent(in) :: u (nlevh,nlonh,nlath)
1067 real(
rp),
intent(in) :: v (nlevh,nlonh,nlath)
1068 real(
rp),
intent(in) :: w (nlevh,nlonh,nlath)
1069 real(
rp),
intent(in) :: temp(nlevh,nlonh,nlath)
1070 real(
rp),
intent(in) :: pres(nlevh,nlonh,nlath)
1071 real(
rp),
intent(in) :: qv (nlevh,nlonh,nlath)
1072 real(
rp),
intent(in) :: qc (nlevh,nlonh,nlath)
1073 real(
rp),
intent(in) :: qr (nlevh,nlonh,nlath)
1074 real(
rp),
intent(in) :: qi (nlevh,nlonh,nlath)
1075 real(
rp),
intent(in) :: qs (nlevh,nlonh,nlath)
1076 real(
rp),
intent(in) :: qg (nlevh,nlonh,nlath)
1077 real(
rp),
intent(in) :: rh (nlevh,nlonh,nlath)
1078 real(
rp),
intent(in) :: hgt (nlevh,nlonh,nlath)
1080 real(
rp),
intent(in) :: topo(nlonh,nlath)
1081 real(
rp),
intent(in) :: ps (nlonh,nlath)
1082 real(
rp),
intent(in) :: rain(nlonh,nlath)
1083 real(
rp),
intent(in) :: u10m(nlonh,nlath)
1084 real(
rp),
intent(in) :: v10m(nlonh,nlath)
1085 real(
rp),
intent(in) :: t2m (nlonh,nlath)
1086 real(
rp),
intent(in) :: q2m (nlonh,nlath)
1088 integer,
optional,
intent(in) :: nobs_extern
1093 integer :: it, im, iof, islot, ierr
1094 integer :: n, nn, nsub, nmod, n1, n2
1098 integer :: nobs_max_per_file
1099 integer :: nobs_max_per_file_sub
1100 integer :: slot_nobsg
1102 integer :: ip, ibufs
1103 integer,
allocatable :: cntr(:), dspr(:)
1104 integer,
allocatable :: cnts(:), dsps(:)
1105 integer,
allocatable :: bsn(:,:), bsna(:,:), bsnext(:,:)
1106 integer :: islot_time_out, islot_domain_out
1108 integer,
allocatable :: obrank_bufs(:)
1109 real(
rp),
allocatable :: ri_bufs(:)
1110 real(
rp),
allocatable :: rj_bufs(:)
1112 integer,
allocatable :: obset_bufs(:)
1113 integer,
allocatable :: obidx_bufs(:)
1115 integer :: slot_id(slot_start:slot_end)
1116 real(
rp) :: slot_lb(slot_start:slot_end)
1117 real(
rp) :: slot_ub(slot_start:slot_end)
1119 real(
rp),
allocatable :: v3dg(:,:,:,:)
1120 real(
rp),
allocatable :: v2dg(:,:,:)
1122 real(
rp) :: ril, rjl, rk, rkz
1124 character(len=4) :: nstr
1127 log_progress(*)
'data-assimilation / LETKF / obs / operator'
1134 nobs_max_per_file = 0
1135 do iof = 1, obs_in_num
1136 if (obs(iof)%nobs > nobs_max_per_file)
then
1137 nobs_max_per_file = obs(iof)%nobs
1139 if (obsda_run(iof))
then
1140 nobs_all = nobs_all + obs(iof)%nobs
1144 nobs_max_per_file_sub = (nobs_max_per_file - 1) / nprc_lcl + 1
1145 allocate( obrank_bufs(nobs_max_per_file_sub) )
1146 allocate( ri_bufs(nobs_max_per_file_sub) )
1147 allocate( rj_bufs(nobs_max_per_file_sub) )
1149 allocate( cntr(nprc_lcl) )
1150 allocate( dspr(nprc_lcl) )
1156 do iof = 1, obs_in_num
1157 if( obs(iof)%nobs > 0 )
then
1159 nsub = obs(iof)%nobs / nprc_lcl
1160 nmod = mod( obs(iof)%nobs, nprc_lcl )
1164 do ip = nmod+1, nprc_lcl
1169 dspr(ip) = dspr(ip-1) + cntr(ip-1)
1173 do ibufs = 1, cntr(rank_lcl+1)
1174 n = dspr(rank_lcl+1) + ibufs
1176 call phys2ij(obs(iof)%lon(n), obs(iof)%lat(n), ri_bufs(ibufs), rj_bufs(ibufs))
1177 call rij_rank(ri_bufs(ibufs), rj_bufs(ibufs), obrank_bufs(ibufs))
1180 call mpi_allgatherv( obrank_bufs, cntr(rank_lcl+1), mpi_integer, obs(iof)%rank, cntr, dspr, mpi_integer, comm_lcl, ierr )
1181 call mpi_allgatherv( ri_bufs, cntr(rank_lcl+1), datatype, obs(iof)%ri, cntr, dspr, datatype, comm_lcl, ierr )
1182 call mpi_allgatherv( rj_bufs, cntr(rank_lcl+1), datatype, obs(iof)%rj, cntr, dspr, datatype, comm_lcl, ierr )
1187 deallocate(cntr, dspr)
1188 deallocate(obrank_bufs, ri_bufs, rj_bufs)
1193 islot_time_out = slot_end + 1
1194 islot_domain_out = slot_end + 2
1196 allocate( bsn( slot_start :slot_end+2, 0:nprc_lcl-1 ) )
1197 allocate( bsna( slot_start-1:slot_end+2, 0:nprc_lcl-1 ) )
1199 if (rank_ens == 0)
then
1200 allocate ( obset_bufs(nobs_all) )
1201 allocate ( obidx_bufs(nobs_all) )
1204 if (rank_ens == 0 .and. rank_lcl == 0)
then
1205 allocate( bsnext( slot_start:slot_end+2, 0:nprc_lcl-1 ) )
1210 do iof = 1, obs_in_num
1211 if (obsda_run(iof) .and. obs(iof)%nobs > 0)
then
1212 do n = 1, obs(iof)%nobs
1213 if (obs(iof)%rank(n) == -1)
then
1215 bsn(islot_domain_out, 0) = bsn(islot_domain_out, 0) + 1
1217 islot = ceiling(obs(iof)%dif(n) / slot_tinterval - 0.5d0) + slot_base
1218 if (islot < slot_start .or. islot > slot_end)
then
1219 islot = islot_time_out
1221 bsn(islot, obs(iof)%rank(n)) = bsn(islot, obs(iof)%rank(n)) + 1
1227 do ip = 0, nprc_lcl-1
1229 bsna(slot_start-1, ip) = bsna(slot_end+2, ip-1)
1231 do islot = slot_start, slot_end+2
1232 bsna(islot, ip) = bsna(islot-1, ip) + bsn(islot, ip)
1234 bsnext(slot_start:slot_end+2, ip) = bsna(slot_start-1:slot_end+1, ip)
1237 do iof = 1, obs_in_num
1238 if (obsda_run(iof) .and. obs(iof)%nobs > 0)
then
1239 do n = 1, obs(iof)%nobs
1240 if (obs(iof)%rank(n) == -1)
then
1242 bsnext(islot_domain_out, 0) = bsnext(islot_domain_out, 0) + 1
1243 obset_bufs(bsnext(islot_domain_out, 0)) = iof
1244 obidx_bufs(bsnext(islot_domain_out, 0)) = n
1246 islot = ceiling(obs(iof)%dif(n) / slot_tinterval - 0.5d0) + slot_base
1247 if (islot < slot_start .or. islot > slot_end)
then
1248 islot = islot_time_out
1250 bsnext(islot, obs(iof)%rank(n)) = bsnext(islot, obs(iof)%rank(n)) + 1
1251 obset_bufs(bsnext(islot, obs(iof)%rank(n))) = iof
1252 obidx_bufs(bsnext(islot, obs(iof)%rank(n))) = n
1258 deallocate( bsnext )
1265 if( rank_lcl == 0 )
then
1266 call mpi_bcast( bsn, (slot_end-slot_start+3)*nprc_lcl, mpi_integer, 0, comm_ens, ierr )
1267 call mpi_bcast( bsna, (slot_end-slot_start+4)*nprc_lcl, mpi_integer, 0, comm_ens, ierr )
1269 call mpi_bcast( bsn, (slot_end-slot_start+3)*nprc_lcl, mpi_integer, 0, comm_lcl, ierr )
1270 call mpi_bcast( bsna, (slot_end-slot_start+4)*nprc_lcl, mpi_integer, 0, comm_lcl, ierr )
1272 do islot = slot_start, slot_end
1273 slot_id(islot) = islot - slot_start + 1
1274 slot_lb(islot) = (real(islot - slot_base,
rp) - 0.5d0) * slot_tinterval
1275 slot_ub(islot) = (real(islot - slot_base,
rp) + 0.5d0) * slot_tinterval
1282 nobs = bsna( slot_end+2, rank_lcl ) - bsna( slot_start-1, rank_lcl )
1284 obsda_tmp%nobs = nobs
1285 call obs_da_value_allocate(obsda_tmp, 0)
1287 if (
present(nobs_extern))
then
1288 obsda%nobs = nobs + nobs_extern
1292 call obs_da_value_allocate(obsda, nitmax)
1294 if (rank_ens == 0)
then
1295 allocate (cnts(nprc_lcl))
1296 allocate (dsps(nprc_lcl))
1297 do ip = 0, nprc_lcl-1
1298 dsps(ip+1) = bsna(slot_start-1, ip)
1299 cnts(ip+1) = bsna(slot_end+2, ip) - dsps(ip+1)
1302 call mpi_scatterv(obset_bufs, cnts, dsps, mpi_integer, obsda_tmp%set, cnts(rank_lcl+1), mpi_integer, 0, comm_lcl, ierr)
1303 call mpi_scatterv(obidx_bufs, cnts, dsps, mpi_integer, obsda_tmp%idx, cnts(rank_lcl+1), mpi_integer, 0, comm_lcl, ierr)
1305 deallocate (cnts, dsps)
1306 deallocate (obset_bufs, obidx_bufs)
1313 call mpi_bcast(obsda_tmp%set, nobs, mpi_integer, 0, comm_ens, ierr)
1314 call mpi_bcast(obsda_tmp%idx, nobs, mpi_integer, 0, comm_ens, ierr)
1316 obsda%set(1:nobs) = obsda_tmp%set
1317 obsda%idx(1:nobs) = obsda_tmp%idx
1323 allocate ( v3dg(nlevh,nlonh,nlath,nv3dd) )
1324 allocate ( v2dg( nlonh,nlath,nv2dd) )
1329 v3dg(k,i,j,iv3dd_u ) = u(k,i,j)
1330 v3dg(k,i,j,iv3dd_v ) = v(k,i,j)
1331 v3dg(k,i,j,iv3dd_w ) = w(k,i,j)
1332 v3dg(k,i,j,iv3dd_t ) = temp(k,i,j)
1333 v3dg(k,i,j,iv3dd_p ) = pres(k,i,j)
1334 v3dg(k,i,j,iv3dd_q ) = qv(k,i,j)
1335 v3dg(k,i,j,iv3dd_qc ) = qc(k,i,j)
1336 v3dg(k,i,j,iv3dd_qr ) = qr(k,i,j)
1337 v3dg(k,i,j,iv3dd_qi ) = qi(k,i,j)
1338 v3dg(k,i,j,iv3dd_qs ) = qs(k,i,j)
1339 v3dg(k,i,j,iv3dd_qg ) = qg(k,i,j)
1340 v3dg(k,i,j,iv3dd_rh ) = rh(k,i,j)
1341 v3dg(k,i,j,iv3dd_hgt) = hgt(k,i,j)
1347 v2dg(i,j,iv2dd_topo) = topo(i,j)
1348 v2dg(i,j,iv2dd_ps ) = ps(i,j)
1349 v2dg(i,j,iv2dd_rain) = rain(i,j)
1350 v2dg(i,j,iv2dd_u10m) = u10m(i,j)
1351 v2dg(i,j,iv2dd_v10m) = v10m(i,j)
1352 v2dg(i,j,iv2dd_t2m ) = t2m(i,j)
1353 v2dg(i,j,iv2dd_q2m ) = q2m(i,j)
1359 obsda_tmp%qc(1:nobs) = iqc_undef
1364 n1 = bsna(islot_time_out-1, rank_lcl) - bsna(slot_start-1, rank_lcl) + 1
1365 n2 = bsna(islot_time_out, rank_lcl) - bsna(slot_start-1, rank_lcl)
1367 obsda_tmp%qc(n1:n2) = iqc_time
1372 n1 = bsna(islot_domain_out-1, rank_lcl) - bsna(slot_start-1, rank_lcl) + 1
1373 n2 = bsna(islot_domain_out, rank_lcl) - bsna(slot_start-1, rank_lcl)
1375 obsda_tmp%qc(n1:n2) = iqc_out_h
1380 do islot = slot_start, slot_end
1381 n1 = bsna(islot-1, rank_lcl) - bsna(slot_start-1, rank_lcl) + 1
1382 n2 = bsna(islot, rank_lcl) - bsna(slot_start-1, rank_lcl)
1383 slot_nobsg = sum(bsn(islot, :))
1385 if (slot_nobsg <= 0)
then
1392 iof = obsda_tmp%set(nn)
1393 n = obsda_tmp%idx(nn)
1395 call rij_g2l(rank_lcl, obs(iof)%ri(n), obs(iof)%rj(n), ril, rjl)
1397 if (.not. use_obs(obs(iof)%typ(n)))
then
1398 obsda_tmp%qc(nn) = iqc_otype
1402 select case (obs_in_format(iof))
1406 call phys2ijk(v3dg(:,:,:,iv3dd_p), obs(iof)%elm(n), ril, rjl, obs(iof)%lev(n), rk, obsda_tmp%qc(nn))
1407 if (obsda_tmp%qc(nn) == iqc_good)
then
1408 call trans_xtoy(obs(iof)%elm(n), ril, rjl, rk, &
1409 obs(iof)%lon(n), obs(iof)%lat(n), v3dg, v2dg, obsda_tmp%val(nn), obsda_tmp%qc(nn))
1412 case (
'RADAR',
'PAWR_TOSHIBA',
'MP_PAWR_TOSHIBA',
'PAWR_JRC',
'HIMAWARI8' )
1414 if ( obs(iof)%lev(n) > radar_zmax .or. obs(iof)%lev(n) < radar_zmin )
then
1415 obsda_tmp%qc(nn) = iqc_radar_vhi
1417 call phys2ijkz(v3dg(:,:,:,iv3dd_hgt), ril, rjl, obs(iof)%lev(n), rkz, obsda_tmp%qc(nn))
1419 if (obsda_tmp%qc(nn) == iqc_good)
then
1420 call trans_xtoy_radar(obs(iof)%elm(n), obs(iof)%meta(1), obs(iof)%meta(2), obs(iof)%meta(3), ril, rjl, rkz, &
1421 obs(iof)%lon(n), obs(iof)%lat(n), obs(iof)%lev(n), v3dg, v2dg, obsda_tmp%val(nn), obsda_tmp%qc(nn))
1424 call itpl_3d( v3dg(:,:,:,iv3dd_p), rkz, ril, rjl, obsda_tmp%pm(nn) )
1425 call itpl_3d( v3dg(:,:,:,iv3dd_t), rkz, ril, rjl, obsda_tmp%tm(nn) )
1426 call itpl_3d( v3dg(:,:,:,iv3dd_q), rkz, ril, rjl, obsda_tmp%qv(nn) )
1434 where( obsda_tmp%qc(:) == iqc_ref_low )
1435 obsda_tmp%qc(:) = iqc_good
1440 call obs_da_value_partial_reduce_iter(obsda, it, 1, nobs, obsda_tmp%val, obsda_tmp%qc, &
1441 obsda_tmp%qv, obsda_tmp%tm, obsda_tmp%pm )
1444 deallocate ( v3dg, v2dg )
1445 deallocate ( bsn, bsna )
1447 call obs_da_value_deallocate( obsda_tmp )
1460 integer,
intent(in) :: obs_in_num
1461 integer,
optional,
intent(in) :: nobs_extern
1463 INTEGER :: n,i,j,ierr,im,iof,iidx
1468 real(
rp) :: qvs, qdry
1471 integer :: ityp,ielm,ielm_u,ictype
1472 real(
rp) :: target_grdspc
1474 integer :: myp_i,myp_j
1475 integer :: ip_i,ip_j
1480 integer :: nobs_elms_sum(
nid_obs)
1482 integer :: nobs_intern
1483 integer :: nobs_extern_
1486 character(len=3) :: use_obs_print
1487 character(4) :: nstr
1490 integer :: cntr(nprc_lcl)
1491 integer :: dspr(nprc_lcl)
1492 integer :: nensobs_div, nensobs_mod
1493 integer :: im_obs_1, im_obs_2, nensobs_part
1495 integer :: ns_ext, ne_ext, ns_bufr, ne_bufr
1496 integer :: ishift, jshift
1499 integer :: imin1,imax1,jmin1,jmax1,imin2,imax2,jmin2,jmax2
1501 real(
rp),
allocatable :: tmpelm(:)
1502 integer :: monit_nobs(
nid_obs)
1510 log_progress(*)
'data-assimilation / LETKF / obs / initialize'
1512 if(
present(nobs_extern) )
then
1513 nobs_extern_ = nobs_extern
1518 nobs_intern = obsda%nobs - nobs_extern_
1519 if( letkf_debug_log )
then
1520 log_info(
"LETKF_debug",
'(1x,A,I8)')
'Internally processed observations: ', nobs_intern
1521 log_info(
"LETKF_debug",
'(1x,A,I8)')
'Externally processed observations: ', nobs_extern_
1522 log_info(
"LETKF_debug",
'(1x,A,I8)')
'Total observations: ', obsda%nobs
1588 call obs_da_value_allreduce( obsda )
1597 ctype_use(:,:) = .false.
1598 do iof = 1, obs_in_num
1599 do n = 1, obs(iof)%nobs
1600 select case (obs(iof)%elm(n))
1601 case (id_radar_ref_obs)
1602 if (obs(iof)%dat(n) >= 0.0d0 .and. obs(iof)%dat(n) < 1.0d10)
then
1603 if (obs(iof)%dat(n) < min_radar_ref)
then
1604 obs(iof)%elm(n) = id_radar_ref_zero_obs
1605 obs(iof)%dat(n) = min_radar_ref_dbz + low_ref_shift
1607 obs(iof)%dat(n) = 10.0d0 * log10(obs(iof)%dat(n))
1610 obs(iof)%dat(n) = undef
1612 if (use_obserr_radar_ref)
then
1613 obs(iof)%err(n) = obserr_radar_ref
1615 case (id_radar_ref_zero_obs)
1616 obs(iof)%dat(n) = min_radar_ref_dbz + low_ref_shift
1617 if (use_obserr_radar_ref)
then
1618 obs(iof)%err(n) = obserr_radar_ref
1620 case (id_radar_vr_obs)
1621 if (use_obserr_radar_vr)
then
1622 obs(iof)%err(n) = obserr_radar_vr
1627 ctype_use(uid_obs(obs(iof)%elm(n)), obs(iof)%typ(n)) = .true.
1632 nctype = count(ctype_use)
1633 if (
allocated(elm_ctype ))
deallocate(elm_ctype )
1634 if (
allocated(elm_u_ctype ))
deallocate(elm_u_ctype )
1635 if (
allocated(typ_ctype ))
deallocate(typ_ctype )
1636 if (
allocated(hori_loc_ctype))
deallocate(hori_loc_ctype)
1637 if (
allocated(vert_loc_ctype))
deallocate(vert_loc_ctype)
1638 allocate (elm_ctype(nctype))
1639 allocate (elm_u_ctype(nctype))
1640 allocate (typ_ctype(nctype))
1641 allocate (hori_loc_ctype(nctype))
1642 allocate (vert_loc_ctype(nctype))
1644 ctype_elmtyp(:,:) = 0
1647 if (ctype_use(ielm_u, ityp))
then
1649 ctype_elmtyp(ielm_u, ityp) = ictype
1651 elm_ctype(ictype) = elem_uid(ielm_u)
1652 elm_u_ctype(ictype) = ielm_u
1653 typ_ctype(ictype) = ityp
1656 if (elm_ctype(ictype) == id_radar_ref_zero_obs)
then
1657 hori_loc_ctype(ictype) = hori_local_radar_obsnoref
1658 else if (elm_ctype(ictype) == id_radar_vr_obs)
then
1659 hori_loc_ctype(ictype) = hori_local_radar_vr
1661 hori_loc_ctype(ictype) = hori_local(ityp)
1664 if (elm_ctype(ictype) == id_radar_vr_obs)
then
1665 vert_loc_ctype(ictype) = vert_local_radar_vr
1667 vert_loc_ctype(ictype) = vert_local(ityp)
1679 allocate(tmpelm(obsda%nobs))
1681 do n = 1, obsda%nobs
1682 IF(obsda%qc(n) > 0) cycle
1687 tmpelm(n) = obs(iof)%elm(iidx)
1690 if (obs(iof)%elm(iidx) == id_radar_ref_obs .or. obs(iof)%elm(iidx) == id_radar_ref_zero_obs)
then
1691 if (.not. use_radar_ref)
then
1692 obsda%qc(n) = iqc_otype
1696 if (obs(iof)%dat(iidx) == undef)
then
1697 obsda%qc(n) = iqc_obs_bad
1704 if (obsda%ensval(i,n) > radar_ref_thres_dbz+1.0d-6 )
then
1705 mem_ref = mem_ref + 1
1710 if (obs(iof)%dat(iidx) > radar_ref_thres_dbz+1.0d-6)
then
1711 if (mem_ref < min_radar_ref_member_obsrain)
then
1712 obsda%qc(n) = iqc_ref_mem
1714 if ( .not. radar_pqv ) cycle
1721 if (mem_ref < min_radar_ref_member_obsnorain)
then
1722 obsda%qc(n) = iqc_ref_mem
1729 if (obs(iof)%elm(iidx) == id_radar_vr_obs)
then
1730 if (.not. use_radar_vr)
then
1731 obsda%qc(n) = iqc_otype
1737 obsda%val(n) = 0.0_rp
1739 obsda%val(n) = obsda%val(n) + obsda%ensval(i,n)
1741 obsda%val(n) = obsda%val(n) / real(nmem,kind=
rp)
1744 obsda%ensval(i,n) = obsda%ensval(i,n) - obsda%val(n)
1746 obsda%val(n) = obs(iof)%dat(iidx) - obsda%val(n)
1747 if( ens_with_mdet )
then
1748 obsda%ensval(mmdetobs,n) = obs(iof)%dat(iidx) - obsda%ensval(mmdetobs,n)
1751 select case (obs(iof)%elm(iidx))
1753 IF(abs(obsda%val(n)) > gross_error_rain * obs(iof)%err(iidx))
THEN
1754 obsda%qc(n) = iqc_gross_err
1756 case (id_radar_ref_obs,id_radar_ref_zero_obs)
1758 if( radar_pqv .and. obsda%val(n) > radar_pqv_omb )
then
1761 obsda%val(n) = 0.0_rp
1763 obsda%val(n) = obsda%val(n) + obsda%eqv(i,n)
1765 obsda%val(n) = obsda%val(n) / real(nmem, kind=
rp)
1768 obsda%ensval(i,n) = obsda%eqv(i,n) - obsda%val(n)
1772 qvs = 611.2d0*exp(17.67d0*(obsda%tm(n)-tem00)/(obsda%tm(n) - tem00 + 243.5d0))
1775 qvs = 0.622d0*qvs / ( obsda%pm(n) - qvs )
1777 obsda%val(n) = qvs - obsda%val(n)
1779 if (ens_with_mdet)
then
1780 obsda%ensval(mmdetobs,n) = qvs - obsda%eqv(mmdetobs,n)
1783 obsda%tm(n) = -1.0d0
1786 IF(abs(obsda%val(n)) > gross_error_radar_ref * obs(iof)%err(iidx))
THEN
1787 obsda%qc(n) = iqc_gross_err
1790 case (id_radar_vr_obs)
1791 IF(abs(obsda%val(n)) > gross_error_radar_vr * obs(iof)%err(iidx))
THEN
1792 obsda%qc(n) = iqc_gross_err
1794 case (id_radar_prh_obs)
1795 IF(abs(obsda%val(n)) > gross_error_radar_prh * obs(iof)%err(iidx))
THEN
1796 obsda%qc(n) = iqc_gross_err
1799 IF(abs(obsda%val(n)) > gross_error * obs(iof)%err(iidx))
THEN
1800 obsda%qc(n) = iqc_gross_err
1806 if( letkf_debug_log )
then
1807 log_info(
"LETKF_debug",
'(1x,A,I6,A)')
'OBSERVATIONAL DEPARTURE STATISTICS (IN THIS SUBDOMAIN #', rank_lcl,
'):'
1809 call monit_dep( obsda%nobs, tmpelm, obsda%val, obsda%qc, monit_nobs, bias, rmse )
1810 call monit_print( monit_nobs, bias, rmse )
1819 if (
allocated(obsgrd))
deallocate(obsgrd)
1820 allocate (obsgrd(nctype))
1825 do ictype = 1, nctype
1826 ityp = typ_ctype(ictype)
1828 if( obs_sort_grid_spacing(ityp) > 0.0_rp )
then
1829 target_grdspc = obs_sort_grid_spacing(ityp)
1830 else if( max_nobs_per_grid(ityp) > 0 )
then
1831 target_grdspc = 0.1_rp * sqrt( real(max_nobs_per_grid(ityp), kind=
rp) ) * obs_min_spacing(ityp)
1833 target_grdspc = hori_loc_ctype(ictype) * dist_zero_fac / 6.0_rp
1835 obsgrd(ictype)%ngrd_i = min( ceiling( dx * real( nlon, kind=
rp ) / target_grdspc), nlon )
1836 obsgrd(ictype)%ngrd_j = min( ceiling( dy * real( nlat, kind=
rp ) / target_grdspc), nlat )
1837 obsgrd(ictype)%grdspc_i = dx * real( nlon, kind=
rp ) / real( obsgrd( ictype )%ngrd_i, kind=
rp )
1838 obsgrd(ictype)%grdspc_j = dy * real( nlat, kind=
rp ) / real( obsgrd( ictype )%ngrd_j, kind=
rp )
1839 if( letkf_entire_grid_search_x )
then
1840 obsgrd(ictype)%ngrdsch_i = ceiling( dx * nlong / obsgrd( ictype )%grdspc_i )
1842 obsgrd(ictype)%ngrdsch_i = ceiling( hori_loc_ctype( ictype ) * dist_zero_fac / obsgrd( ictype )%grdspc_i )
1844 if( letkf_entire_grid_search_y )
then
1845 obsgrd(ictype)%ngrdsch_j = ceiling( dy * nlatg / obsgrd( ictype )%grdspc_j )
1847 obsgrd(ictype)%ngrdsch_j = ceiling( hori_loc_ctype( ictype ) * dist_zero_fac / obsgrd( ictype )%grdspc_j )
1849 obsgrd(ictype)%ngrdext_i = obsgrd( ictype )%ngrd_i + obsgrd( ictype )%ngrdsch_i * 2
1850 obsgrd(ictype)%ngrdext_j = obsgrd( ictype )%ngrd_j + obsgrd( ictype )%ngrdsch_j * 2
1852 allocate (obsgrd(ictype)%n ( obsgrd(ictype)%ngrd_i, obsgrd(ictype)%ngrd_j, 0:nprc_lcl-1))
1853 allocate (obsgrd(ictype)%ac(0:obsgrd(ictype)%ngrd_i, obsgrd(ictype)%ngrd_j, 0:nprc_lcl-1))
1854 allocate (obsgrd(ictype)%tot(0:nprc_lcl-1))
1855 allocate (obsgrd(ictype)%n_ext ( obsgrd(ictype)%ngrdext_i, obsgrd(ictype)%ngrdext_j))
1856 allocate (obsgrd(ictype)%ac_ext(0:obsgrd(ictype)%ngrdext_i, obsgrd(ictype)%ngrdext_j))
1858 obsgrd(ictype)%n (:,:,:) = 0
1859 obsgrd(ictype)%ac(:,:,:) = 0
1860 obsgrd(ictype)%tot(:) = 0
1861 obsgrd(ictype)%n_ext (:,:) = 0
1862 obsgrd(ictype)%ac_ext(:,:) = 0
1863 obsgrd(ictype)%tot_ext = 0
1864 obsgrd(ictype)%tot_sub(:) = 0
1865 obsgrd(ictype)%tot_g(:) = 0
1867 allocate (obsgrd(ictype)%next(obsgrd(ictype)%ngrd_i, obsgrd(ictype)%ngrd_j))
1873 do n = 1, obsda%nobs
1876 ictype = ctype_elmtyp(uid_obs(obs(iof)%elm(iidx)), obs(iof)%typ(iidx))
1878 if (obsda%qc(n) == iqc_good)
then
1879 call ij_obsgrd( ictype, obs(iof)%ri(iidx), obs(iof)%rj(iidx), i, j )
1881 if (i > obsgrd(ictype)%ngrd_i) i = obsgrd(ictype)%ngrd_i
1883 if (j > obsgrd(ictype)%ngrd_j) j = obsgrd(ictype)%ngrd_j
1885 obsgrd(ictype)%n(i,j,rank_lcl) = obsgrd(ictype)%n(i,j,rank_lcl) + 1
1895 do ictype = 1, nctype
1896 if (ictype > 1)
then
1897 obsgrd(ictype)%ac(0,1,rank_lcl) = obsgrd(ictype-1)%ac(obsgrd(ictype-1)%ngrd_i,obsgrd(ictype-1)%ngrd_j,rank_lcl)
1899 do j = 1, obsgrd(ictype)%ngrd_j
1901 obsgrd(ictype)%ac(0,j,rank_lcl) = obsgrd(ictype)%ac(obsgrd(ictype)%ngrd_i,j-1,rank_lcl)
1903 do i = 1, obsgrd(ictype)%ngrd_i
1904 obsgrd(ictype)%ac(i,j,rank_lcl) = obsgrd(ictype)%ac(i-1,j,rank_lcl) + obsgrd(ictype)%n(i,j,rank_lcl)
1907 obsgrd(ictype)%next(1:obsgrd(ictype)%ngrd_i,:) = obsgrd(ictype)%ac(0:obsgrd(ictype)%ngrd_i-1,:,rank_lcl)
1913 do n = 1, obsda%nobs
1914 if (obsda%qc(n) == iqc_good)
then
1917 ictype = ctype_elmtyp(uid_obs(obs(iof)%elm(iidx)), obs(iof)%typ(iidx))
1919 call ij_obsgrd( ictype, obs(iof)%ri(iidx), obs(iof)%rj(iidx), i, j )
1921 if (i > obsgrd(ictype)%ngrd_i) i = obsgrd(ictype)%ngrd_i
1923 if (j > obsgrd(ictype)%ngrd_j) j = obsgrd(ictype)%ngrd_j
1925 obsgrd(ictype)%next(i,j) = obsgrd(ictype)%next(i,j) + 1
1926 obsda%key(obsgrd(ictype)%next(i,j)) = n
1935 do ictype = 1, nctype
1936 if (nprc_lcl > 1)
then
1937 call mpi_allreduce(mpi_in_place, obsgrd(ictype)%n, obsgrd(ictype)%ngrd_i*obsgrd(ictype)%ngrd_j*nprc_lcl, &
1938 mpi_integer, mpi_sum, comm_lcl, ierr)
1939 call mpi_allreduce(mpi_in_place, obsgrd(ictype)%ac(0:obsgrd(ictype)%ngrd_i,:,:), (obsgrd(ictype)%ngrd_i+1)*obsgrd(ictype)%ngrd_j*nprc_lcl, &
1940 mpi_integer, mpi_sum, comm_lcl, ierr)
1942 call mpi_allreduce(obsgrd(ictype)%tot_sub, obsgrd(ictype)%tot_g,
n_qc_steps, mpi_integer, mpi_sum, comm_lcl, ierr)
1944 if (ictype == 1)
then
1945 obsgrd(ictype)%tot(:) = obsgrd(ictype)%ac(obsgrd(ictype)%ngrd_i,obsgrd(ictype)%ngrd_j,:)
1947 obsgrd(ictype)%tot(:) = obsgrd(ictype)%ac(obsgrd(ictype)%ngrd_i,obsgrd(ictype)%ngrd_j,:) &
1948 - obsgrd(ictype-1)%ac(obsgrd(ictype-1)%ngrd_i,obsgrd(ictype-1)%ngrd_j,:)
1951 nobs_sub(:) = nobs_sub(:) + obsgrd(ictype)%tot_sub(:)
1952 nobs_g(:) = nobs_g(:) + obsgrd(ictype)%tot_g(:)
1954 deallocate (obsgrd(ictype)%next)
1963 call rank_1d_2d(rank_lcl, myp_i, myp_j)
1965 do ictype = 1, nctype
1966 imin1 = myp_i*obsgrd(ictype)%ngrd_i+1 - obsgrd(ictype)%ngrdsch_i
1967 imax1 = (myp_i+1)*obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
1968 jmin1 = myp_j*obsgrd(ictype)%ngrd_j+1 - obsgrd(ictype)%ngrdsch_j
1969 jmax1 = (myp_j+1)*obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
1971 do ip = 0, nprc_lcl-1
1972 call rank_1d_2d(ip, ip_i, ip_j)
1973 imin2 = max(1, imin1 - ip_i*obsgrd(ictype)%ngrd_i)
1974 imax2 = min(obsgrd(ictype)%ngrd_i, imax1 - ip_i*obsgrd(ictype)%ngrd_i)
1975 jmin2 = max(1, jmin1 - ip_j*obsgrd(ictype)%ngrd_j)
1976 jmax2 = min(obsgrd(ictype)%ngrd_j, jmax1 - ip_j*obsgrd(ictype)%ngrd_j)
1977 if (imin2 > imax2 .or. jmin2 > jmax2) cycle
1979 ishift = (ip_i - myp_i) * obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
1980 jshift = (ip_j - myp_j) * obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
1981 obsgrd(ictype)%n_ext(imin2+ishift:imax2+ishift, jmin2+jshift:jmax2+jshift) = obsgrd(ictype)%n(imin2:imax2, jmin2:jmax2, ip)
1984 if (ictype > 1)
then
1985 obsgrd(ictype)%ac_ext(0,1) = obsgrd(ictype-1)%ac_ext(obsgrd(ictype-1)%ngrdext_i,obsgrd(ictype-1)%ngrdext_j)
1987 do j = 1, obsgrd(ictype)%ngrdext_j
1989 obsgrd(ictype)%ac_ext(0,j) = obsgrd(ictype)%ac_ext(obsgrd(ictype)%ngrdext_i,j-1)
1991 do i = 1, obsgrd(ictype)%ngrdext_i
1992 obsgrd(ictype)%ac_ext(i,j) = obsgrd(ictype)%ac_ext(i-1,j) + obsgrd(ictype)%n_ext(i,j)
1996 if (ictype == 1)
then
1997 obsgrd(ictype)%tot_ext = obsgrd(ictype)%ac_ext(obsgrd(ictype)%ngrdext_i,obsgrd(ictype)%ngrdext_j)
1999 obsgrd(ictype)%tot_ext = obsgrd(ictype)%ac_ext(obsgrd(ictype)%ngrdext_i,obsgrd(ictype)%ngrdext_j) &
2000 - obsgrd(ictype-1)%ac_ext(obsgrd(ictype-1)%ngrdext_i,obsgrd(ictype-1)%ngrdext_j)
2004 if (nctype > 0)
then
2005 nobstotal = obsgrd(nctype)%ac_ext(obsgrd(nctype)%ngrdext_i,obsgrd(nctype)%ngrdext_j)
2007 maxnobs_per_ctype = obsgrd(1)%tot_ext
2008 do ictype = 2, nctype
2009 maxnobs_per_ctype = max(maxnobs_per_ctype, obsgrd(ictype)%tot_ext)
2013 maxnobs_per_ctype = 0
2022 if (nctype > 0)
then
2023 cntr(:) = obsgrd(nctype)%ac(obsgrd(nctype)%ngrd_i,obsgrd(nctype)%ngrd_j,:)
2024 cnts = cntr(rank_lcl+1)
2027 dspr(ip) = dspr(ip-1) + cntr(ip-1)
2030 nensobs_mod = mod(nensobs, nprc_ens)
2031 nensobs_div = (nensobs - nensobs_mod) / nprc_ens
2032 if (rank_ens < nensobs_mod)
then
2033 im_obs_1 = (nensobs_div+1) * rank_ens + 1
2034 im_obs_2 = (nensobs_div+1) * (rank_ens+1)
2035 nensobs_part = nensobs_div + 1
2037 im_obs_1 = (nensobs_div+1) * nensobs_mod + nensobs_div * (rank_ens-nensobs_mod) + 1
2038 im_obs_2 = (nensobs_div+1) * nensobs_mod + nensobs_div * (rank_ens-nensobs_mod+1)
2039 nensobs_part = nensobs_div
2043 call obs_da_value_allocate(obsbufs, nensobs_part)
2046 obsbufs%set(n) = obsda%set(obsda%key(n))
2047 obsbufs%idx(n) = obsda%idx(obsda%key(n))
2048 obsbufs%val(n) = obsda%val(obsda%key(n))
2049 obsbufs%tm(n) = obsda%tm(obsda%key(n))
2050 if (nensobs_part > 0)
then
2051 obsbufs%ensval(1:nensobs_part,n) = obsda%ensval(im_obs_1:im_obs_2,obsda%key(n))
2053 obsbufs%qc(n) = obsda%qc(obsda%key(n))
2057 call obs_da_value_deallocate(obsda)
2063 if (nctype > 0)
then
2065 call obs_da_value_allocate(obsbufr, nensobs_part)
2067 call mpi_allgatherv( obsbufs%set, cnts, mpi_integer, obsbufr%set, cntr, dspr, mpi_integer, comm_lcl, ierr )
2068 call mpi_allgatherv( obsbufs%idx, cnts, mpi_integer, obsbufr%idx, cntr, dspr, mpi_integer, comm_lcl, ierr )
2069 call mpi_allgatherv( obsbufs%val, cnts, datatype, obsbufr%val, cntr, dspr, datatype, comm_lcl, ierr )
2070 call mpi_allgatherv( obsbufs%tm, cnts, datatype, obsbufr%tm, cntr, dspr, datatype, comm_lcl, ierr )
2071 if (nensobs_part > 0)
then
2072 call mpi_allgatherv(obsbufs%ensval, cnts*nensobs_part, datatype, obsbufr%ensval, cntr*nensobs_part, dspr*nensobs_part, datatype, comm_lcl, ierr)
2074 call mpi_allgatherv(obsbufs%qc, cnts, mpi_integer, obsbufr%qc, cntr, dspr, mpi_integer, comm_lcl, ierr)
2076 call obs_da_value_deallocate(obsbufs)
2082 obsda_sort%nobs = nobstotal
2083 call obs_da_value_allocate(obsda_sort, nensobs)
2085 do ip = 0, nprc_lcl-1
2086 call rank_1d_2d(ip, ip_i, ip_j)
2088 do ictype = 1, nctype
2089 imin1 = myp_i*obsgrd(ictype)%ngrd_i+1 - obsgrd(ictype)%ngrdsch_i
2090 imax1 = (myp_i+1)*obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
2091 jmin1 = myp_j*obsgrd(ictype)%ngrd_j+1 - obsgrd(ictype)%ngrdsch_j
2092 jmax1 = (myp_j+1)*obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
2094 imin2 = max(1, imin1 - ip_i*obsgrd(ictype)%ngrd_i)
2095 imax2 = min(obsgrd(ictype)%ngrd_i, imax1 - ip_i*obsgrd(ictype)%ngrd_i)
2096 jmin2 = max(1, jmin1 - ip_j*obsgrd(ictype)%ngrd_j)
2097 jmax2 = min(obsgrd(ictype)%ngrd_j, jmax1 - ip_j*obsgrd(ictype)%ngrd_j)
2098 if (imin2 > imax2 .or. jmin2 > jmax2) cycle
2100 ishift = (ip_i - myp_i) * obsgrd(ictype)%ngrd_i + obsgrd(ictype)%ngrdsch_i
2101 jshift = (ip_j - myp_j) * obsgrd(ictype)%ngrd_j + obsgrd(ictype)%ngrdsch_j
2104 ns_ext = obsgrd(ictype)%ac_ext(imin2+ishift-1,j+jshift) + 1
2105 ne_ext = obsgrd(ictype)%ac_ext(imax2+ishift ,j+jshift)
2106 if (ns_ext > ne_ext) cycle
2108 ns_bufr = dspr(ip+1) + obsgrd(ictype)%ac(imin2-1,j,ip) + 1
2109 ne_bufr = dspr(ip+1) + obsgrd(ictype)%ac(imax2 ,j,ip)
2111 obsda_sort%set(ns_ext:ne_ext) = obsbufr%set(ns_bufr:ne_bufr)
2112 obsda_sort%idx(ns_ext:ne_ext) = obsbufr%idx(ns_bufr:ne_bufr)
2113 obsda_sort%val(ns_ext:ne_ext) = obsbufr%val(ns_bufr:ne_bufr)
2114 obsda_sort%tm(ns_ext:ne_ext) = obsbufr%tm(ns_bufr:ne_bufr)
2115 if (nensobs_part > 0)
then
2116 obsda_sort%ensval(im_obs_1:im_obs_2,ns_ext:ne_ext) = obsbufr%ensval(1:nensobs_part,ns_bufr:ne_bufr)
2118 obsda_sort%qc(ns_ext:ne_ext) = obsbufr%qc(ns_bufr:ne_bufr)
2124 obsda_sort%nobs_in_key = 0
2125 do ictype = 1, nctype
2126 imin1 = obsgrd(ictype)%ngrdsch_i + 1
2127 imax1 = obsgrd(ictype)%ngrdsch_i + obsgrd(ictype)%ngrd_i
2128 jmin1 = obsgrd(ictype)%ngrdsch_j + 1
2129 jmax1 = obsgrd(ictype)%ngrdsch_j + obsgrd(ictype)%ngrd_j
2130 call obs_choose_ext(ictype, imin1, imax1, jmin1, jmax1, obsda_sort%nobs_in_key, obsda_sort%key)
2132 deallocate (obsgrd(ictype)%n)
2133 deallocate (obsgrd(ictype)%ac)
2134 deallocate (obsgrd(ictype)%n_ext)
2137 if (nctype > 0)
then
2138 call obs_da_value_deallocate(obsbufr)
2144 if (nprc_ens > 1)
then
2145 call mpi_allreduce(mpi_in_place, obsda_sort%ensval, nensobs*nobstotal, datatype, mpi_sum, comm_ens, ierr)
2148 if( letkf_debug_log )
then
2149 log_info(
"LETKF_debug",
'(1x,A,I6,A)')
'OBSERVATION COUNTS (GLOABL AND IN THIS SUBDOMAIN #', rank_lcl,
'):'
2150 log_info(
"LETKF_debug",
'(1x,A)')
'====================================================================='
2151 log_info(
"LETKF_debug",
'(1x,A)')
'TYPE VAR GLOBAL GLOBAL SUBDOMAIN SUBDOMAIN EXT_SUBDOMAIN'
2152 log_info(
"LETKF_debug",
'(1x,A)')
' before QC after QC before QC after QC after QC'
2153 log_info(
"LETKF_debug",
'(1x,A)')
'---------------------------------------------------------------------'
2154 do ictype = 1, nctype
2155 ityp = typ_ctype(ictype)
2156 ielm_u = elm_u_ctype(ictype)
2157 log_info(
"LETKF_debug",
'(1x,A6,1x,A3,1x,4I11,I14)') obtypelist(ityp), obelmlist(ielm_u), &
2162 obsgrd(ictype)%tot_ext
2164 log_info(
"LETKF_debug",
'(1x,A)')
'---------------------------------------------------------------------'
2166 log_info(
"LETKF_debug",
'(1x,A)')
'====================================================================='
2176 OBS_IN_NUM, & ! [IN]
2177 OBS_IN_FORMAT, & ! [IN]
2194 average => statistics_average
2201 integer,
intent(in) :: obs_in_num
2202 character(len=H_LONG),
intent(in) :: obs_in_format(:)
2204 real(
rp),
intent(inout) :: u (nlev,nlon,nlat)
2205 real(
rp),
intent(inout) :: v (nlev,nlon,nlat)
2206 real(
rp),
intent(inout) :: w (nlev,nlon,nlat)
2207 real(
rp),
intent(inout) :: temp(nlev,nlon,nlat)
2208 real(
rp),
intent(inout) :: pres(nlev,nlon,nlat)
2209 real(
rp),
intent(inout) :: qv (nlev,nlon,nlat)
2210 real(
rp),
intent(inout) :: qc (nlev,nlon,nlat)
2211 real(
rp),
intent(inout) :: qr (nlev,nlon,nlat)
2212 real(
rp),
intent(inout) :: qi (nlev,nlon,nlat)
2213 real(
rp),
intent(inout) :: qs (nlev,nlon,nlat)
2214 real(
rp),
intent(inout) :: qg (nlev,nlon,nlat)
2216 real(
rp) :: gues3d(nij1,nlev,nens+1,
nv3d)
2217 real(
rp) :: gues2d(nij1, nens+1,
nv2d)
2219 real(
rp) :: anal3d(nij1,nlev,nens+1,
nv3d)
2220 real(
rp) :: anal2d(nij1, nens+1,
nv2d)
2222 real(
rp) :: addi3d(nij1,nlev,nens+1,
nv3d)
2223 real(
rp) :: addi2d(nij1, nens+1,
nv2d)
2225 real(
rp) :: v3dg(nlev,nlon,nlat,
nv3d)
2226 real(
rp) :: v2dg( nlon,nlat,
nv2d)
2228 real(
rp) :: work3d(nij1,nlev,
nv3d)
2229 real(
rp) :: work2d(nij1,
nv2d)
2230 real(
rp),
allocatable :: work3da(:,:,:)
2231 real(
rp),
allocatable :: work2da(:,:)
2232 real(
rp),
allocatable :: work3dn(:,:,:,:)
2233 real(
rp),
allocatable :: work2dn(:,:,:)
2234 real(
rp),
allocatable :: work3dg(:,:,:,:)
2235 real(
rp),
allocatable :: work2dg(:,:,:)
2237 real(
rp),
allocatable :: hdxf(:,:)
2238 real(
rp),
allocatable :: rdiag(:)
2239 real(
rp),
allocatable :: rloc(:)
2240 real(
rp),
allocatable :: dep(:)
2241 real(
rp),
allocatable :: depd(:)
2245 integer :: var_local_n2nc_max
2246 integer :: var_local_n2nc(
nv3d+
nv2d)
2247 integer :: var_local_n2n(
nv3d+
nv2d)
2248 logical :: var_local_n2n_found
2249 integer :: n2n, n2nc
2252 real(
rp),
allocatable :: trans(:,:,:)
2253 real(
rp),
allocatable :: transm(:,:)
2254 real(
rp),
allocatable :: transmd(:,:)
2255 real(
rp),
allocatable :: pa(:,:,:)
2256 real(
rp) :: transrlx(nmem,nmem)
2259 integer :: ij,ilev,n,m,i,j,k,nobsl
2264 real(
rp) :: q_mean,q_sprd
2265 real(
rp) :: q_anal(nmem)
2267 integer :: mshuf,ierr
2268 integer :: ishuf(nmem)
2269 real(
rp),
allocatable :: addinfl_weight(:)
2270 real(
rp) :: rdx,rdy,rdxy,ref_min_dist
2271 integer :: ic,ic2,iob
2273 integer,
allocatable :: search_q0(:,:,:)
2275 log_progress(*)
'data-assimilation / LETKF / system'
2281 v3dg(k,i,j,iv3d_u ) = u(k,i,j)
2282 v3dg(k,i,j,iv3d_v ) = v(k,i,j)
2283 v3dg(k,i,j,iv3d_w ) = w(k,i,j)
2284 v3dg(k,i,j,iv3d_t ) = temp(k,i,j)
2285 v3dg(k,i,j,iv3d_p ) = pres(k,i,j)
2286 v3dg(k,i,j,iv3d_q ) = qv(k,i,j)
2287 v3dg(k,i,j,iv3d_qc) = qc(k,i,j)
2288 v3dg(k,i,j,iv3d_qr) = qr(k,i,j)
2289 v3dg(k,i,j,iv3d_qi) = qi(k,i,j)
2290 v3dg(k,i,j,iv3d_qs) = qs(k,i,j)
2291 v3dg(k,i,j,iv3d_qg) = qg(k,i,j)
2301 if( letkf_debug_log )
then
2302 call monit_obs_mpi( obs_in_num, obs_in_format, v3dg, v2dg, monit_step=1 )
2305 call scatter_grd_mpi_all2all( 1, nens, v3dg, v2dg, gues3d(:,:,1:nens,:), gues2d(:,1:nens,:) )
2311 gues3d(ij,k,mmean,n) = average( nmem, gues3d(ij,k,1:nmem,n), undef )
2317 gues2d(ij,mmean,n) = average( nmem, gues2d(ij,1:nmem,n), undef )
2324 var_local(:,1) = var_local_uv(:)
2325 var_local(:,2) = var_local_t(:)
2326 var_local(:,3) = var_local_q(:)
2327 var_local(:,4) = var_local_ps(:)
2328 var_local(:,5) = var_local_rain(:)
2329 var_local(:,6) = var_local_tc(:)
2330 var_local(:,7) = var_local_radar_ref(:)
2331 var_local(:,8) = var_local_radar_vr(:)
2333 var_local_n2nc_max = 1
2334 var_local_n2nc(1) = 1
2335 var_local_n2n(1) = 1
2338 var_local_n2n_found = .false.
2339 do i = 1, var_local_n2nc_max
2341 if (all(var_local(var_local_n2nc(i),:) == var_local(n,:)))
then
2342 var_local_n2nc(n) = var_local_n2nc(i)
2343 var_local_n2n(n) = var_local_n2n(var_local_n2nc(n))
2344 var_local_n2n_found = .true.
2348 if( .NOT. var_local_n2n_found )
then
2349 var_local_n2nc_max = var_local_n2nc_max + 1
2350 var_local_n2nc(n) = var_local_n2nc_max
2351 var_local_n2n(n) = n
2358 ctype_merge(:,:) = 0
2359 ctype_merge(uid_obs(id_radar_ref_obs),22) = 1
2360 ctype_merge(uid_obs(id_radar_ref_zero_obs),22) = 1
2362 allocate (n_merge(nctype))
2366 if (n_merge(ic) > 0)
then
2368 if (ctype_merge(elm_u_ctype(ic),typ_ctype(ic)) > 0)
then
2369 do ic2 = ic+1, nctype
2370 if (ctype_merge(elm_u_ctype(ic2),typ_ctype(ic2)) == ctype_merge(elm_u_ctype(ic),typ_ctype(ic)))
then
2371 n_merge(ic) = n_merge(ic) + 1
2372 ic_merge(n_merge(ic),ic) = ic2
2379 n_merge_max = maxval(n_merge)
2381 allocate (search_q0(nctype,
nv3d+1,nij1))
2382 search_q0(:,:,:) = 1
2386 if (obtypelist(typ_ctype(ic)) /=
'PHARAD')
then
2387 radar_only = .false.
2400 gues3d(i,k,m,n) = gues3d(i,k,m,n) - gues3d(i,k,mmean,n)
2408 gues2d(i,m,n) = gues2d(i,m,n) - gues2d(i,mmean,n)
2416 IF(infl_mul > 0.0d0)
THEN
2420 log_error(
"LETKF_system",*)
'This system has not been implemented yet. INFL_MUL must be greather than 0.0.'
2432 IF(infl_mul_min > 0.0d0)
THEN
2433 work3d = max(work3d, infl_mul_min)
2434 work2d = max(work2d, infl_mul_min)
2438 allocate( hdxf( nobstotal, nmem ) )
2439 allocate (rdiag(nobstotal))
2440 allocate (rloc(nobstotal))
2441 allocate (dep(nobstotal))
2442 if( ens_with_mdet )
then
2443 allocate (depd(nobstotal))
2445 allocate (trans(nmem,nmem,var_local_n2nc_max))
2446 allocate (transm(nmem, var_local_n2nc_max))
2447 allocate (transmd(nmem, var_local_n2nc_max))
2448 allocate (pa(nmem,nmem,var_local_n2nc_max))
2456 trans_done(:) = .false.
2460 call relax_beta(rig1(ij),rjg1(ij),hgt1(ij,ilev),beta)
2462 if (beta == 0.0d0)
then
2465 anal3d(ij,ilev,m,n) = gues3d(ij,ilev,mmean,n) + gues3d(ij,ilev,m,n)
2467 if (ens_with_mdet)
then
2468 anal3d(ij,ilev,mmdet,n) = gues3d(ij,ilev,mmdet,n)
2474 anal2d(ij,m,n) = gues2d(ij,mmean,n) + gues2d(ij,m,n)
2476 if (ens_with_mdet)
then
2477 anal2d(ij,mmdet,n) = gues2d(ij,mmdet,n)
2488 n2nc = var_local_n2nc(n)
2489 n2n = var_local_n2n(n)
2491 if (gues3d(ij,ilev,mmean,iv3d_p) < q_update_top .and. n >= iv3d_q .and. n <= iv3d_qg)
then
2493 anal3d(ij,ilev,m,n) = gues3d(ij,ilev,mmean,n) + gues3d(ij,ilev,m,n)
2495 if (ens_with_mdet)
then
2496 anal3d(ij,ilev,mmdet,n) = gues3d(ij,ilev,mmdet,n)
2502 if (relax_to_inflated_prior)
then
2503 parm = work3d(ij,ilev,n)
2509 if (trans_done(n2nc))
then
2511 if (infl_mul_adaptive)
then
2512 work3d(ij,ilev,n) = work3d(ij,ilev,n2n)
2515 work3dn(:,ij,ilev,n) = work3dn(:,ij,ilev,n2n)
2520 if (ens_with_mdet)
then
2521 CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),n,&
2522 hdxf,rdiag,rloc,dep,nobsl,depd=depd,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,n,ij))
2524 CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),n,&
2525 hdxf,rdiag,rloc,dep,nobsl,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,n,ij))
2527 IF(relax_alpha_spread /= 0.0d0)
THEN
2528 if (ens_with_mdet)
then
2529 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), &
2530 trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), &
2531 rdiag_wloc=.true.,infl_update=infl_mul_adaptive, &
2532 depd=depd,transmd=transmd(:,n2nc))
2534 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), &
2535 trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), &
2536 rdiag_wloc=.true.,infl_update=infl_mul_adaptive)
2539 if (ens_with_mdet)
then
2540 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), &
2541 trans(:,:,n2nc),transm=transm(:,n2nc), &
2542 rdiag_wloc=.true.,infl_update=infl_mul_adaptive, &
2543 depd=depd,transmd=transmd(:,n2nc))
2545 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work3d(ij,ilev,n), &
2546 trans(:,:,n2nc),transm=transm(:,n2nc), &
2547 rdiag_wloc=.true.,infl_update=infl_mul_adaptive)
2550 trans_done(n2nc) = .true.
2552 work3dn(:,ij,ilev,n) = real(sum(nobsl_t, dim=1),
rp)
2553 work3dn(
nobtype+1,ij,ilev,n) = real(nobsl_t(9,22),
rp)
2554 work3dn(
nobtype+2,ij,ilev,n) = real(nobsl_t(10,22),
rp)
2555 work3dn(
nobtype+3,ij,ilev,n) = real(nobsl_t(11,22),
rp)
2556 work3dn(
nobtype+4,ij,ilev,n) = real(cutd_t(9,22),
rp)
2557 work3dn(
nobtype+5,ij,ilev,n) = real(cutd_t(10,22),
rp)
2558 work3dn(
nobtype+6,ij,ilev,n) = real(cutd_t(11,22),
rp)
2564 IF(relax_alpha /= 0.0d0)
THEN
2565 CALL weight_rtpp(trans(:,:,n2nc),parm,transrlx)
2566 ELSE IF(relax_alpha_spread /= 0.0d0)
THEN
2567 IF(relax_spread_out)
THEN
2568 CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues3d(ij,ilev,:,n), &
2569 parm,transrlx,work3da(ij,ilev,n))
2571 CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues3d(ij,ilev,:,n), &
2572 parm,transrlx,tmpinfl)
2575 transrlx = trans(:,:,n2nc)
2581 transrlx(k,m) = (transrlx(k,m) + transm(k,n2nc)) * beta
2583 transrlx(m,m) = transrlx(m,m) + (1.0_rp-beta)
2587 anal3d(ij,ilev,m,n) = gues3d(ij,ilev,mmean,n)
2589 anal3d(ij,ilev,m,n) = anal3d(ij,ilev,m,n) + gues3d(ij,ilev,k,n) * transrlx(k,m)
2593 if( ens_with_mdet )
then
2595 anal3d(ij,ilev,mmdet,n) = 0.0_rp
2597 anal3d(ij,ilev,mmdet,n) = anal3d(ij,ilev,mmdet,n) + gues3d(ij,ilev,k,n) * transmd(k,n2nc)
2599 anal3d(ij,ilev,mmdet,n) = gues3d(ij,ilev,mmdet,n) + anal3d(ij,ilev,mmdet,n) * beta
2603 IF(q_sprd_max > 0.0d0 .and. n == iv3d_q)
THEN
2604 q_mean = sum(anal3d(ij,ilev,1:nmem,n)) / real(nmem,
rp)
2607 q_anal(m) = anal3d(ij,ilev,m,n) - q_mean
2608 q_sprd = q_sprd + q_anal(m)**2
2611 if ( q_mean > 0.0_rp )
then
2612 q_sprd = sqrt(q_sprd / real(nmem-1,
rp)) / q_mean
2613 IF(q_sprd > q_sprd_max)
THEN
2615 anal3d(ij,ilev,m,n) = q_mean + q_anal(m) * q_sprd_max / q_sprd
2628 n2nc = var_local_n2nc(
nv3d+n)
2629 n2n = var_local_n2n(
nv3d+n)
2631 if (relax_to_inflated_prior)
then
2638 if (trans_done(n2nc))
then
2640 IF(n2n <=
nv3d)
then
2641 if (infl_mul_adaptive)
then
2642 work2d(ij,n) = work3d(ij,ilev,n2n)
2645 work2dn(:,ij,n) = work3dn(:,ij,ilev,n2n)
2648 if (infl_mul_adaptive)
then
2649 work2d(ij,n) = work2d(ij,n2n-
nv3d)
2652 work2dn(:,ij,n) = work2dn(:,ij,n2n-
nv3d)
2658 if (ens_with_mdet)
then
2659 CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),
nv3d+n,hdxf,rdiag,rloc,dep,nobsl,depd=depd,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,
nv3d+1,ij))
2661 CALL obs_local(obs(:),rig1(ij),rjg1(ij),gues3d(ij,ilev,mmean,iv3d_p),hgt1(ij,ilev),
nv3d+n,hdxf,rdiag,rloc,dep,nobsl,nobsl_t=nobsl_t,cutd_t=cutd_t,srch_q0=search_q0(:,
nv3d+1,ij))
2663 IF(relax_alpha_spread /= 0.0d0)
THEN
2664 if (ens_with_mdet)
then
2665 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), &
2666 trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), &
2667 rdiag_wloc=.true.,infl_update=infl_mul_adaptive, &
2668 depd=depd,transmd=transmd(:,n2nc))
2670 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), &
2671 trans(:,:,n2nc),transm=transm(:,n2nc),pao=pa(:,:,n2nc), &
2672 rdiag_wloc=.true.,infl_update=infl_mul_adaptive)
2675 if (ens_with_mdet)
then
2676 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), &
2677 trans(:,:,n2nc),transm=transm(:,n2nc), &
2678 rdiag_wloc=.true.,infl_update=infl_mul_adaptive, &
2679 depd=depd,transmd=transmd(:,n2nc))
2681 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work2d(ij,n), &
2682 trans(:,:,n2nc),transm=transm(:,n2nc), &
2683 rdiag_wloc=.true.,infl_update=infl_mul_adaptive)
2686 trans_done(n2nc) = .true.
2688 work2dn(:,ij,n) = real(sum(nobsl_t,dim=1),
rp)
2694 IF(relax_alpha /= 0.0d0)
THEN
2695 CALL weight_rtpp(trans(:,:,n2nc),parm,transrlx)
2696 ELSE IF(relax_alpha_spread /= 0.0d0)
THEN
2697 IF(relax_spread_out)
THEN
2698 CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues2d(ij,:,n), &
2699 parm,transrlx,work2da(ij,n))
2701 CALL weight_rtps(trans(:,:,n2nc),pa(:,:,n2nc),gues2d(ij,:,n), &
2702 parm,transrlx,tmpinfl)
2705 transrlx = trans(:,:,n2nc)
2711 transrlx(k,m) = (transrlx(k,m) + transm(k,n2nc)) * beta
2713 transrlx(m,m) = transrlx(m,m) + (1.0_rp-beta)
2718 anal2d(ij,m,n) = gues2d(ij,mmean,n)
2720 anal2d(ij,m,n) = anal2d(ij,m,n) + gues2d(ij,k,n) * transrlx(k,m)
2725 if (ens_with_mdet)
then
2726 anal2d(ij,mmdet,n) = 0.0_rp
2728 anal2d(ij,mmdet,n) = anal2d(ij,mmdet,n) + gues2d(ij,k,n) * transmd(k,n2nc)
2730 anal2d(ij,mmdet,n) = gues2d(ij,mmdet,n) + anal2d(ij,mmdet,n) * beta
2738 deallocate (hdxf,rdiag,rloc,dep)
2739 if (ens_with_mdet)
then
2742 deallocate (trans,transm,transmd,pa)
2744 deallocate (n_merge,ic_merge)
2745 deallocate (search_q0)
2747 IF (
allocated(work3dg))
deallocate (work3dg)
2748 IF (
allocated(work2dg))
deallocate (work2dg)
2753 IF(infl_add > 0.0d0)
THEN
2755 if (infl_add_q_ratio)
then
2756 work3d(:,:,:) = gues3d(:,:,mmean,:)
2758 work3d(:,:,:) = 1.0d0
2761 allocate (addinfl_weight(nij1))
2762 if (infl_add_ref_only)
then
2763 addinfl_weight(:) = 0.0d0
2764 ic = ctype_elmtyp(uid_obs(id_radar_ref_obs), 22)
2767 ref_min_dist = 1.0d33
2769 do iob = obsgrd(ic)%ac_ext(0, 1) + 1, obsgrd(ic)%ac_ext(obsgrd(ic)%ngrdext_i, obsgrd(ic)%ngrdext_j)
2770 rdx = (rig1(ij) - obs(obsda_sort%set(iob))%ri(obsda_sort%idx(iob))) * dx
2771 rdy = (rjg1(ij) - obs(obsda_sort%set(iob))%rj(obsda_sort%idx(iob))) * dy
2772 rdxy = rdx*rdx + rdy*rdy
2773 if (rdxy < ref_min_dist)
then
2778 ref_min_dist = ref_min_dist / (hori_loc_ctype(ic) * hori_loc_ctype(ic))
2779 if (ref_min_dist <= dist_zero_fac_square)
then
2780 addinfl_weight(ij) = exp(-0.5d0 * ref_min_dist)
2785 addinfl_weight(:) = 1.0d0
2788 log_info(
"LETKF_system",*)
'Additive covariance inflation, parameter:', infl_add
2789 if (infl_add_shuffle)
then
2790 if (rank_ens== 0)
then
2791 call knuth_shuffle(nens, ishuf)
2793 call mpi_bcast(ishuf, nens, mpi_integer, 0, comm_ens, ierr)
2794 log_info(
"LETKF_system",*)
'shuffle members: on'
2795 log_info(
"LETKF_system",*)
'shuffle sequence:', ishuf
2800 if (infl_add_shuffle)
then
2805 if (n == iv3d_q .or. n == iv3d_qc .or. n == iv3d_qr .or. n == iv3d_qi .or. n == iv3d_qs .or. n == iv3d_qg)
then
2808 anal3d(i,k,m,n) = anal3d(i,k,m,n) &
2809 & + addi3d(i,k,mshuf,n) * infl_add * addinfl_weight(i) * work3d(i,k,n)
2815 anal3d(i,k,m,n) = anal3d(i,k,m,n) &
2816 & + addi3d(i,k,mshuf,n) * infl_add * addinfl_weight(i)
2824 if (infl_add_shuffle)
then
2830 anal2d(i,m,n) = anal2d(i,m,n) + addi2d(i,mshuf,n) * infl_add * addinfl_weight(i)
2835 deallocate (addinfl_weight)
2840 call gather_grd_mpi_all2all( 1, nens, anal3d(:,:,1:nens,:), anal2d(:,1:nens,:), v3dg, v2dg )
2845 u(k,i,j) = v3dg(k,i,j,iv3d_u )
2846 v(k,i,j) = v3dg(k,i,j,iv3d_v )
2847 w(k,i,j) = v3dg(k,i,j,iv3d_w )
2848 temp(k,i,j) = v3dg(k,i,j,iv3d_t )
2849 pres(k,i,j) = v3dg(k,i,j,iv3d_p )
2850 qv(k,i,j) = v3dg(k,i,j,iv3d_q )
2851 qc(k,i,j) = v3dg(k,i,j,iv3d_qc)
2852 qr(k,i,j) = v3dg(k,i,j,iv3d_qr)
2853 qi(k,i,j) = v3dg(k,i,j,iv3d_qi)
2854 qs(k,i,j) = v3dg(k,i,j,iv3d_qs)
2855 qg(k,i,j) = v3dg(k,i,j,iv3d_qg)
2860 if( letkf_debug_log )
then
2861 call monit_obs_mpi( obs_in_num, obs_in_format, v3dg, v2dg, monit_step=2 )
2876 average => statistics_average
2879 integer,
intent(in) :: pest_pmax
2881 real(
rp),
intent(inout) :: pest_var0(nens,pest_pmax)
2883 real(
rp) :: gues0d(nens+1,pest_pmax)
2884 real(
rp) :: anal0d(nens+1,pest_pmax)
2886 real(
rp) :: work0d(pest_pmax)
2888 real(
rp),
allocatable :: hdxf(:,:)
2889 real(
rp),
allocatable :: rdiag(:)
2890 real(
rp),
allocatable :: rloc(:)
2891 real(
rp),
allocatable :: dep(:)
2892 real(
rp),
allocatable :: depd(:)
2895 real(
rp),
allocatable :: trans(:,:)
2896 real(
rp),
allocatable :: transm(:)
2897 real(
rp),
allocatable :: transmd(:)
2898 real(
rp),
allocatable :: pa(:,:)
2899 real(
rp) :: transrlx(nmem,nmem)
2901 integer :: n,m,k,nobsl
2906 log_progress(*)
'data-assimilation / LETKF / parameter estimation'
2911 gues0d(m,n) = pest_var0(m,n)
2917 gues0d(mmean,n) = average( nmem, gues0d(1:nmem,n), undef )
2930 gues0d(m,n) = gues0d(m,n) - gues0d(mmean,n)
2934 allocate (hdxf(nobstotal,nmem))
2935 allocate (rdiag(nobstotal))
2936 allocate (rloc(nobstotal))
2937 allocate (dep(nobstotal))
2938 if (ens_with_mdet)
then
2939 allocate (depd(nobstotal))
2940 allocate (transmd(nmem))
2942 allocate (trans(nmem,nmem))
2943 allocate (transm(nmem))
2944 allocate (pa(nmem,nmem))
2957 if (ens_with_mdet)
then
2958 CALL obs_pest_etkf(hdxf, rdiag, rloc, dep, nobsl, depd=depd)
2960 CALL obs_pest_etkf(hdxf, rdiag, rloc, dep, nobsl)
2963 if (ens_with_mdet)
then
2964 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work0d(n), &
2965 trans(:,:),transm=transm(:),pao=pa(:,:), &
2967 depd=depd,transmd=transmd(:))
2969 CALL letkf_core(nmem,nobstotal,nobsl,hdxf,rdiag,rloc,dep,work0d(n), &
2970 trans(:,:),transm=transm(:),pao=pa(:,:), &
2975 CALL weight_rtps_const(trans(:,:),pa(:,:),gues0d(:,n),transrlx,tmpinfl)
2980 transrlx(k,m) = (transrlx(k,m) + transm(k)) * beta
2982 transrlx(m,m) = transrlx(m,m) + (1.0_rp-beta)
2987 anal0d(m,n) = gues0d(mmean,n)
2989 anal0d(m,n) = anal0d(m,n) + gues0d(k,n) * transrlx(k,m)
2994 if (ens_with_mdet)
then
2995 anal0d(mmdet,n) = 0.0_rp
2997 anal0d(mmdet,n) = anal0d(mmdet,n) + gues0d(k,n) * transmd(k)
2999 anal0d(mmdet,n) = gues0d(mmdet,n) + anal0d(mmdet,n) * beta
3006 pest_var0(m,n) = anal0d(m,n)
3010 deallocate (hdxf,rdiag,rloc,dep)
3011 if (ens_with_mdet)
then
3012 deallocate (depd,transmd)
3014 deallocate (trans,transm,pa)
3027 real(
rp),
intent(out) :: addi3d(:,:,:,:)
3028 real(
rp),
intent(out) :: addi2d(:,:,:)
3029 integer :: i, k, m, n
3031 call read_ens_mpi_addiinfl(addi3d, addi2d)
3033 call ensmean_grd(nens, nens, nij1, addi3d, addi2d)
3039 addi3d(i,k,m,n) = addi3d(i,k,m,n) - addi3d(i,k,mmean,n)
3048 addi2d(i,m,n) = addi2d(i,m,n) - addi2d(i,mmean,n)
3092 subroutine letkf_core(ne,nobs,nobsl,hdxb,rdiag,rloc,dep,parm_infl,trans,transm,pao,rdiag_wloc,infl_update,depd,transmd)
3099 INTEGER,
INTENT(IN) :: ne
3100 INTEGER,
INTENT(IN) :: nobs
3101 INTEGER,
INTENT(IN) :: nobsl
3102 REAL(
rp),
INTENT(IN) :: hdxb(1:nobs,1:ne)
3103 REAL(
rp),
INTENT(IN) :: rdiag(1:nobs)
3104 REAL(
rp),
INTENT(IN) :: rloc(1:nobs)
3105 REAL(
rp),
INTENT(IN) :: dep(1:nobs)
3106 REAL(
rp),
INTENT(INOUT) :: parm_infl
3107 REAL(
rp),
INTENT(OUT) :: trans(ne,ne)
3108 REAL(
rp),
INTENT(OUT),
OPTIONAL :: transm(ne)
3109 REAL(
rp),
INTENT(OUT),
OPTIONAL :: pao(ne,ne)
3110 LOGICAL,
INTENT(IN),
OPTIONAL :: rdiag_wloc
3111 LOGICAL,
INTENT(IN),
OPTIONAL :: infl_update
3112 REAL(
rp),
INTENT(IN),
OPTIONAL :: depd(1:nobs)
3113 REAL(
rp),
INTENT(OUT),
OPTIONAL :: transmd(ne)
3115 REAL(
rp) :: hdxb_rinv(nobsl,ne)
3116 REAL(
rp) :: eivec(ne,ne)
3117 REAL(
rp) :: eival(ne)
3118 REAL(
rp) :: pa(ne,ne)
3119 REAL(
rp) :: work1(ne,ne)
3120 REAL(
rp) :: work2(ne)
3121 REAL(
rp) :: work2d(ne)
3122 REAL(
rp) :: work3(ne)
3124 REAL(
rp) :: parm(4),sigma_o,gain
3125 REAL(
rp),
PARAMETER :: sigma_b = 0.04d0
3126 LOGICAL :: rdiag_wloc_
3127 LOGICAL :: infl_update_
3130 real(
rp) :: work1_evp(ne,ne)
3131 real(
rp) :: eivec_evp(ne,ne)
3132 real(
rp) :: eival_evp(ne)
3135 rdiag_wloc_ = .false.
3136 IF(
present(rdiag_wloc)) rdiag_wloc_ = rdiag_wloc
3137 infl_update_ = .false.
3138 IF(
present(infl_update)) infl_update_ = infl_update
3143 trans(i,i) = sqrt(parm_infl)
3145 IF (
PRESENT(transm))
THEN
3148 IF (
PRESENT(transmd))
THEN
3151 IF (
PRESENT(pao))
THEN
3154 pao(i,i) = parm_infl / real(ne-1,kind=
rp)
3162 IF(rdiag_wloc_)
THEN
3165 hdxb_rinv(i,j) = hdxb(i,j) / rdiag(i)
3171 hdxb_rinv(i,j) = hdxb(i,j) / rdiag(i) * rloc(i)
3180 CALL sgemm(
't',
'n',ne,ne,nobsl,1.0e0,hdxb_rinv,nobsl,hdxb(1:nobsl,:),nobsl,0.0e0,work1,ne)
3182 CALL dgemm(
't',
'n',ne,ne,nobsl,1.0d0,hdxb_rinv,nobsl,hdxb(1:nobsl,:),nobsl,0.0d0,work1,ne)
3195 rho = 1.0d0 / parm_infl
3197 work1(i,i) = work1(i,i) + real(ne-1,kind=
rp) * rho
3202 work1_evp = real( work1, kind=
rp )
3204 eival = real( eival_evp, kind=
rp )
3205 eivec = real( eivec_evp, kind=
rp )
3211 work1(i,j) = eivec(i,j) / eival(j)
3215 CALL sgemm(
'n',
't',ne,ne,ne,1.0e0,work1,ne,eivec,ne,0.0e0,pa,ne)
3217 CALL dgemm(
'n',
't',ne,ne,ne,1.0d0,work1,ne,eivec,ne,0.0d0,pa,ne)
3231 CALL sgemv(
't',nobsl,ne,1.0e0,hdxb_rinv,nobsl,dep,1,0.0e0,work2,1)
3233 CALL dgemv(
't',nobsl,ne,1.0d0,hdxb_rinv,nobsl,dep,1,0.0d0,work2,1)
3241 IF (
PRESENT(depd) .AND.
PRESENT(transmd))
THEN
3243 CALL sgemv(
't',nobsl,ne,1.0e0,hdxb_rinv,nobsl,depd,1,0.0e0,work2d,1)
3245 CALL dgemv(
't',nobsl,ne,1.0d0,hdxb_rinv,nobsl,depd,1,0.0d0,work2d,1)
3258 CALL sgemv(
'n',ne,ne,1.0e0,pa,ne,work2,1,0.0e0,work3,1)
3260 CALL dgemv(
'n',ne,ne,1.0d0,pa,ne,work2,1,0.0d0,work3,1)
3271 IF (
PRESENT(depd) .AND.
PRESENT(transmd))
THEN
3273 CALL sgemv(
'n',ne,ne,1.0e0,pa,ne,work2d,1,0.0e0,transmd,1)
3275 CALL dgemv(
'n',ne,ne,1.0d0,pa,ne,work2d,1,0.0d0,transmd,1)
3327 rho = sqrt( real(ne-1,kind=
rp) / eival(j) )
3329 work1(i,j) = eivec(i,j) * rho
3333 CALL sgemm(
'n',
't',ne,ne,ne,1.0e0,work1,ne,eivec,ne,0.0e0,trans,ne)
3335 CALL dgemm(
'n',
't',ne,ne,ne,1.0d0,work1,ne,eivec,ne,0.0d0,trans,ne)
3348 IF (
PRESENT(transm))
THEN
3353 trans(i,j) = trans(i,j) + work3(i)
3357 IF (
PRESENT(pao)) pao = pa
3359 IF (.NOT. infl_update_)
RETURN
3364 IF(rdiag_wloc_)
THEN
3366 parm(1) = parm(1) + dep(i)*dep(i)/rdiag(i)
3370 parm(1) = parm(1) + dep(i)*dep(i)/rdiag(i) * rloc(i)
3375 parm(2) = parm(2) + hdxb_rinv(i,j) * hdxb(i,j)
3378 parm(2) = parm(2) / real(ne-1,kind=
rp)
3379 parm(3) = sum(rloc(1:nobsl))
3380 parm(4) = (parm(1)-parm(3))/parm(2) - parm_infl
3382 sigma_o = 2.0d0/parm(3)*((parm_infl*parm(2)+parm(3))/parm(2))**2
3383 gain = sigma_b**2 / (sigma_o + sigma_b**2)
3384 parm_infl = parm_infl + gain * parm(4)
3388 end subroutine letkf_core
3393 subroutine get_nobs(cfile,nrec,nn)
3396 character(*),
intent(in) :: cfile
3397 integer,
intent(in) :: nrec
3398 integer,
intent(out) :: nn
3406 inquire(file=cfile,exist=ex)
3408 open(unit=iunit,file=cfile,form=
'unformatted',access=
'stream')
3409 inquire(unit=iunit, size=sz)
3410 if (mod(sz,
sp * (nrec+2) ) /= 0)
then
3411 log_error(
'get_nobs',*)
': reading error -- skipped: ', cfile
3415 nn = sz / (
sp * (nrec+2) )
3418 log_info(
'get_nobs',*)
'file not found -- skipped: ', cfile
3422 end subroutine get_nobs
3424 subroutine get_nobs_radar(cfile,nn,radarlon,radarlat,radarz)
3426 character(*),
intent(in) :: cfile
3427 integer,
intent(out) :: nn
3428 real(
rp),
intent(out) :: radarlon,radarlat,radarz
3436 if(radar_obs_4d)
then
3448 inquire(file=cfile,exist=ex)
3450 open(unit=iunit,file=cfile,form=
'unformatted',access=
'stream')
3451 read(unit=iunit,iostat=ios) tmp
3452 read(unit=iunit,iostat=ios) tmp
3454 log_error(
'get_nobs_radar',*)
': reading error -- skipped: ', cfile
3458 radarlon = real(tmp,kind=
rp)
3459 read(unit=iunit,iostat=ios) tmp
3460 read(unit=iunit,iostat=ios) tmp
3461 read(unit=iunit,iostat=ios) tmp
3463 log_error(
'get_nobs_radar',*)
': reading error -- skipped: ', cfile
3467 radarlat = real(tmp,kind=
rp)
3468 read(unit=iunit,iostat=ios) tmp
3469 read(unit=iunit,iostat=ios) tmp
3470 read(unit=iunit,iostat=ios) tmp
3472 log_error(
'get_nobs_radar',*)
': reading error -- skipped: ', cfile
3476 radarz = real(tmp,kind=
rp)
3477 read(unit=iunit,iostat=ios) tmp
3480 inquire(unit=iunit, size=sz)
3481 sz = sz -
sp * (1+2) * 3
3482 if (mod(sz,
sp * (nrec+2)) /= 0)
then
3483 log_error(
'get_nobs_radar',*)
': reading error -- skipped: ', cfile
3487 nn = sz / (
sp * (nrec+2))
3491 log_info(
'get_nobs_radar',*)
'file not found -- skipped: ', cfile
3495 end subroutine get_nobs_radar
3497 subroutine read_obs(cfile,obs)
3501 mapprojection_lonlat2xy
3504 character(*),
intent(in) :: cfile
3505 type(
obs_info),
intent(inout) :: obs
3507 real(
rp) :: x_rp, y_rp
3512 open(unit=iunit,file=cfile,form=
'unformatted',access=
'stream')
3514 read(unit=iunit) tmp
3515 read(unit=iunit) wk(:)
3516 read(unit=iunit) tmp
3517 select case(nint(wk(1)))
3519 wk(4) = wk(4) * 100.0
3521 wk(4) = wk(4) * 100.0
3523 wk(4) = wk(4) * 100.0
3525 wk(4) = wk(4) * 100.0
3527 wk(4) = wk(4) * 100.0
3529 wk(5) = wk(5) * 100.0
3530 wk(6) = wk(6) * 100.0
3532 wk(4) = wk(4) * 100.0
3533 wk(5) = wk(5) * 0.01
3534 wk(6) = wk(6) * 0.01
3536 wk(4) = wk(4) * 100.0
3537 wk(5) = wk(5) * 100.0
3538 wk(6) = real( obserr_tcp,kind=
sp )
3540 call mapprojection_lonlat2xy( real(wk(2),kind=
rp)*pi/180.0_rp,&
3541 real(wk(3),kind=
rp)*pi/180.0_rp,&
3543 wk(4) = wk(4) * 100.0
3544 wk(5) = real( x_rp, kind=
sp )
3545 wk(6) = real( obserr_tcx, kind=
sp )
3547 call mapprojection_lonlat2xy( real(wk(2),kind=
rp)*pi/180.0_rp,&
3548 real(wk(3),kind=
rp)*pi/180.0_rp,&
3550 wk(4) = wk(4) * 100.0
3551 wk(5) = real( y_rp, kind=
sp )
3552 wk(6) = real( obserr_tcy, kind=
sp )
3554 obs%elm(n) = nint( wk(1) )
3555 obs%lon(n) = real( wk(2), kind=
rp )
3556 obs%lat(n) = real( wk(3), kind=
rp )
3557 obs%lev(n) = real( wk(4), kind=
rp )
3558 obs%dat(n) = real( wk(5), kind=
rp )
3559 obs%err(n) = real( wk(6), kind=
rp )
3560 obs%typ(n) = nint( wk(7) )
3561 obs%dif(n) = real( wk(8), kind=
rp )
3566 end subroutine read_obs
3568 subroutine read_obs_radar(cfile,obs)
3570 character(*),
intent(in) :: cfile
3571 type(
obs_info),
intent(inout) :: obs
3575 integer :: n,iunit,ios
3577 if(radar_obs_4d)
then
3583 open(unit=iunit,file=cfile,form=
'unformatted',access=
'stream')
3584 read(unit=iunit,iostat=ios) tmp
3585 read(unit=iunit,iostat=ios) tmp
3586 read(unit=iunit,iostat=ios) tmp
3587 read(unit=iunit,iostat=ios) tmp
3588 read(unit=iunit,iostat=ios) tmp
3589 read(unit=iunit,iostat=ios) tmp
3590 read(unit=iunit,iostat=ios) tmp
3591 read(unit=iunit,iostat=ios) tmp
3592 read(unit=iunit,iostat=ios) tmp
3594 read(unit=iunit) tmp
3595 read(unit=iunit) wk(1:nrec)
3596 read(unit=iunit) tmp
3597 obs%elm(n) = nint(wk(1))
3598 obs%lon(n) = real(wk(2),kind=
rp)
3599 obs%lat(n) = real(wk(3),kind=
rp)
3600 obs%lev(n) = real(wk(4),kind=
rp)
3601 obs%dat(n) = real(wk(5),kind=
rp)
3602 obs%err(n) = real(wk(6),kind=
rp)
3604 if(radar_obs_4d)
then
3605 obs%dif(n) = real(wk(8),kind=
rp)
3613 end subroutine read_obs_radar
3615 subroutine read_obs_radar_toshiba_pawr( obs, cfile )
3625 character(len=*),
intent(in) :: cfile
3630 integer,
parameter :: n_type = 3
3631 character(len=4),
parameter :: file_type_sfx(n_type) = &
3632 (/
'.ze ',
'.vr ',
'.qcf'/)
3633 logical,
parameter :: input_is_dbz = .true.
3636 real(kind=c_float),
allocatable,
save :: rtdat(:, :, :, :)
3637 real(kind=c_float),
allocatable,
save :: az(:, :, :)
3638 real(kind=c_float),
allocatable,
save :: el(:, :, :)
3639 integer :: j, ierr, ierr2
3640 character(len=3) :: fname
3643 real(
rp),
allocatable :: ze(:, :, :), vr(:, :, :), qcflag(:, :, :), attenuation(:, :, :), rrange(:)
3644 real(
rp),
allocatable :: radlon(:, :, :), radlat(:, :, :), radz(:, :, :)
3645 real(
rp),
allocatable :: lon(:), lat(:), z(:)
3646 integer(8),
allocatable :: grid_index(:), grid_count_ze(:), grid_count_vr(:)
3647 real(
rp),
allocatable :: grid_ze(:), grid_vr(:)
3648 real(
rp),
allocatable :: grid_lon_ze(:), grid_lat_ze(:), grid_z_ze(:)
3649 real(
rp),
allocatable :: grid_lon_vr(:), grid_lat_vr(:), grid_z_vr(:)
3651 character(len=1024) :: input_fname(n_type)
3653 real(
rp) :: dlon, dlat
3656 integer,
save :: na, nr, ne
3657 real(
rp),
save :: lon0, lat0, z0
3658 real(
rp),
save :: missing
3659 integer,
save :: range_res
3661 real(
rp) :: max_obs_ze , min_obs_ze , max_obs_vr , min_obs_vr
3662 integer :: nobs_ze, nobs_vr
3663 integer(8) :: idx, n, n_ref
3665 integer,
parameter :: int1 = selected_int_kind(1)
3666 integer(kind = int1) :: tmp_qcf, valid_qcf
3668 integer,
parameter :: qcf_mask(8)=(/ 0, 1, 1, 1, 1, 0, 0, 0 /)
3670 integer::qcf_count(0:255)
3672 character(len=90) :: plotname
3673 character(len=19) :: timelabel
3675 integer :: ii, jj, kk
3677 real(
sp),
allocatable :: ref3d(:,:,:)
3678 real(
sp),
allocatable :: vr3d(:,:,:)
3679 character(len=255) :: filename
3680 integer :: irec, iunit, iolen
3683 integer :: nlon, nlat, nlev
3685 radar_so_size_hori = max( real( dx, kind=
rp ), radar_so_size_hori )
3686 radar_so_size_hori = max( real( dy, kind=
rp ), radar_so_size_hori )
3688 if( .not.
allocated(rtdat) )
allocate(rtdat(
rdim,
azdim,
eldim, n_type))
3689 if( .not.
allocated(az) )
allocate(az(
azdim,
eldim, n_type))
3690 if( .not.
allocated(el) )
allocate(el(
azdim,
eldim, n_type))
3693 input_fname(j) = trim(cfile)
3694 call str_replace(input_fname(j),
'.<type>', trim(file_type_sfx(j)), pos)
3695 ierr =
read_toshiba( input_fname(j), hd(j), az(:,:,j), el(:,:,j), rtdat(:,:,:,j) )
3696 if( ierr /= 0 )
then
3697 log_info(
"LETKF_obs_readfile",*)
'Warning: File (',trim(input_fname(j)),
') cannot be read. Skip.'
3704 lon0 = hd(1)%longitude
3705 lat0 = hd(1)%latitude
3708 missing = real( hd(1)%mesh_offset, kind=
rp )
3709 range_res = hd(1)%range_res
3711 nr = hd(1)%range_num
3712 na = hd(1)%sector_num
3715 call jst2utc( hd(1)%s_yr, hd(1)%s_mn, hd(1)%s_dy, hd(1)%s_hr, hd(1)%s_mi, hd(1)%s_sc, 0.0_dp, utime_obs )
3717 allocate( ze(na, nr, ne) )
3718 allocate( vr(na, nr, ne) )
3719 allocate( qcflag(na, nr, ne) )
3720 allocate( attenuation(na, nr, ne) )
3724 if(qcf_mask(j) > 0) valid_qcf = ibset(valid_qcf, j - 1)
3730 ze(ia,ir,ie) = rtdat(ir,ia,ie,1)
3731 vr(ia,ir,ie) = rtdat(ir,ia,ie,2)
3733 if( vr(ia,ir,ie) > radar_max_abs_vr .or. vr(ia,ir,ie) < -radar_max_abs_vr ) vr(ia,ir,ie) = missing
3735 tmp_qcf = int(rtdat(ir, ia, ie, 3), int1)
3736 if(iand(valid_qcf, tmp_qcf) == 0)
then
3737 qcflag(ia, ir, ie) = 0.0_rp
3739 qcflag(ia, ir, ie) = 1000.0_rp
3741 attenuation(ia, ir, ie) = 1.0_rp
3747 allocate(rrange(nr))
3749 rrange(ir) = (dble(ir) - 0.5_rp) * range_res
3752 allocate( radlon(na, nr, ne) )
3753 allocate( radlat(na, nr, ne) )
3754 allocate( radz(na, nr, ne) )
3756 call radar_georeference( lon0, lat0, z0, na, nr, ne, &
3757 real( az(:, 1, 1), kind=
rp ), &
3759 real( el(1, :, 1), kind=
rp ), &
3760 radlon, radlat, radz )
3762 call define_grid( lon0, lat0, nr, rrange, rrange(nr), radar_zmax, &
3763 radar_so_size_hori, radar_so_size_hori, radar_so_size_vert, &
3764 dlon, dlat, nlon, nlat, nlev, lon, lat, z )
3767 qcflag, attenuation, &
3768 nlon, nlat, nlev, lon, lat, z, dlon, dlat, radar_so_size_vert, &
3769 missing, input_is_dbz, &
3771 nobs_sp, grid_index, &
3772 grid_ze, grid_lon_ze, grid_lat_ze, grid_z_ze, grid_count_ze, &
3773 grid_vr, grid_lon_vr, grid_lat_vr, grid_z_vr, grid_count_vr )
3775 if(
allocated( ze ) )
deallocate(ze)
3776 if(
allocated( vr ) )
deallocate(vr)
3777 if(
allocated( qcflag ) )
deallocate(qcflag)
3778 if(
allocated( attenuation ) )
deallocate(attenuation)
3779 if(
allocated( rrange ) )
deallocate(rrange)
3790 ii = nint( ( grid_lon_ze(idx) - lon(1) ) / dlon ) + 1
3791 jj = nint( ( grid_lat_ze(idx) - lat(1) ) / dlat ) + 1
3792 kk = nint( ( grid_z_ze(idx) - z(1) ) / radar_so_size_vert ) + 1
3794 if ( mod(ii, radar_thin_hori) /= 0 .or. mod(jj, radar_thin_hori) /= 0 .or. &
3795 mod(kk, radar_thin_vert) /= 0 ) cycle
3797 if (grid_count_ze(idx) > 0)
then
3798 obs%nobs = obs%nobs + 1
3801 if ( grid_ze(idx) > min_radar_ref .and. grid_z_ze(idx) < radar_zmax )
then
3802 obs_ref%nobs = obs_ref%nobs + 1
3805 if (grid_count_vr(idx) > 0)
then
3806 obs%nobs = obs%nobs + 1
3810 call obs_info_allocate( obs, extended=.true. )
3811 call obs_info_allocate( obs_ref, extended=.true. )
3817 min_obs_ze = huge(1.0_rp)
3818 max_obs_ze = -huge(1.0_rp)
3819 min_obs_vr = huge(1.0_rp)
3820 max_obs_vr = -huge(1.0_rp)
3823 ii = nint( ( grid_lon_ze(idx) - lon(1) ) / dlon ) + 1
3824 jj = nint( ( grid_lat_ze(idx) - lat(1) ) / dlat ) + 1
3825 kk = nint( ( grid_z_ze(idx) - z(1) ) / radar_so_size_vert ) + 1
3828 if ( mod(ii, radar_thin_hori) /= 0 .or. mod(jj, radar_thin_hori) /= 0 .or. &
3829 mod(kk, radar_thin_vert) /= 0 ) cycle
3831 if (grid_count_ze(idx) > 0)
then
3833 obs%elm(n) = id_radar_ref_obs
3834 obs%lon(n) = grid_lon_ze(idx)
3835 obs%lat(n) = grid_lat_ze(idx)
3836 obs%lev(n) = grid_z_ze(idx)
3837 obs%dat(n) = grid_ze(idx)
3839 if ( radar_bias_cor_rain .and. grid_ze(idx) > min_radar_ref )
then
3840 obs%dat(n) = grid_ze(idx) * radar_bias_rain_const
3841 elseif ( radar_bias_cor_clr .and. grid_ze(idx) < min_radar_ref )
then
3842 obs%dat(n) = grid_ze(idx) * radar_bias_clr_const
3844 obs%err(n) = obserr_radar_ref
3847 nobs_ze = nobs_ze + 1
3848 if (grid_ze(idx) > max_obs_ze) max_obs_ze = grid_ze(idx)
3849 if (grid_ze(idx) < min_obs_ze) min_obs_ze = grid_ze(idx)
3851 if ( grid_ze(idx) > min_radar_ref .and. grid_z_ze(idx) < radar_zmax )
then
3853 obs_ref%elm(n_ref) = id_radar_ref_obs
3854 obs_ref%lon(n_ref) = grid_lon_ze(idx)
3855 obs_ref%lat(n_ref) = grid_lat_ze(idx)
3856 obs_ref%lev(n_ref) = grid_z_ze(idx)
3857 if ( radar_bias_cor_rain )
then
3859 obs_ref%dat(n_ref) = grid_ze(idx) * radar_bias_rain_const
3861 obs_ref%dat(n_ref) = grid_ze(idx)
3866 if (grid_count_vr(idx) > 0)
then
3868 obs%elm(n) = id_radar_vr_obs
3869 obs%lon(n) = grid_lon_ze(idx)
3870 obs%lat(n) = grid_lat_ze(idx)
3871 obs%lev(n) = grid_z_ze(idx)
3872 obs%dat(n) = grid_vr(idx)
3873 obs%err(n) = obserr_radar_vr
3876 nobs_vr = nobs_vr + 1
3877 if (grid_vr(idx) > max_obs_vr) max_obs_vr = grid_vr(idx)
3878 if (grid_vr(idx) < min_obs_vr) min_obs_vr = grid_vr(idx)
3882 if(
allocated(radlon) )
deallocate(radlon)
3883 if(
allocated(radlat) )
deallocate(radlat)
3884 if(
allocated(radz ) )
deallocate(radz)
3887 call obs_info_deallocate( obs_ref )
3889 if(
allocated( grid_index ) )
deallocate(grid_index)
3890 if(
allocated( grid_ze ) )
deallocate(grid_ze)
3891 if(
allocated( grid_lon_ze ) )
deallocate(grid_lon_ze)
3892 if(
allocated( grid_lat_ze ) )
deallocate(grid_lat_ze)
3893 if(
allocated( grid_z_ze ) )
deallocate(grid_z_ze)
3894 if(
allocated( grid_count_ze ) )
deallocate(grid_count_ze)
3895 if(
allocated( grid_vr ) )
deallocate(grid_vr)
3896 if(
allocated( grid_lon_vr ) )
deallocate(grid_lon_vr)
3897 if(
allocated( grid_lat_vr ) )
deallocate(grid_lat_vr)
3898 if(
allocated( grid_z_vr ) )
deallocate(grid_z_vr)
3899 if(
allocated( grid_count_vr ) )
deallocate(grid_count_vr)
3902 end subroutine read_obs_radar_toshiba_pawr
3904 subroutine read_obs_radar_toshiba_mp_pawr( obs, cfile, maskfile )
3914 character(len=*),
intent(in) :: cfile
3915 character(len=*),
intent(in) :: maskfile
3920 integer,
parameter :: n_type = 2
3921 character(len=4),
parameter :: file_type_sfx(n_type) = (/
'.ze',
'.vr'/)
3922 logical,
parameter :: input_is_dbz = .true.
3923 integer,
parameter :: opt_verbose = 0
3927 integer(4),
save :: shadow_na, shadow_ne
3928 integer(4),
allocatable,
save :: tmpshadow(:)
3929 integer(2),
allocatable,
save :: shadow(:,:)
3930 real(8),
save :: shadow_del_az
3933 real(kind=c_float),
allocatable,
save :: rtdat(:, :, :, :)
3934 real(kind=c_float),
allocatable,
save :: az(:, :, :)
3935 real(kind=c_float),
allocatable,
save :: el(:, :, :)
3936 integer :: j, ierr, ierr2
3937 character(len=3) :: fname
3940 real(
rp),
allocatable :: ze(:, :, :), vr(:, :, :), qcflag(:, :, :), attenuation(:, :, :), rrange(:)
3941 real(
rp),
allocatable :: radlon(:, :, :), radlat(:, :, :), radz(:, :, :)
3942 real(
rp),
allocatable :: lon(:), lat(:), z(:)
3943 integer(8),
allocatable :: grid_index(:), grid_count_ze(:), grid_count_vr(:)
3944 real(
rp),
allocatable :: grid_ze(:), grid_vr(:)
3945 real(
rp),
allocatable :: grid_lon_ze(:), grid_lat_ze(:), grid_z_ze(:)
3946 real(
rp),
allocatable :: grid_lon_vr(:), grid_lat_vr(:), grid_z_vr(:)
3948 character(len=1024) :: input_fname(n_type)
3950 real(
rp) :: dlon, dlat
3953 integer,
save :: na, nr, ne
3954 real(
rp),
save :: lon0, lat0, z0
3955 real(
rp),
save :: missing
3956 integer,
save :: range_res
3958 real(
rp) :: max_obs_ze , min_obs_ze , max_obs_vr , min_obs_vr
3959 integer :: nobs_ze, nobs_vr
3960 integer(8) :: idx, n, n_ref
3962 integer,
parameter :: int1 = selected_int_kind(1)
3963 integer(kind = int1) :: tmp_qcf, valid_qcf
3965 integer,
parameter :: qcf_mask(8)=(/ 0, 1, 1, 1, 1, 0, 0, 0 /)
3967 integer::qcf_count(0:255)
3969 character(len=90) :: plotname
3970 character(len=19) :: timelabel
3972 integer :: ii, jj, kk
3974 real(
sp),
allocatable :: ref3d(:,:,:)
3975 real(
sp),
allocatable :: vr3d(:,:,:)
3976 character(len=255) :: filename
3977 integer :: irec, iunit, iolen
3980 integer :: nlon, nlat, nlev
3982 radar_so_size_hori = max( real( dx, kind=
rp ), radar_so_size_hori )
3983 radar_so_size_hori = max( real( dy, kind=
rp ), radar_so_size_hori )
3985 if( .not.
allocated(rtdat) )
allocate(rtdat(
rdim,
azdim,
eldim, n_type))
3986 if( .not.
allocated(az) )
allocate(az(
azdim,
eldim, n_type))
3987 if( .not.
allocated(el) )
allocate(el(
azdim,
eldim, n_type))
3989 if( trim(maskfile) /=
'' .and. .NOT.
allocated( shadow ) )
then
3991 open(iunit, file=trim(maskfile), status=
"old", access=
"stream", form=
"unformatted")
3992 read(iunit,iostat=ios) shadow_na, shadow_ne
3993 allocate( shadow( shadow_na, shadow_ne ) )
3995 read(iunit,iostat=ios) shadow
3998 write(6,
'(3A)')
'file ',trim(maskfile) ,
' not found or unsupported format.'
4004 input_fname(j) = trim(cfile)
4005 call str_replace(input_fname(j),
'.<type>', trim(file_type_sfx(j)), pos)
4006 ierr =
read_toshiba_mpr( input_fname(j), opt_verbose, hd(j), az(:,:,j), el(:,:,j), rtdat(:,:,:,j) )
4007 if( ierr /= 0 )
then
4008 log_info(
"LETKF_obs_readfile",*)
'Warning: File (',trim(input_fname(j)),
') cannot be read. Skip.'
4015 lon0 = hd(1)%longitude
4016 lat0 = hd(1)%latitude
4019 missing = real( hd(1)%mesh_offset, kind=
rp )
4020 range_res = hd(1)%range_res
4022 nr = hd(1)%range_num
4026 call jst2utc( hd(1)%s_yr, hd(1)%s_mn, hd(1)%s_dy, hd(1)%s_hr, hd(1)%s_mi, hd(1)%s_sc, 0.0_dp, utime_obs )
4028 allocate( ze(na, nr, ne) )
4029 allocate( vr(na, nr, ne) )
4030 allocate( qcflag(na, nr, ne) )
4031 allocate( attenuation(na, nr, ne) )
4033 shadow_del_az = 360.0_rp / shadow_na
4034 if ( trim(maskfile) /=
'' .and. .not.
allocated(tmpshadow) )
then
4035 allocate(tmpshadow(na))
4041 ze(ia,ir,ie) = rtdat(ir,ia,ie,1)
4042 vr(ia,ir,ie) = rtdat(ir,ia,ie,2)
4044 if( vr(ia,ir,ie) > radar_max_abs_vr .or. vr(ia,ir,ie) < -radar_max_abs_vr ) vr(ia,ir,ie) = missing
4050 qcflag(ia,ir,ie) = 0.0_rp
4053 if( trim(maskfile) /=
'' .and.
allocated(shadow) )
then
4054 if( ie <= shadow_ne )
then
4056 tmpshadow(ia) = shadow( min(shadow_na,nint(az(ia, ie, 1) / shadow_del_az) + 1), ie )
4060 if( tmpshadow(ia) /= 0 .and. ir >= tmpshadow(ia) )
then
4061 qcflag(ia, ir, ie) = 1000.0_rp
4070 attenuation(ia, ir, ie) = 1.0_rp
4076 allocate(rrange(nr))
4078 rrange(ir) = (dble(ir) - 0.5_rp) * range_res
4081 allocate( radlon(na, nr, ne) )
4082 allocate( radlat(na, nr, ne) )
4083 allocate( radz(na, nr, ne) )
4085 call radar_georeference( lon0, lat0, z0, na, nr, ne, &
4086 real( az(:, 1, 1), kind=
rp ), &
4088 real( el(1, :, 1), kind=
rp ), &
4089 radlon, radlat, radz )
4091 call define_grid( lon0, lat0, nr, rrange, rrange(nr), radar_zmax, &
4092 radar_so_size_hori, radar_so_size_hori, radar_so_size_vert, &
4093 dlon, dlat, nlon, nlat, nlev, lon, lat, z )
4096 qcflag, attenuation, &
4097 nlon, nlat, nlev, lon, lat, z, dlon, dlat, radar_so_size_vert, &
4098 missing, input_is_dbz, &
4100 nobs_sp, grid_index, &
4101 grid_ze, grid_lon_ze, grid_lat_ze, grid_z_ze, grid_count_ze, &
4102 grid_vr, grid_lon_vr, grid_lat_vr, grid_z_vr, grid_count_vr )
4104 if(
allocated( ze ) )
deallocate(ze)
4105 if(
allocated( vr ) )
deallocate(vr)
4106 if(
allocated( qcflag ) )
deallocate(qcflag)
4107 if(
allocated( attenuation ) )
deallocate(attenuation)
4108 if(
allocated( rrange ) )
deallocate(rrange)
4119 ii = nint( ( grid_lon_ze(idx) - lon(1) ) / dlon ) + 1
4120 jj = nint( ( grid_lat_ze(idx) - lat(1) ) / dlat ) + 1
4121 kk = nint( ( grid_z_ze(idx) - z(1) ) / radar_so_size_vert ) + 1
4123 if ( mod(ii, radar_thin_hori) /= 0 .or. mod(jj, radar_thin_hori) /= 0 .or. &
4124 mod(kk, radar_thin_vert) /= 0 ) cycle
4126 if (grid_count_ze(idx) > 0)
then
4127 obs%nobs = obs%nobs + 1
4130 if ( grid_ze(idx) > min_radar_ref .and. grid_z_ze(idx) < radar_zmax )
then
4131 obs_ref%nobs = obs_ref%nobs + 1
4134 if (grid_count_vr(idx) > 0)
then
4135 obs%nobs = obs%nobs + 1
4139 call obs_info_allocate( obs, extended=.true. )
4140 call obs_info_allocate( obs_ref, extended=.true. )
4146 min_obs_ze = huge(1.0_rp)
4147 max_obs_ze = -huge(1.0_rp)
4148 min_obs_vr = huge(1.0_rp)
4149 max_obs_vr = -huge(1.0_rp)
4152 ii = nint( ( grid_lon_ze(idx) - lon(1) ) / dlon ) + 1
4153 jj = nint( ( grid_lat_ze(idx) - lat(1) ) / dlat ) + 1
4154 kk = nint( ( grid_z_ze(idx) - z(1) ) / radar_so_size_vert ) + 1
4157 if ( mod(ii, radar_thin_hori) /= 0 .or. mod(jj, radar_thin_hori) /= 0 .or. &
4158 mod(kk, radar_thin_vert) /= 0 ) cycle
4160 if (grid_count_ze(idx) > 0)
then
4162 obs%elm(n) = id_radar_ref_obs
4163 obs%lon(n) = grid_lon_ze(idx)
4164 obs%lat(n) = grid_lat_ze(idx)
4165 obs%lev(n) = grid_z_ze(idx)
4166 obs%dat(n) = grid_ze(idx)
4168 if ( radar_bias_cor_rain .and. grid_ze(idx) > min_radar_ref )
then
4169 obs%dat(n) = grid_ze(idx) * radar_bias_rain_const
4170 elseif ( radar_bias_cor_clr .and. grid_ze(idx) < min_radar_ref )
then
4171 obs%dat(n) = grid_ze(idx) * radar_bias_clr_const
4173 obs%err(n) = obserr_radar_ref
4176 nobs_ze = nobs_ze + 1
4177 if (grid_ze(idx) > max_obs_ze) max_obs_ze = grid_ze(idx)
4178 if (grid_ze(idx) < min_obs_ze) min_obs_ze = grid_ze(idx)
4180 if ( grid_ze(idx) > min_radar_ref .and. grid_z_ze(idx) < radar_zmax )
then
4182 obs_ref%elm(n_ref) = id_radar_ref_obs
4183 obs_ref%lon(n_ref) = grid_lon_ze(idx)
4184 obs_ref%lat(n_ref) = grid_lat_ze(idx)
4185 obs_ref%lev(n_ref) = grid_z_ze(idx)
4186 if ( radar_bias_cor_rain )
then
4188 obs_ref%dat(n_ref) = grid_ze(idx) * radar_bias_rain_const
4190 obs_ref%dat(n_ref) = grid_ze(idx)
4195 if (grid_count_vr(idx) > 0)
then
4197 obs%elm(n) = id_radar_vr_obs
4198 obs%lon(n) = grid_lon_ze(idx)
4199 obs%lat(n) = grid_lat_ze(idx)
4200 obs%lev(n) = grid_z_ze(idx)
4201 obs%dat(n) = grid_vr(idx)
4202 obs%err(n) = obserr_radar_vr
4205 nobs_vr = nobs_vr + 1
4206 if (grid_vr(idx) > max_obs_vr) max_obs_vr = grid_vr(idx)
4207 if (grid_vr(idx) < min_obs_vr) min_obs_vr = grid_vr(idx)
4211 if(
allocated(radlon) )
deallocate(radlon)
4212 if(
allocated(radlat) )
deallocate(radlat)
4213 if(
allocated(radz ) )
deallocate(radz)
4216 call obs_info_deallocate( obs_ref )
4218 if(
allocated( grid_index ) )
deallocate(grid_index)
4219 if(
allocated( grid_ze ) )
deallocate(grid_ze)
4220 if(
allocated( grid_lon_ze ) )
deallocate(grid_lon_ze)
4221 if(
allocated( grid_lat_ze ) )
deallocate(grid_lat_ze)
4222 if(
allocated( grid_z_ze ) )
deallocate(grid_z_ze)
4223 if(
allocated( grid_count_ze ) )
deallocate(grid_count_ze)
4224 if(
allocated( grid_vr ) )
deallocate(grid_vr)
4225 if(
allocated( grid_lon_vr ) )
deallocate(grid_lon_vr)
4226 if(
allocated( grid_lat_vr ) )
deallocate(grid_lat_vr)
4227 if(
allocated( grid_z_vr ) )
deallocate(grid_z_vr)
4228 if(
allocated( grid_count_vr ) )
deallocate(grid_count_vr)
4231 end subroutine read_obs_radar_toshiba_mp_pawr
4245 subroutine str_replace(str, oldsub, newsub, pos)
4249 character(len=*),
intent(inout) :: str
4250 character(len=*),
intent(in) :: oldsub
4251 character(len=*),
intent(in) :: newsub
4252 integer,
intent(out) :: pos
4253 integer :: str_lent, oldsub_len, newsub_len, shift
4256 str_lent = len_trim(str)
4257 oldsub_len = len(oldsub)
4258 newsub_len = len(newsub)
4260 pos = index(str, oldsub)
4262 shift = newsub_len - oldsub_len
4264 if (str_lent+shift > len(str))
then
4265 log_error(
'str_replace', *)
"The length of 'str' string is not enough for substitution."
4268 str(pos+oldsub_len:str_lent+shift) = adjustr(str(pos+oldsub_len:str_lent+shift))
4269 else if (shift < 0)
then
4270 str(pos+newsub_len:pos+oldsub_len-1) = repeat(
' ', 0-shift)
4271 str(pos+newsub_len:str_lent) = adjustl(str(pos+newsub_len:str_lent))
4273 str(pos:pos+newsub_len-1) = newsub
4277 end subroutine str_replace
4279 subroutine jst2utc(jyear, jmonth, jday, jhour, jminute, jsecond, jtime_ms, utime)
4286 integer,
intent(in) :: jyear, jmonth, jday
4287 integer,
intent(in) :: jhour, jminute, jsecond
4288 real(
dp) :: jtime_ms
4289 integer,
intent(out) :: utime(6)
4292 real(
dp) :: abssec, utime_ms
4307 abssec = abssec - real(3600*9, kind=
dp)
4319 end subroutine jst2utc
4332 subroutine radar_georeference(lon0, lat0, z0, na, nr, ne, azimuth, rrange, elevation, radlon, radlat, radz)
4339 real(
rp),
parameter :: ke = 4.0_rp / 3.0_rp
4343 real(
rp),
intent(in) :: lon0, lat0, z0
4344 integer,
intent(in) :: na, nr, ne
4345 real(
rp),
intent(in) :: azimuth(na), rrange(nr), elevation(ne)
4346 real(
rp),
intent(out) :: radlon(na, nr, ne), radlat(na, nr, ne)
4347 real(
rp),
intent(out) :: radz(na, nr, ne)
4349 real(
rp) sin_elev_ke_re_2, cos_elev_div_ke_re, tmpdist
4350 real(
rp) cdist, sdist, sinll1, cosll1, sinll1_cdist, cosll1_sdist, cosll1_cdist, sinll1_sdist
4351 real(
rp) :: azimuth_rad, sin_azim(na), cos_azim(na)
4352 integer :: ia, ir, ie
4357 sinll1 = sin(lat0 * d2r)
4358 cosll1 = cos(lat0 * d2r)
4360 azimuth_rad = azimuth(ia) * d2r
4361 sin_azim(ia) = sin(azimuth_rad)
4362 cos_azim(ia) = cos(azimuth_rad)
4366 sin_elev_ke_re_2 = sin(elevation(ie) * d2r) * ke_re * 2
4367 cos_elev_div_ke_re = cos(elevation(ie) * d2r) / ke_re
4370 radz(:, ir, ie) = z0 + sqrt(rrange(ir) * (rrange(ir) + sin_elev_ke_re_2) + ke_re * ke_re) - ke_re
4371 tmpdist = ke * asin(rrange(ir) * cos_elev_div_ke_re)
4372 if (tmpdist .eq. 0d0)
then
4373 radlon(1:na, ir, ie) = lon0
4374 radlat(1:na, ir, ie) = lat0
4376 cdist = cos(tmpdist)
4377 sdist = sin(tmpdist)
4378 sinll1_cdist = sinll1 * cdist
4379 cosll1_sdist = cosll1 * sdist
4380 cosll1_cdist = cosll1 * cdist
4381 sinll1_sdist = sinll1 * sdist
4383 radlat(ia, ir, ie) = asin(sinll1_cdist + cosll1_sdist * cos_azim(ia)) * r2d
4384 radlon(ia, ir, ie) = lon0 + atan2(sdist * sin_azim(ia), cosll1_cdist - sinll1_sdist * cos_azim(ia)) * r2d
4391 end subroutine radar_georeference
4396 subroutine radar_superobing(na, nr, ne, radlon, radlat, radz, ze, vr, & ! input spherical
4397 & qcflag, attenuation, & ! input spherical 2
4398 & nlon, nlat, nlev, lon, lat, z, dlon, dlat, dz, & ! input cartesian
4399 & missing, input_is_dbz, & ! input param
4400 & lon0, lat0, & ! input param
4401 & nobs_sp, grid_index, & ! output array info
4402 & grid_ref, grid_lon_ref, grid_lat_ref, grid_z_ref, grid_count_ref, & ! output ze
4403 & grid_vr, grid_lon_vr, grid_lat_vr, grid_z_vr, grid_count_vr )
4408 integer,
intent(in) :: na, nr, ne
4409 real(RP),
intent(in),
dimension(na, nr, ne) :: radlon, radlat, radz
4410 real(RP),
intent(in) :: ze(na, nr, ne), vr(na, nr, ne)
4411 real(RP),
intent(in) :: qcflag(na, nr, ne), attenuation(na, nr, ne)
4412 integer,
intent(in) :: nlon, nlat, nlev
4413 real(RP),
intent(in) :: lon(nlon), lat(nlat), z(nlev)
4414 real(RP),
intent(in) :: dlon, dlat, dz, missing
4415 logical,
intent(in) :: input_is_dbz
4416 real(RP),
intent(in) :: lon0, lat0
4418 integer(8),
intent(out) :: nobs_sp
4419 integer(8),
allocatable,
intent(out) :: grid_index(:), grid_count_ref(:), grid_count_vr(:)
4420 REAL(RP),
allocatable,
intent(out) :: grid_ref(:), grid_vr(:)
4421 REAL(RP),
allocatable,
intent(out) :: grid_lon_ref(:), grid_lat_ref(:), grid_z_ref(:)
4422 REAL(RP),
allocatable,
intent(out) :: grid_lon_vr(:), grid_lat_vr(:), grid_z_vr(:)
4424 REAL(RP),
ALLOCATABLE :: grid_w_vr(:)
4425 REAL(RP),
ALLOCATABLE :: grid_meanvr(:), grid_stdvr(:)
4427 integer(8),
allocatable :: packed_grid(:), pack2uniq(:), nobs_each_elev(:)
4428 real(RP),
allocatable :: packed_data(:, :)
4429 logical,
allocatable :: packed_attn(:)
4431 REAL(RP) :: count_inv, tmpstdvr, tmpmeanvr
4432 integer(8) :: idx, jdx, nobs, sidx(ne), eidx(ne)
4456 allocate(nobs_each_elev(ne))
4479 allocate(packed_grid(na * nr * ne), &
4480 & packed_data(5, na * nr * ne), &
4481 & packed_attn(na * nr * ne))
4482 call qc_indexing_and_packing( &
4483 & na, nr, ne, ze, vr, radlon, radlat, radz, &
4484 & qcflag, input_is_dbz, attenuation, &
4485 & nlon, nlat, nlev, lon, lat, z, dlon, dlat, dz, &
4488 & nobs_each_elev, packed_grid, packed_data, packed_attn)
4490 nobs = sum(nobs_each_elev)
4493 sidx(i + 1) = sidx(i) + nobs_each_elev(i)
4495 eidx = sidx + nobs_each_elev - 1
4496 deallocate(nobs_each_elev)
4521 allocate(grid_index(nobs))
4523 grid_index(idx) = packed_grid(idx)
4529 call merge_sort_parallel(nobs, grid_index)
4536 allocate(pack2uniq(nobs))
4542 pack2uniq(idx) = binary_search_i8(nobs_sp, grid_index, packed_grid(idx))
4554 allocate(grid_ref(nobs_sp), grid_vr(nobs_sp))
4555 allocate(grid_count_ref(nobs_sp), grid_count_vr(nobs_sp))
4556 allocate(grid_lon_ref(nobs_sp), grid_lat_ref(nobs_sp), grid_z_ref(nobs_sp))
4557 allocate(grid_lon_vr(nobs_sp) , grid_lat_vr(nobs_sp) , grid_z_vr(nobs_sp))
4558 allocate(grid_w_vr(nobs_sp))
4559 allocate(grid_meanvr(nobs_sp), grid_stdvr(nobs_sp))
4565 grid_count_ref(idx) = 0
4566 grid_count_vr(idx) = 0
4567 grid_ref(idx) = 0.0d0
4568 grid_vr(idx) = 0.0d0
4569 grid_w_vr(idx) = 0.0d0
4570 grid_lon_ref(idx) = 0.0d0
4571 grid_lat_ref(idx) = 0.0d0
4572 grid_z_ref(idx) = 0.0d0
4573 grid_lon_vr(idx) = 0.0d0
4574 grid_lat_vr(idx) = 0.0d0
4575 grid_z_vr(idx) = 0.0d0
4576 grid_meanvr(idx) = 0.0d0
4577 grid_stdvr(idx) = 0.0d0
4581 idx = pack2uniq(jdx)
4584 if(packed_attn(jdx))
then
4585 grid_ref(idx) = grid_ref(idx) + packed_data(1, jdx)
4586 grid_count_ref(idx) = grid_count_ref(idx) + 1
4587 grid_lon_ref(idx) = grid_lon_ref(idx) + packed_data(3, jdx)
4588 grid_lat_ref(idx) = grid_lat_ref(idx) + packed_data(4, jdx)
4589 grid_z_ref(idx) = grid_z_ref(idx) + packed_data(5, jdx)
4595 IF(packed_data(2, jdx) .GT. missing)
THEN
4596 grid_w_vr(idx) = grid_w_vr(idx) + packed_data(1, jdx)
4597 grid_count_vr(idx) = grid_count_vr(idx) + 1
4598 grid_meanvr(idx) = grid_meanvr(idx) + packed_data(2, jdx)
4599 grid_stdvr(idx) = grid_stdvr(idx) + packed_data(2, jdx) ** 2
4600 grid_vr(idx) = grid_vr(idx) + packed_data(2, jdx) * packed_data(1, jdx)
4601 grid_lon_vr(idx) = grid_lon_vr(idx) + packed_data(3, jdx) * packed_data(1, jdx)
4602 grid_lat_vr(idx) = grid_lat_vr(idx) + packed_data(4, jdx) * packed_data(1, jdx)
4603 grid_z_vr(idx) = grid_z_vr(idx) + packed_data(5, jdx) * packed_data(1, jdx)
4609 IF(grid_count_ref(idx) > 0)
THEN
4610 count_inv = 1.0d0 / dble(grid_count_ref(idx))
4611 grid_ref(idx) = grid_ref(idx) * count_inv
4612 grid_lon_ref(idx) = grid_lon_ref(idx) * count_inv
4613 grid_lat_ref(idx) = grid_lat_ref(idx) * count_inv
4614 grid_z_ref(idx) = grid_z_ref(idx) * count_inv
4617 IF(grid_count_vr(idx) > 0)
THEN
4618 count_inv = 1.0d0 / grid_w_vr(idx)
4619 grid_vr(idx) = grid_vr(idx) * count_inv
4620 grid_lon_vr(idx) = grid_lon_vr(idx) * count_inv
4621 grid_lat_vr(idx) = grid_lat_vr(idx) * count_inv
4622 grid_z_vr(idx) = grid_z_vr(idx) * count_inv
4629 IF( radar_use_vr_std )
THEN
4630 count_inv = 1.0d0 / dble(grid_count_vr(idx))
4631 tmpmeanvr = grid_meanvr(idx) * count_inv
4632 tmpstdvr = grid_stdvr(idx) * count_inv
4633 tmpstdvr = sqrt(tmpstdvr - (tmpmeanvr ** 2))
4634 IF(tmpstdvr > vr_std_threshold) grid_count_vr(idx) = 0
4645 subroutine define_grid(lon0, lat0, nr, rrange, maxrange, maxz, dx, dy, dz, & ! input
4646 & dlon, dlat, nlon, nlat, nlev, lon, lat, z)
4652 real(RP),
intent(in) :: lon0, lat0
4653 integer,
intent(in) :: nr
4654 real(RP),
intent(in) :: rrange(nr)
4655 real(RP),
intent(in) :: maxrange, maxz, dx, dy, dz
4656 real(RP),
intent(out) :: dlon, dlat
4657 integer,
intent(out) :: nlon, nlat, nlev
4658 real(RP),
allocatable,
intent(out) :: lon(:), lat(:), z(:)
4659 real(RP) :: tmpmaxrange
4665 dlon = r2d * dx / (cos(lat0 * d2r) * radius)
4666 dlat = r2d * dx / radius
4669 tmpmaxrange = maxrange
4670 tmpmaxrange = min(tmpmaxrange, rrange(nr))
4671 tmpmaxrange = 2.0 * tmpmaxrange
4672 nlon = ceiling(tmpmaxrange / dx)
4673 nlat = ceiling(tmpmaxrange / dy)
4674 nlev = ceiling(maxz / dz)
4676 allocate(lon(nlon), lat(nlat), z(nlev))
4679 IF( mod(nlon,2) == 0 ) nlon = nlon + 1
4680 IF( mod(nlat,2) == 0 ) nlat = nlat + 1
4683 lon(i) = lon0 + dlon * (-1.0 - (nlon - 1.0) / 2.0 + i)
4687 lat(j) = lat0 + dlat * (-1.0 - (nlat - 1.0) / 2.0 + j)
4694 end subroutine define_grid
4702 subroutine calc_ref_vr(qv,qc,qr,qci,qs,qg,u,v,w,t,p,az,elev,ref,vr)
4710 real(RP),
intent(in) :: qv
4711 real(RP),
intent(in) :: qc,qr
4712 real(RP),
intent(in) :: qci,qs,qg
4713 real(RP),
intent(in) :: t,p
4714 real(RP),
intent(in) :: u,v,w
4715 real(RP),
intent(inout) :: ref
4717 real(RP),
intent(in) :: az
4718 real(RP),
intent(in) :: elev
4719 real(RP),
intent(inout) :: vr
4720 real(RP) :: qms , qmg
4722 real(RP) :: zr , zs , zg
4723 real(RP) :: wr , ws , wg
4725 real(RP) :: nor, nos, nog
4726 real(RP) :: ror, ros, rog , roi
4727 real(RP) :: a,b,c,d,Cd
4728 real(RP) :: cr_t08,cs_t08,cg_t08,dr_t08,ds_t08,dg_t08
4729 real(RP) :: cf, pip , roo
4730 real(RP) :: ki2 , kr2
4731 real(RP) :: lr , ls , lg
4732 real(RP) :: tmp_factor , rofactor
4734 real(RP) :: Fs, Fg , zms , zmg , fws , fwg
4735 real(RP) :: qrp , qsp , qgp
4738 real(RP) ,
parameter :: qeps = 1.0d-20
4740 real(RP) ,
parameter :: as_RS14 = 6.9d-2
4741 real(RP) :: tc , MOMs_0bs
4761 ro = p / ( rdry * t)
4766 if (method_ref_calc == 1)
then
4789 IF( qt .GT. qeps )
THEN
4790 ref = cf * ( ( ro * qt )**1.75 )
4791 ref = ref / ( pip * ( nor ** 0.75 ) * ( ror ** 1.75 ) )
4799 IF ( qt .GT. qeps )
THEN
4801 wt = 5.40d0 * a * ( qt ** 0.125 )
4807 else if (method_ref_calc == 2)
then
4827 IF( qr .GT. qeps )
THEN
4828 zr= cf * ( ( ro * qr )**1.75 )
4829 zr= zr / ( pip * ( nor ** 0.75 ) * ( ror ** 1.75 ) )
4832 IF( qs .GT. qeps )
THEN
4833 IF ( t <= 273.16 )
THEN
4834 zs = cf * ki2 * ( ros ** 0.25 ) * ( ( ro * qs ) ** 1.75 )
4835 zs = zs / ( pip * kr2 * ( nos ** 0.75 ) * ( roi ** 2 ) )
4840 zs = cf * ( ( ro * qs ) ** 1.75 )
4841 zs = zs / ( pip * ( nos ** 0.75 ) * ( roi ** 1.75 ) )
4846 IF( qg .GT. qeps )
THEN
4847 zg= ( cf / ( pip * ( nog ** 0.75) * ( rog ** 1.75 ) ) ) ** 0.95
4848 zg= zg * ( ( ro * qg ) ** 1.6625 )
4855 IF( ref > 0.0d0 )
THEN
4874 rofactor= ( roo / ro ) ** 0.25
4876 CALL com_gamma( 4.0_rp + b , tmp_factor )
4877 lr= ( pi * ror * nor / ( ro * qr ) ) ** 0.25
4878 wr= a * tmp_factor / ( 6.0d0 * ( lr ** b ) )
4879 wr= 1.0d-2*wr * rofactor
4885 CALL com_gamma( 4.0_rp + d , tmp_factor )
4886 ls= ( pi * ros * nos / ( ro * qs ) ) ** 0.25
4887 ws= c * tmp_factor / ( 6.0d0 * ( ls ** d ) )
4888 ws= 1.0d-2*ws * rofactor
4895 lg= ( pi * rog * nog / ( ro * qg ) ) ** 0.25
4896 wg= tmp_factor * ( ( ( 4.0d0 * grav * 100.0d0 * rog )/( 3.0d0 * cd * ro ) ) ** 0.5 )
4897 wg= 1.0d-2*wg / ( 6.0d0 * ( lg ** 0.5 ) )
4903 wt = ( wr * zr + ws * zs + wg * zg )/ ( zr + zs + zg )
4911 else if (method_ref_calc == 3)
then
4926 IF( qr .GT. qeps .AND. qg .GT. qeps)
THEN
4927 fg=maxf * ( min( qr/qg , qg/qr ) )**(1.0d0/3.0d0)
4928 fwg= qr / ( qr + qg )
4930 IF( qr .GT. qeps .AND. qs .GT. qeps)
THEN
4931 fs=maxf * ( min( qr/qs , qs/qr ) )**(1.0d0/3.0d0)
4932 fws= qr / ( qr + qs )
4935 if ( .not. use_method3_ref_melt )
then
4943 qrp=(1.0d0-fs-fg)*qr
4957 IF( qrp .GT. qeps)
THEN
4958 zr= 2.53d4 * ( ro * qrp * 1.0d3 )**1.84
4960 IF( qsp .GT. qeps)
THEN
4961 zs= 3.48d3 * ( ro * qsp * 1.0d3 )**1.66
4963 IF( qgp .GT. qeps)
THEN
4964 zg= 5.54d3 * ( ro * qgp * 1.0d3 )**1.70
4966 IF( qms .GT. qeps )
THEN
4967 zms=( 0.00491 + 5.75*fws - 5.588*(fws**2) )*1.0d5
4968 zms= zms * ( ro * qms * 1.0d3 )**( 1.67 - 0.202*fws + 0.398*(fws**2) )
4971 IF( qmg .GT. qeps )
THEN
4972 zmg=( 0.0358 + 5.27*fwg -9.51*(fwg**2) + 4.68 *(fwg**3) )*1.0d5
4973 zmg= zmg * ( ro * qmg * 1.0d3 )**( 1.70 + 0.020*fwg + 0.287 * (fwg**2) - 0.186*(fwg**3) )
4976 ref = zr + zg + zs + zms + zmg
4985 IF( ref > 0.0d0 )
THEN
5005 rofactor= ( roo / ro ) ** 0.5
5007 IF ( qr .GT. qeps )
THEN
5008 lr= ( pi * ror * nor / ( ro * qr ) ) ** 0.25
5009 CALL com_gamma( 4.0_rp + dr_t08 , tmp_factor )
5010 wr= cr_t08 * tmp_factor / ( 6.0d0 * ( ( lr * 1.0e2 ) ** dr_t08 ) )
5016 IF( qs .GT. qeps )
THEN
5017 IF ( use_t08_rs2014 )
then
5018 tc=min(-0.1_rp, t-273.15_rp)
5019 moms_0bs = exp(log(10.0_rp) * loga_(tc,2.0_rp) + log(ro * qs / as_rs14) * b_(tc,2.0_rp) )
5020 ws = cs_t08 * rofactor * exp(log(10.0_rp) * loga_(tc,2.25_rp) + log(ro * qs / as_rs14) * b_(tc,2.25_rp) ) / moms_0bs
5022 ls= ( pi * ros * nos / ( ro * qs ) ) ** 0.25
5023 CALL com_gamma( 4.0_rp + ds_t08 , tmp_factor )
5024 ws= cs_t08 * tmp_factor / ( 6.0d0 * ( ( ls * 1.0e2 ) ** ds_t08 ) )
5032 IF ( qg .GT. qeps )
THEN
5033 lg= ( pi * rog * nog / ( ro * qg ) ) ** 0.25
5034 CALL com_gamma( 4.0_rp + dg_t08 , tmp_factor )
5035 wg = tmp_factor * ( ( ( 4.0d0 * grav * rog )/( 3.0d0 * cd * roo ) ) ** 0.5 )
5036 wg = wg / ( 6.0d0 * ( ( lg * 1.0e2 ) ** dg_t08 ) )
5045 wt = ( wr * zr + ws * zs + ws * zms + wg * zg + wg * zmg ) / ( zr + zs + zg + zms + zmg )
5055 WRITE(6,*)
'[Error] Not recognized method for radar reflectivity and wind computation'
5063 vr = u * cos(elev*d2r) * sin(az*d2r)
5064 vr = vr + v * cos(elev*d2r) * cos(az*d2r)
5065 IF( use_terminal_velocity )
THEN
5066 vr = vr + (w - wt)*sin(elev*d2r)
5068 vr = vr + (w)*sin(elev*d2r)
5072 if (.not.(abs(vr) .lt. 1.0e20))
then
5073 write(*,*)
'***ERROR*** vr is NaN'
5077 write(*,*) zr,zg,zs,zms,zmg
5079 WRITE(6,*) elev , az , d2r
5081 elseif (.not. (abs(ref) .lt. 1.0e20))
then
5082 write(*,*)
'***ERROR*** ref is NaN'
5084 write(*,*) zr, zg, zs, zms, zmg
5086 elseif (ref < 0.0)
then
5087 write(*,*)
'***ERROR*** ref is negative'
5089 write(*,*) zr, zg, zs, zms, zmg
5112 real(RP) :: z , r , gr
5113 integer :: m1, k , m
5116 if (x.eq.int(x))
then
5117 if (x.gt.0.0d0)
then
5128 if (abs(x).gt.1.0_rp)
then
5141 0.5772156649015329e0, &
5142 -0.6558780715202538e0, &
5143 -0.420026350340952e-1, &
5144 0.1665386113822915e0, &
5145 -0.421977345555443e-1, &
5146 -0.96219715278770e-2, &
5147 0.72189432466630e-2, &
5148 -0.11651675918591e-2, &
5149 -0.2152416741149e-3, &
5150 0.1280502823882e-3, &
5151 -0.201348547807e-4, &
5152 -0.12504934821e-5, &
5171 if (abs(x).gt.1.0_rp)
then
5173 if (x.lt.0.0_rp) ga=-pi/(x*ga*sin(pi*x))
5180 function loga_(tems, nm)
5184 real(RP),
parameter :: coef_a01 = 5.065339_rp
5185 real(RP),
parameter :: coef_a02 = -0.062659_rp
5186 real(RP),
parameter :: coef_a03 = -3.032362_rp
5187 real(RP),
parameter :: coef_a04 = 0.029469_rp
5188 real(RP),
parameter :: coef_a05 = -0.000285_rp
5189 real(RP),
parameter :: coef_a06 = 0.31255_rp
5190 real(RP),
parameter :: coef_a07 = 0.000204_rp
5191 real(RP),
parameter :: coef_a08 = 0.003199_rp
5192 real(RP),
parameter :: coef_a09 = 0.0_rp
5193 real(RP),
parameter :: coef_a10 = -0.015952_rp
5194 real(RP) :: coef_at(4)
5196 coef_at(1) = coef_a01 + tems * ( coef_a02 + tems * ( coef_a05 + tems * coef_a09 ) )
5197 coef_at(2) = coef_a03 + tems * ( coef_a04 + tems * coef_a07 )
5198 coef_at(3) = coef_a06 + tems * coef_a08
5199 coef_at(4) = coef_a10
5200 loga_ = coef_at(1) + nm * ( coef_at(2) + nm * ( coef_at(3) + nm * coef_at(4) ) )
5204 function b_(tems, nm)
5208 real(RP),
parameter :: coef_b01 = 0.476221_rp
5209 real(RP),
parameter :: coef_b02 = -0.015896_rp
5210 real(RP),
parameter :: coef_b03 = 0.165977_rp
5211 real(RP),
parameter :: coef_b04 = 0.007468_rp
5212 real(RP),
parameter :: coef_b05 = -0.000141_rp
5213 real(RP),
parameter :: coef_b06 = 0.060366_rp
5214 real(RP),
parameter :: coef_b07 = 0.000079_rp
5215 real(RP),
parameter :: coef_b08 = 0.000594_rp
5216 real(RP),
parameter :: coef_b09 = 0.0_rp
5217 real(RP),
parameter :: coef_b10 = -0.003577_rp
5218 real(RP) :: coef_bt(4)
5220 coef_bt(1) = coef_b01 + tems * ( coef_b02 + tems * ( coef_b05 + tems * coef_b09 ) )
5221 coef_bt(2) = coef_b03 + tems * ( coef_b04 + tems * coef_b07 )
5222 coef_bt(3) = coef_b06 + tems * coef_b08
5223 coef_bt(4) = coef_b10
5224 b_ = coef_bt(1) + nm * ( coef_bt(2) + nm * ( coef_bt(3) + nm * coef_bt(4) ) )
5228 end subroutine calc_ref_vr
5251 subroutine obs_local(obs, ri, rj, rlev, rz, nvar, hdxf, rdiag, rloc, dep, nobsl, depd, nobsl_t, cutd_t, srch_q0)
5257 type(
obs_info),
intent(in) :: obs(:)
5258 real(RP),
intent(in) :: ri, rj, rlev, rz
5259 integer,
intent(in) :: nvar
5260 real(RP),
intent(out) :: hdxf(nobstotal,nmem)
5261 real(RP),
intent(out) :: rdiag(nobstotal)
5262 real(RP),
intent(out) :: rloc(nobstotal)
5263 real(RP),
intent(out) :: dep(nobstotal)
5264 integer,
intent(out) :: nobsl
5265 real(RP),
intent(out),
optional :: depd(nobstotal)
5266 integer,
intent(out),
optional :: nobsl_t(nid_obs,nobtype)
5267 real(RP),
intent(out),
optional :: cutd_t(nid_obs,nobtype)
5268 integer,
intent(inout),
optional :: srch_q0(nctype)
5270 integer,
allocatable :: nobs_use(:)
5271 integer,
allocatable :: nobs_use2(:)
5272 real(RP),
allocatable :: dist_tmp(:)
5273 real(RP),
allocatable :: rloc_tmp(:)
5274 real(RP),
allocatable :: rdiag_tmp(:)
5276 real(RP) :: nrloc, nrdiag
5277 real(RP) :: ndist_dummy
5278 integer :: iob, ityp, ielm
5279 integer :: imin, imax, jmin, jmax
5280 integer :: ic, ic2, icm
5281 integer :: n, nn, nn_prev
5282 integer :: nobsl_prev, nobsl_incr
5283 integer :: nobsl_max_master
5284 integer :: ielm_u_master, ityp_master
5288 integer :: nn_steps(n_merge_max+1)
5289 logical :: reach_cutoff
5290 integer :: imin_cutoff(n_merge_max), imax_cutoff(n_merge_max)
5291 integer :: jmin_cutoff(n_merge_max), jmax_cutoff(n_merge_max)
5292 real(RP) :: search_incr(n_merge_max)
5293 real(RP) :: search_incr_i(n_merge_max), search_incr_j(n_merge_max)
5294 real(RP) :: dist_cutoff_fac, dist_cutoff_fac_square
5297 integer :: nobsl_cm(n_merge_max)
5304 if (
present(nobsl_t))
then
5307 if (
present(cutd_t))
then
5309 if (max_nobs_per_grid_criterion == 1)
then
5311 cutd_t(elm_u_ctype(ic),typ_ctype(ic)) = hori_loc_ctype(ic) * dist_zero_fac
5316 if (nobstotal == 0)
then
5320 if (maxval(max_nobs_per_grid(:)) > 0)
then
5321 allocate (nobs_use(nobstotal))
5322 allocate (nobs_use2(nobstotal))
5323 allocate (rloc_tmp(nobstotal))
5324 allocate (rdiag_tmp(nobstotal))
5325 if (max_nobs_per_grid_criterion == 1)
then
5326 allocate (dist_tmp(nobstotal))
5329 rloc_tmp(:) = -1.0d6
5332 allocate (nobs_use(maxnobs_per_ctype))
5342 if (n_merge(ic) == 0)
then
5343 if (
present(cutd_t))
then
5344 cutd_t(elm_u_ctype(ic),typ_ctype(ic)) = 0.0d0
5349 nobsl_max_master = max_nobs_per_grid(typ_ctype(ic))
5350 ielm_u_master = elm_u_ctype(ic)
5351 ityp_master = typ_ctype(ic)
5353 if (nobsl_max_master <= 0)
then
5361 do icm = 1, n_merge(ic)
5362 ic2 = ic_merge(icm,ic)
5363 ielm = elm_ctype(ic2)
5364 ityp = typ_ctype(ic2)
5366 if (obsgrd(ic2)%tot_ext > 0)
then
5368 call obs_local_range(ic2, ri, rj, imin, imax, jmin, jmax)
5369 call obs_choose_ext(ic2, imin, imax, jmin, jmax, nn, nobs_use)
5374 call obs_local_cal(obs(:),ri, rj, rlev, rz, nvar, iob, ic2, ndist_dummy, nrloc, nrdiag)
5375 if (nrloc == 0.0d0) cycle
5378 hdxf(nobsl,:) = obsda_sort%ensval(1:nmem,iob)
5379 rdiag(nobsl) = nrdiag
5381 dep(nobsl) = obsda_sort%val(iob)
5382 if (
present(depd))
then
5383 depd(nobsl) = obsda_sort%ensval(mmdetobs,iob)
5388 if (
present(nobsl_t))
then
5389 nobsl_t(elm_u_ctype(ic2),ityp) = nobsl - nobsl_prev
5394 else if (max_nobs_per_grid_criterion == 1)
then
5402 do icm = 1, n_merge(ic)
5403 ic2 = ic_merge(icm,ic)
5405 if (obsgrd(ic2)%tot_ext > 0)
then
5406 call obs_local_range(ic2, ri, rj, imin, imax, jmin, jmax)
5407 call obs_choose_ext(ic2, imin, imax, jmin, jmax, nn, nobs_use)
5415 search_incr(1) = hori_loc_ctype(ic) * dist_zero_fac / real(
n_search_incr, rp)
5416 search_incr(1) = max(search_incr(1), obsgrd(ic)%grdspc_i, obsgrd(ic)%grdspc_j)
5418 do icm = 1, n_merge(ic)
5419 ic2 = ic_merge(icm,ic)
5423 search_incr(icm) = search_incr(1) / hori_loc_ctype(ic) * hori_loc_ctype(ic2)
5425 search_incr_i(icm) = search_incr(icm) / dx
5426 search_incr_j(icm) = search_incr(icm) / dy
5427 call obs_local_range(ic2, ri, rj, imin_cutoff(icm), imax_cutoff(icm), jmin_cutoff(icm), jmax_cutoff(icm))
5431 if (
present(srch_q0))
then
5441 reach_cutoff = .true.
5443 do icm = 1, n_merge(ic)
5444 ic2 = ic_merge(icm,ic)
5447 if (obsgrd(ic2)%tot_ext > 0)
then
5448 call ij_obsgrd_ext(ic2, ri-search_incr_i(icm)*q, rj-search_incr_j(icm)*q, imin, jmin)
5449 call ij_obsgrd_ext(ic2, ri+search_incr_i(icm)*q, rj+search_incr_j(icm)*q, imax, jmax)
5451 if (imin <= imin_cutoff(icm) .and. imax >= imax_cutoff(icm) .and. &
5452 jmin <= jmin_cutoff(icm) .and. jmax >= jmax_cutoff(icm))
then
5453 imin = imin_cutoff(icm)
5454 imax = imax_cutoff(icm)
5455 jmin = jmin_cutoff(icm)
5456 jmax = jmax_cutoff(icm)
5458 reach_cutoff = .false.
5461 call obs_choose_ext(ic2, imin, imax, jmin, jmax, nn, nobs_use)
5465 nn_steps(n_merge(ic)+1) = nn
5467 if ((.not. reach_cutoff) .and. nn < nobsl_max_master) cycle
5468 if (reach_cutoff)
then
5477 dist_cutoff_fac = search_incr(1) * q / hori_loc_ctype(ic)
5478 dist_cutoff_fac_square = dist_cutoff_fac * dist_cutoff_fac
5480 do icm = 1, n_merge(ic)
5481 ic2 = ic_merge(icm,ic)
5483 if (nn_steps(icm+1) > nn_steps(icm))
then
5484 do n = nn_steps(icm)+1, nn_steps(icm+1)
5487 if (rloc_tmp(iob) == 0.0d0) cycle
5489 if (rloc_tmp(iob) < 0.0d0)
then
5490 call obs_local_cal(obs(:),ri, rj, rlev, rz, nvar, iob, ic2, dist_tmp(iob), rloc_tmp(iob), rdiag_tmp(iob))
5491 if (rloc_tmp(iob) == 0.0d0) cycle
5494 if (.not. reach_cutoff)
then
5495 if (dist_tmp(iob) > dist_cutoff_fac_square) cycle
5498 nobsl_incr = nobsl_incr + 1
5499 nobs_use2(nobsl_incr) = iob
5504 if (nobsl_incr >= nobsl_max_master) loop = .false.
5507 if (
present(srch_q0))
then
5508 if (q == srch_q0(ic) .and. nobsl_incr > nobsl_max_master * 3)
then
5510 else if (q > srch_q0(ic))
then
5515 if (nobsl_incr == 0) cycle
5517 if (nobsl_incr > nobsl_max_master)
then
5519 nobsl_incr = nobsl_max_master
5522 do n = 1, nobsl_incr
5525 hdxf(nobsl,:) = obsda_sort%ensval(1:nmem,iob)
5526 rdiag(nobsl) = rdiag_tmp(iob)
5527 rloc(nobsl) = rloc_tmp(iob)
5528 dep(nobsl) = obsda_sort%val(iob)
5529 if (
present(depd))
then
5530 depd(nobsl) = obsda_sort%ensval(mmdetobs,iob)
5534 if (
present(nobsl_t))
then
5535 nobsl_t(ielm_u_master,ityp_master) = nobsl_incr
5537 if (
present(cutd_t))
then
5538 if (nobsl_incr == nobsl_max_master)
then
5539 cutd_t(ielm_u_master,ityp_master) = hori_loc_ctype(ic) * sqrt(dist_tmp(nobs_use2(nobsl_incr)))
5553 do icm = 1, n_merge(ic)
5554 ic2 = ic_merge(icm,ic)
5556 if (obsgrd(ic2)%tot_ext > 0)
then
5558 call obs_local_range(ic2, ri, rj, imin, imax, jmin, jmax)
5559 call obs_choose_ext(ic2, imin, imax, jmin, jmax, nn, nobs_use)
5561 do n = nn_prev+1, nn
5564 call obs_local_cal(obs(:),ri, rj, rlev, rz, nvar, iob, ic2, ndist_dummy, rloc_tmp(iob), rdiag_tmp(iob))
5565 if (rloc_tmp(iob) == 0.0d0) cycle
5567 nobsl_incr = nobsl_incr + 1
5568 nobs_use2(nobsl_incr) = iob
5573 if (nobsl_incr == 0) cycle
5575 if (nobsl_incr > nobsl_max_master)
then
5576 if (max_nobs_per_grid_criterion == 2)
then
5578 else if (max_nobs_per_grid_criterion == 3)
then
5581 write (6,
'(A,I6)')
"[Error] Unsupported 'MAX_NOBS_PER_GRID_CRITERION':", max_nobs_per_grid_criterion
5584 nobsl_incr = nobsl_max_master
5587 do n = 1, nobsl_incr
5590 hdxf(nobsl,:) = obsda_sort%ensval(1:nmem,iob)
5591 rdiag(nobsl) = rdiag_tmp(iob)
5592 rloc(nobsl) = rloc_tmp(iob)
5593 dep(nobsl) = obsda_sort%val(iob)
5594 if (
present(depd))
then
5595 depd(nobsl) = obsda_sort%ensval(mmdetobs,iob)
5599 if (
present(nobsl_t))
then
5600 nobsl_t(ielm_u_master,ityp_master) = nobsl_incr
5602 if (
present(cutd_t))
then
5603 if (nobsl_incr == nobsl_max_master)
then
5604 if (max_nobs_per_grid_criterion == 2)
then
5605 cutd_t(ielm_u_master,ityp_master) = rloc_tmp(nobs_use2(nobsl_incr))
5606 else if (max_nobs_per_grid_criterion == 3)
then
5607 cutd_t(ielm_u_master,ityp_master) = rdiag_tmp(nobs_use2(nobsl_incr))
5621 deallocate (nobs_use)
5622 if (maxval(max_nobs_per_grid(:)) > 0)
then
5623 deallocate (nobs_use2)
5624 deallocate (rloc_tmp)
5625 deallocate (rdiag_tmp)
5626 if (max_nobs_per_grid_criterion == 1)
then
5627 deallocate (dist_tmp)
5632 end subroutine obs_local
5638 subroutine obs_local_range(ctype, ri, rj, imin, imax, jmin, jmax)
5640 integer,
intent(in) :: ctype
5641 real(RP),
intent(in) :: ri, rj
5642 integer,
intent(out) :: imin, imax, jmin, jmax
5644 real(RP) :: dist_zero_i, dist_zero_j
5648 dist_zero_i = 0.0_rp
5650 dist_zero_i = hori_loc_ctype(ctype) * dist_zero_fac / dx
5652 dist_zero_j = hori_loc_ctype(ctype) * dist_zero_fac / dy
5653 call ij_obsgrd_ext(ctype, ri - dist_zero_i, rj - dist_zero_j, imin, jmin)
5654 call ij_obsgrd_ext(ctype, ri + dist_zero_i, rj + dist_zero_j, imax, jmax)
5657 end subroutine obs_local_range
5662 subroutine obs_local_cal(obs, ri, rj, rlev, rz, nvar, iob, ic, ndist, nrloc, nrdiag)
5664 type(
obs_info),
intent(in) :: obs(:)
5665 real(RP),
intent(in) :: ri, rj, rlev, rz
5666 integer,
intent(in) :: nvar
5667 integer,
intent(in) :: iob
5668 integer,
intent(in) :: ic
5669 real(RP),
intent(out) :: ndist
5670 real(RP),
intent(out) :: nrloc
5671 real(RP),
intent(out) :: nrdiag
5677 real(RP) :: rdx, rdy, rdz
5678 real(RP) :: nd_h, nd_v
5680 integer :: di, dj, dk
5686 obelm = elm_ctype(ic)
5687 obtyp = typ_ctype(ic)
5688 obset = obsda_sort%set(iob)
5689 obidx = obsda_sort%idx(iob)
5694 nrloc = var_local(nvar,uid_obs_varlocal(obelm))
5697 if (nrloc < tiny(var_local))
then
5705 if (vert_loc_ctype(ic) == 0.0d0)
then
5707 else if (obelm == id_ps_obs)
then
5708 nd_v = abs(log(obs(obset)%dat(obidx)) - log(rlev)) / vert_loc_ctype(ic)
5709 else if (obelm == id_rain_obs)
then
5710 nd_v = abs(log(vert_local_rain_base) - log(rlev)) / vert_loc_ctype(ic)
5711 else if (obtyp == 22)
then
5712 nd_v = abs(obs(obset)%lev(obidx) - rz) / vert_loc_ctype(ic)
5714 nd_v = abs(log(obs(obset)%lev(obidx)) - log(rlev)) / vert_loc_ctype(ic)
5719 if (nd_v > dist_zero_fac)
then
5726 rdx = (ri - obs(obset)%ri(obidx)) * dx
5727 rdy = (rj - obs(obset)%rj(obidx)) * dy
5728 nd_h = sqrt(rdx*rdx + rdy*rdy) / hori_loc_ctype(ic)
5731 if (nd_h > dist_zero_fac)
then
5738 ndist = nd_h * nd_h + nd_v * nd_v
5741 if (ndist > dist_zero_fac_square)
then
5747 if ( obtyp == 22 .and. ( radar_thin_letkf_method > 0 ) )
then
5748 rdz = obs(obset)%lev(obidx) - rz
5751 select case( radar_thin_letkf_method )
5755 di = int( abs( rdx / radar_so_size_hori ) )
5756 dj = int( abs( rdy / radar_so_size_hori ) )
5757 dk = int( abs( obs(obset)%lev(obidx) - rz ) / radar_so_size_vert )
5759 if ( ( mod( di, radar_thin_letkf_hgrid ) /= 0 .or. &
5760 mod( dj, radar_thin_letkf_hgrid ) /= 0 .or. &
5761 mod( dk, radar_thin_letkf_vgrid ) /= 0 ) .and. &
5762 ( ( di >= radar_thin_letkf_hnear ) .or. &
5763 ( dj >= radar_thin_letkf_hnear ) .or. &
5764 ( dk >= radar_thin_letkf_vnear ) ) )
then
5772 di = nint( rdx / radar_so_size_hori )
5773 dj = nint( rdy / radar_so_size_hori )
5774 dk = nint( obs(obset)%lev(obidx) - rz ) / radar_so_size_vert
5776 if ( mod( di, radar_thin_letkf_hgrid ) /= 0 .or. &
5777 mod( dj, radar_thin_letkf_hgrid ) /= 0 .or. &
5778 mod( dk, radar_thin_letkf_vgrid ) /= 0 )
then
5794 nrloc = nrloc * exp(-0.5d0 * ndist)
5798 nrdiag = obs(obset)%err(obidx) * obs(obset)%err(obidx) / nrloc
5799 if ( radar_pqv .and. obelm == id_radar_ref_obs .and. obsda_sort%tm(iob) < 0.0d0 )
then
5800 nrdiag = obserr_pq**2 / nrloc
5804 end subroutine obs_local_cal
5819 subroutine obs_pest_etkf(hdxf, rdiag, rloc, dep, nobsl, depd)
5822 real(RP),
intent(out) :: hdxf(nobstotal,nmem)
5823 real(RP),
intent(out) :: rdiag(nobstotal)
5824 real(RP),
intent(out) :: rloc(nobstotal)
5825 real(RP),
intent(out) :: dep(nobstotal)
5826 integer,
intent(out) :: nobsl
5827 real(RP),
intent(out),
optional :: depd(nobstotal)
5832 real(RP) :: nrloc, nrdiag
5844 nn = obsda_sort%nobs
5848 obset = obsda_sort%set(n)
5849 obidx = obsda_sort%idx(n)
5850 nrdiag = obs(obset)%err(obidx) * obs(obset)%err(obidx) / nrloc
5851 if (nrloc == 0.0) cycle
5854 hdxf(nobsl,:) = obsda_sort%ensval(1:nmem,iob)
5855 rdiag(nobsl) = nrdiag
5857 dep(nobsl) = obsda_sort%val(iob)
5858 if (
present(depd))
then
5859 depd(nobsl) = obsda_sort%ensval(mmdetobs,iob)
5864 end subroutine obs_pest_etkf
5869 subroutine relax_beta(ri, rj, rz, beta)
5871 real(RP),
intent(in) :: ri, rj, rz
5872 real(RP),
intent(out) :: beta
5873 real(RP) :: dist_bdy
5879 if (radar_only .and. rz > radar_zmax + max(vert_local(22), vert_local_radar_vr) * dist_zero_fac)
then
5886 if (boundary_buffer_width > 0.0d0)
then
5887 dist_bdy = min(min(ri-xhalo, nlong+xhalo+1-ri) * dx, &
5888 min(rj-yhalo, nlatg+yhalo+1-rj) * dy) / boundary_buffer_width
5889 if (dist_bdy < 1.0d0)
then
5890 beta = max(dist_bdy, 0.0_rp)
5895 end subroutine relax_beta
5900 subroutine weight_rtpp(w, infl, wrlx)
5902 real(RP),
intent(in) :: w(nmem,nmem)
5903 real(RP),
intent(in) :: infl
5904 real(RP),
intent(out) :: wrlx(nmem,nmem)
5907 wrlx = (1.0d0 - relax_alpha) * w
5909 wrlx(m,m) = wrlx(m,m) + relax_alpha * sqrt(infl)
5913 end subroutine weight_rtpp
5918 subroutine weight_rtps(w, pa, xb, infl, wrlx, infl_out)
5920 real(RP),
intent(in) :: w(nmem,nmem)
5921 real(RP),
intent(in) :: pa(nmem,nmem)
5922 real(RP),
intent(in) :: xb(nmem)
5923 real(RP),
intent(in) :: infl
5924 real(RP),
intent(out) :: wrlx(nmem,nmem)
5925 real(RP),
intent(out) :: infl_out
5926 real(RP) :: var_g, var_a
5932 var_g = var_g + xb(m) * xb(m)
5934 var_a = var_a + xb(k) * pa(k,m) * xb(m)
5937 if (var_g > 0.0d0 .and. var_a > 0.0d0)
then
5938 infl_out = relax_alpha_spread * sqrt(var_g * infl / (var_a * real(nmem-1,rp))) &
5939 - relax_alpha_spread + 1.0d0
5947 end subroutine weight_rtps
5952 subroutine weight_rtps_const(w, pa, xb, wrlx, infl_out)
5954 real(RP),
intent(in) :: w(nmem,nmem)
5955 real(RP),
intent(in) :: pa(nmem,nmem)
5956 real(RP),
intent(in) :: xb(nmem)
5957 real(RP),
intent(out) :: wrlx(nmem,nmem)
5958 real(RP),
intent(out) :: infl_out
5959 real(RP) :: var_g, var_a
5962 real(RP),
parameter :: RTPS_const = 1.0d0
5967 var_g = var_g + xb(m) * xb(m)
5969 var_a = var_a + xb(k) * pa(k,m) * xb(m)
5972 if (var_g > 0.0d0 .and. var_a > 0.0d0)
then
5973 infl_out = rtps_const * sqrt(var_g * 1.0d0 / (var_a * real(nmem-1,kind=rp))) &
5974 - rtps_const + 1.0d0
5982 end subroutine weight_rtps_const
5987 subroutine com_distll_1(alon,alat,blon,blat,dist)
5992 real(RP),
intent(in) :: alon
5993 real(RP),
intent(in) :: alat
5994 real(RP),
intent(in) :: blon
5995 real(RP),
intent(in) :: blat
5996 real(RP),
intent(out) :: dist
5997 real(RP),
parameter :: r180 = 1.0_rp/180.0_rp
5998 real(RP) :: lon1,lon2,lat1,lat2
6001 lon1 = alon * pi * r180
6002 lon2 = blon * pi * r180
6003 lat1 = alat * pi * r180
6004 lat2 = blat * pi * r180
6006 cosd = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(lon2-lon1)
6007 cosd = min( 1._rp,cosd)
6008 cosd = max(-1._rp,cosd)
6010 dist = acos( cosd ) * re
6013 end subroutine com_distll_1
6032 subroutine ensmean_grd(mem, nens, nij, v3d, v2d)
6034 integer,
intent(in) :: mem
6035 integer,
intent(in) :: nens
6036 integer,
intent(in) :: nij
6037 real(RP),
intent(inout) :: v3d(nij,nlev,nens,nv3d)
6038 real(RP),
intent(inout) :: v2d(nij,nens,nv2d)
6040 integer :: i, k, m, n, mmean
6048 v3d(i,k,mmean,n) = v3d(i,k,1,n)
6050 v3d(i,k,mmean,n) = v3d(i,k,mmean,n) + v3d(i,k,m,n)
6052 v3d(i,k,mmean,n) = v3d(i,k,mmean,n) / real(mem, kind=rp)
6058 v2d(i,mmean,n) = v2d(i,1,n)
6060 v2d(i,mmean,n) = v2d(i,mmean,n) + v2d(i,m,n)
6062 v2d(i,mmean,n) = v2d(i,mmean,n) / real(mem, kind=rp)
6067 end subroutine ensmean_grd
6076 subroutine trans_xtoy(elm,ri,rj,rk,lon,lat,v3d,v2d,yobs,qc,stggrd)
6081 integer,
intent(in) :: elm
6082 real(RP),
intent(in) :: ri,rj,rk
6083 real(RP),
intent(in) :: lon,lat
6084 real(RP),
intent(in) :: v3d(:,:,:,:)
6085 real(RP),
intent(in) :: v2d( :,:,:)
6086 real(RP),
intent(out) :: yobs
6087 integer,
intent(out) :: qc
6088 integer,
intent(in),
optional :: stggrd
6089 real(RP) :: u,v,t,q,topo
6091 integer :: stggrd_ = 0
6092 if (
present(stggrd)) stggrd_ = stggrd
6098 case(id_u_obs,id_v_obs)
6099 if (stggrd_ == 1)
then
6100 call itpl_3d(v3d(:,:,:,iv3dd_u),rk,ri-0.5_rp,rj,u)
6101 call itpl_3d(v3d(:,:,:,iv3dd_v),rk,ri,rj-0.5_rp,v)
6103 call itpl_3d(v3d(:,:,:,iv3dd_u),rk,ri,rj,u)
6104 call itpl_3d(v3d(:,:,:,iv3dd_v),rk,ri,rj,v)
6106 if (elm == id_u_obs)
then
6112 call itpl_3d(v3d(:,:,:,iv3dd_t),rk,ri,rj,yobs)
6114 call itpl_3d(v3d(:,:,:,iv3dd_t),rk,ri,rj,yobs)
6115 call itpl_3d(v3d(:,:,:,iv3dd_q),rk,ri,rj,q)
6116 yobs = yobs * (1.0d0 + fvirt * q)
6118 call itpl_3d(v3d(:,:,:,iv3dd_q),rk,ri,rj,yobs)
6120 call itpl_2d(v2d(:,:,iv2dd_t2m),ri,rj,t)
6121 call itpl_2d(v2d(:,:,iv2dd_q2m),ri,rj,q)
6122 call itpl_2d(v2d(:,:,iv2dd_topo),ri,rj,topo)
6123 call itpl_2d(v2d(:,:,iv2dd_ps),ri,rj,yobs)
6124 call prsadj(yobs,rk-topo,t,q)
6125 if( abs(rk-topo) > ps_adjust_thres )
then
6129 call itpl_3d(v3d(:,:,:,iv3dd_rh),rk,ri,rj,yobs)
6135 end subroutine trans_xtoy
6137 subroutine trans_xtoy_radar(elm,radar_lon,radar_lat,radar_z,ri,rj,rk,lon,lat,lev,v3d,v2d,yobs,qc,stggrd)
6143 mapprojection_rotcoef
6146 integer,
intent(in) :: elm
6147 real(RP),
intent(in) :: ri,rj,rk,radar_lon,radar_lat,radar_z
6148 real(RP),
intent(in) :: lon,lat,lev
6149 real(RP),
intent(in) :: v3d(:,:,:,:)
6150 real(RP),
intent(in) :: v2d(:,:,:)
6151 real(RP),
intent(out) :: yobs
6152 integer,
intent(out) :: qc
6153 integer,
intent(in),
optional :: stggrd
6155 integer :: stggrd_ = 0
6157 real(RP) :: qvr,qcr,qrr,qir,qsr,qgr,ur,vr,wr,tr,pr
6158 real(RP) :: dist , dlon , dlat , az , elev , radar_ref,radar_rv
6160 real(RP) :: utmp, vtmp
6161 real(RP) :: rotc(1,1,2)
6162 real(RP) :: lon_tmp(1,1),lat_tmp(1,1)
6164 if (
present(stggrd)) stggrd_ = stggrd
6169 if (stggrd_ == 1)
then
6170 CALL itpl_3d(v3d(:,:,:,iv3dd_u),rk,ri-0.5_rp,rj,ur)
6171 CALL itpl_3d(v3d(:,:,:,iv3dd_v),rk,ri,rj-0.5_rp,vr)
6172 CALL itpl_3d(v3d(:,:,:,iv3dd_w),rk-0.5_rp,ri,rj,wr)
6174 CALL itpl_3d(v3d(:,:,:,iv3dd_u),rk,ri,rj,ur)
6175 CALL itpl_3d(v3d(:,:,:,iv3dd_v),rk,ri,rj,vr)
6176 CALL itpl_3d(v3d(:,:,:,iv3dd_w),rk,ri,rj,wr)
6178 CALL itpl_3d(v3d(:,:,:,iv3dd_t),rk,ri,rj,tr)
6179 CALL itpl_3d(v3d(:,:,:,iv3dd_p),rk,ri,rj,pr)
6180 CALL itpl_3d(v3d(:,:,:,iv3dd_q),rk,ri,rj,qvr)
6181 CALL itpl_3d(v3d(:,:,:,iv3dd_qc),rk,ri,rj,qcr)
6182 CALL itpl_3d(v3d(:,:,:,iv3dd_qr),rk,ri,rj,qrr)
6183 CALL itpl_3d(v3d(:,:,:,iv3dd_qi),rk,ri,rj,qir)
6184 CALL itpl_3d(v3d(:,:,:,iv3dd_qs),rk,ri,rj,qsr)
6185 CALL itpl_3d(v3d(:,:,:,iv3dd_qg),rk,ri,rj,qgr)
6191 dlon = lon-radar_lon
6192 dlat = lat-radar_lat
6193 if ( dlon == 0.0_rp .and. dlat == 0.0_rp )
then
6197 az = r2d*atan2(dlon*cos(radar_lat*d2r),dlat)
6199 if( az < 0 ) az = 360.0_rp + az
6201 call com_distll_1(lon,lat,radar_lon,radar_lat,dist)
6202 elev = r2d*atan2(lev-radar_z,dist)
6204 lon_tmp(1,1) = lon*d2r
6205 lat_tmp(1,1) = lat*d2r
6206 call mapprojection_rotcoef(1, 1, 1, 1, 1, 1, &
6207 lon_tmp(:,:),lat_tmp(:,:),rotc(:,:,1),rotc(:,:,2))
6211 ur = utmp * rotc(1,1,1) - vtmp * rotc(1,1,2)
6212 vr = utmp * rotc(1,1,2) + vtmp * rotc(1,1,1)
6217 CALL calc_ref_vr(qvr,qcr,qrr,qir,qsr,qgr,ur,vr,wr,tr,pr,az,elev,radar_ref,radar_rv)
6220 case(id_radar_ref_obs,id_radar_ref_zero_obs)
6221 if (radar_ref < min_radar_ref)
then
6223 yobs = min_radar_ref_dbz + low_ref_shift
6225 yobs = 10.0_rp * log10(radar_ref)
6227 case(id_radar_vr_obs)
6228 if (radar_ref < min_radar_ref)
then
6237 end subroutine trans_xtoy_radar
6242 subroutine prsadj(p,dz,t,q)
6247 real(RP),
intent(inout) :: p
6248 real(RP),
intent(in) :: dz
6249 real(RP),
intent(in) :: t
6250 real(RP),
intent(in) :: q
6251 real(RP),
parameter :: gamma=5.0d-3
6255 tv = t * (1.0d0 + 0.608d0 * q)
6256 p = p * ((-gamma*dz+tv)/tv)**(gg/(gamma*rd))
6260 end subroutine prsadj
6267 subroutine phys2ijk(p_full,elem,ri,rj,rlev,rk,qc)
6272 real(RP),
intent(in) :: p_full(:,:,:)
6273 integer,
intent(in) :: elem
6274 real(RP),
intent(in) :: ri
6275 real(RP),
intent(in) :: rj
6276 real(RP),
intent(in) :: rlev
6277 real(RP),
intent(out) :: rk
6278 integer,
intent(out) :: qc
6281 real(RP) :: lnps(nlevh,nlonh,nlath)
6282 real(RP) :: plev(nlevh)
6284 integer :: i,j,k, ii, jj, ks
6290 if( ri < 1.0_rp .or. ri > nlonh .or. rj < 1.0_rp .or. rj > nlath )
then
6296 if(elem > 9999)
then
6309 do ii = max(i-1,1), min(i,nlong+xhalo)
6310 do k=1+zhalo,nlev+zhalo
6311 if( p_full(k,ii,jj) >= 0.0_rp )
exit
6317 lnps(:,i-1:i,j-1:j) = log(p_full(:,i-1:i,j-1:j))
6318 call itpl_2d_column(lnps,ri,rj,plev)
6326 if(rk < plev(nlev+zhalo))
then
6327 call itpl_2d(p_full(nlev+zhalo,:,:),ri,rj,ptmp)
6332 if(rk > plev(ks))
then
6333 call itpl_2d(p_full(ks,:,:),ri,rj,ptmp)
6341 do k=ks+1,nlev+zhalo
6342 if(plev(k) < rk)
exit
6345 if (k == nlev+zhalo+1)
then
6347 rk = real(k-2,kind=rp) + ak
6349 ak = (rk - plev(k-1)) / (plev(k) - plev(k-1))
6350 rk = real(k-1,kind=rp) + ak
6355 end subroutine phys2ijk
6362 subroutine phys2ijkz(z_full,ri,rj,rlev,rk,qc)
6367 real(RP),
intent(in) :: z_full(:,:,:)
6368 real(RP),
intent(in) :: ri
6369 real(RP),
intent(in) :: rj
6370 real(RP),
intent(in) :: rlev
6371 real(RP),
intent(out) :: rk
6372 integer,
intent(out) :: qc
6375 real(RP) :: zlev(nlevh)
6377 integer :: i,j,k, ii, jj, ks
6383 if( ri < 1.0_rp .or. ri > nlonh .or. rj < 1.0_rp .or. rj > nlath )
then
6398 do ii = max(i-1,1), min(i,nlong+xhalo)
6399 do k=1+zhalo,nlev+zhalo
6400 if( z_full(k,ii,jj) > -300.0_rp .and. z_full(k,ii,jj) < 10000.0_rp )
exit
6406 call itpl_2d_column(z_full,ri,rj,zlev)
6411 if( rlev > zlev(nlev+zhalo) )
then
6412 call itpl_2d(z_full(nlev+zhalo,:,:),ri,rj,ztmp)
6417 if(rlev < zlev(ks))
then
6418 call itpl_2d(z_full(ks,:,:),ri,rj,ztmp)
6426 do k=ks+1,nlev+zhalo
6427 if(zlev(k) > rlev)
exit
6430 if (k == nlev+zhalo+1)
then
6432 rk = real(k-2,kind=rp) + ak
6434 ak = (rlev - zlev(k-1)) / (zlev(k) - zlev(k-1))
6435 rk = real(k-1,kind=rp) + ak
6439 end subroutine phys2ijkz
6444 subroutine phys2ij(rlon,rlat,rig,rjg)
6451 mapprojection_lonlat2xy
6453 real(RP),
intent(in) :: rlon
6454 real(RP),
intent(in) :: rlat
6455 real(RP),
intent(out) :: rig
6456 real(RP),
intent(out) :: rjg
6462 call mapprojection_lonlat2xy( real(rlon*pi/180.0_rp, kind=rp), &
6463 real(rlat*PI/180.0_RP, kind=rp), rig_rp, rjg_rp )
6464 rig = real((rig_rp - cxg(1)) / dx, kind=rp) + 1.0_rp
6465 rjg = real((rjg_rp - cyg(1)) / dy, kind=rp) + 1.0_rp
6472 end subroutine phys2ij
6477 subroutine itpl_2d(var,ri,rj,var5)
6479 real(RP),
intent(in) :: var(:,:)
6480 real(RP),
intent(in) :: ri
6481 real(RP),
intent(in) :: rj
6482 real(RP),
intent(out) :: var5
6487 ai = ri - real(i-1,kind=rp)
6489 aj = rj - real(j-1,kind=rp)
6492 var5 = var(i ,j-1) * (1-aj) &
6495 var5 = var(i-1,j-1) * (1-ai) * (1-aj) &
6496 & + var(i ,j-1) * ai * (1-aj) &
6497 & + var(i-1,j ) * (1-ai) * aj &
6498 & + var(i ,j ) * ai * aj
6502 end subroutine itpl_2d
6504 subroutine itpl_2d_column(var,ri,rj,var5)
6506 real(RP),
intent(in) :: var(:,:,:)
6507 real(RP),
intent(in) :: ri
6508 real(RP),
intent(in) :: rj
6509 real(RP),
intent(out) :: var5(:)
6515 ai = ri - real(i-1,kind=rp)
6517 aj = rj - real(j-1,kind=rp)
6520 var5(:) = var(:,i ,j-1) * (1-aj) &
6521 & + var(:,i ,j ) * aj
6523 var5(:) = var(:,i-1,j-1) * (1-ai) * (1-aj) &
6524 & + var(:,i ,j-1) * ai * (1-aj) &
6525 & + var(:,i-1,j ) * (1-ai) * aj &
6526 & + var(:,i ,j ) * ai * aj
6530 end subroutine itpl_2d_column
6532 subroutine itpl_3d(var,rk,ri,rj,var5)
6534 real(RP),
intent(in) :: var(:,:,:)
6535 real(RP),
intent(in) :: ri
6536 real(RP),
intent(in) :: rj
6537 real(RP),
intent(in) :: rk
6538 real(RP),
intent(out) :: var5
6539 real(RP) :: ai,aj,ak
6543 ai = ri - real(i-1,kind=rp)
6545 aj = rj - real(j-1,kind=rp)
6547 ak = rk - real(k-1,kind=rp)
6550 var5 = var(k-1,i ,j-1) * (1-aj) * (1-ak) &
6551 & + var(k-1,i ,j ) * aj * (1-ak) &
6552 & + var(k, i ,j-1) * (1-aj) * ak &
6553 & + var(k, i ,j ) * aj * ak
6555 var5 = var(k-1,i-1,j-1) * (1-ai) * (1-aj) * (1-ak) &
6556 & + var(k-1,i ,j-1) * ai * (1-aj) * (1-ak) &
6557 & + var(k-1,i-1,j ) * (1-ai) * aj * (1-ak) &
6558 & + var(k-1,i ,j ) * ai * aj * (1-ak) &
6559 & + var(k, i-1,j-1) * (1-ai) * (1-aj) * ak &
6560 & + var(k, i ,j-1) * ai * (1-aj) * ak &
6561 & + var(k, i-1,j ) * (1-ai) * aj * ak &
6562 & + var(k, i ,j ) * ai * aj * ak
6566 end subroutine itpl_3d
6568 subroutine merge_sort_parallel(n, array)
6570 integer(8),
intent(in) :: n
6571 integer(8),
intent(inout) :: array(n)
6576 omp_nested = omp_get_nested()
6577 maxnest = floor(log(dble(omp_get_max_threads())) / log(2.0d0))
6578 call omp_set_nested(.true.)
6579 call merge_sort_2threads(n, array, 0, maxnest)
6580 call omp_set_nested(omp_nested)
6583 end subroutine merge_sort_parallel
6585 recursive subroutine merge_sort_2threads(n, array, nest, maxnest)
6589 integer(8),
intent(in) :: n
6590 integer,
intent(in) :: nest, maxnest
6591 integer(8),
intent(inout) :: array(n)
6592 integer(8) :: asize(2)
6593 integer(8),
allocatable :: tmpary(:)
6594 integer,
parameter:: nmin = 4
6597 if(mod(n,2_8) .ne. 0) asize(1) = asize(1) + 1
6599 if(nest < maxnest)
then
6601 tmpary(1:asize(1)) = array(1:asize(1))
6602 if(asize(1) > nmin)
then
6603 call merge_sort_2threads(asize(1), tmpary(1:asize(1)), nest + 1, maxnest)
6607 tmpary((asize(1) + 1):n) = array((asize(1) + 1):n)
6608 if(asize(2) > nmin)
then
6609 call merge_sort_2threads(asize(2), tmpary((asize(1) + 1):n), nest + 1, maxnest)
6613 call merge_2threads(asize(1), tmpary(1:asize(1)), asize(2), tmpary((asize(1) + 1):n), n, array, nest, maxnest)
6619 call merge(asize(1), tmpary(1:asize(1)), asize(2), tmpary((asize(1) + 1):n), n, array)
6621 end subroutine merge_sort_2threads
6623 recursive subroutine merge_2threads(n1, ary1, n2, ary2, n3, ary3, nest, maxnest)
6624 integer(8),
intent(in) :: n1, ary1(n1), n2, ary2(n2), n3
6625 integer,
intent(in) :: nest, maxnest
6626 integer(8),
intent(inout) :: ary3(n3)
6628 integer,
parameter :: threshold = 10
6630 if(nest >= maxnest .or. n1 < threshold .or. n2 < threshold)
then
6631 call merge(n1, ary1, n2, ary2, n3, ary3)
6636 m = binary_search_i8(n2, ary2, ary1(k))
6637 call merge_2threads(k, ary1(1:k), m, ary2(1:m), k + m, ary3(1:(k + m)), nest + 1, maxnest)
6638 call merge_2threads(n1 - k, ary1((k + 1):n1), n2 - m, ary2((m + 1):n2), n3 - (k + m), ary3((k + m + 1):n3), nest + 1, maxnest)
6640 end subroutine merge_2threads
6642 subroutine merge(n1, ary1, n2, ary2, n3, ary3)
6643 integer(8),
intent(in) :: n1, ary1(n1), n2, ary2(n2), n3
6644 integer(8),
intent(inout) :: ary3(n3)
6649 if(ary1(i) < ary2(j))
then
6652 ary3((k + 1):n3) = ary2(j:n2)
6659 ary3((k + 1):n3) = ary1(i:n1)
6665 end subroutine merge
6667 recursive subroutine merge_sort_mpi(n, array, comm)
6672 integer(8),
intent(in) :: n
6673 integer,
intent(in) :: comm
6674 integer(8),
intent(inout) :: array(n)
6675 integer(8) :: asize(2)
6676 integer(8),
allocatable :: tmpary(:)
6677 integer :: parent, child, procs(2), newcomm, tag, nprocs, rank
6680 call mpi_comm_size(comm, nprocs, ierr)
6681 call mpi_comm_rank(comm, rank, ierr)
6683 procs(1) = nprocs / 2
6684 procs(2) = nprocs - procs(1)
6686 child = parent + procs(1)
6689 if(mod(n,2_8) .ne. 0) asize(1) = asize(1) + 1
6692 if(rank < child)
then
6693 call mpi_comm_split(comm, 0, rank, newcomm, ierr)
6694 tmpary(1:asize(1)) = array(1:asize(1))
6695 if(procs(1) > 1)
then
6696 call merge_sort_mpi(asize(1), tmpary(1:asize(1)), newcomm)
6701 call mpi_comm_split(comm, 1, rank, newcomm, ierr)
6702 tmpary((asize(1) + 1):n) = array((asize(1) + 1):n)
6703 if(procs(2) > 1)
then
6704 call merge_sort_mpi(asize(2), tmpary((asize(1) + 1):n), newcomm)
6710 call merge_mpi_no_nest(asize(1), tmpary(1:asize(1)), asize(2), tmpary((asize(1) + 1):n), n, array, parent, child, comm)
6712 end subroutine merge_sort_mpi
6714 recursive subroutine merge_mpi(n1, ary1, n2, ary2, n3, ary3, parent, child, nprocs, comm)
6717 integer(8),
intent(in) :: n1, ary1(n1), n2, ary2(n2), n3
6718 integer,
intent(in) :: parent, child, nprocs, comm
6719 integer(8),
intent(inout) :: ary3(n3)
6720 integer(8) k, m, pivot
6721 integer,
parameter :: threshold = 10
6722 integer procs(2), grandchild1, grandchild2, rank, newcomm
6723 integer,
parameter :: tag_p2c = 12345
6724 integer,
parameter :: tag_p2g = 23451
6725 integer,
parameter :: tag_c2g = 34512
6726 integer,
parameter :: tag_p2h = 45123
6727 integer,
parameter :: tag_c2h = 51234
6731 call mpi_comm_rank(comm, rank, ierr)
6732 procs(1) = child - parent
6733 procs(2) = nprocs - procs(1)
6736 if(rank == parent) pivot = ary1(k)
6737 call mpi_bcast(pivot, 1, mpi_integer8, parent, comm, ierr)
6738 if(rank == child) m = binary_search_i8(n2, ary2, pivot)
6739 call mpi_bcast(m, 1, mpi_integer8, child, comm, ierr)
6741 if(procs(1) > 1 .and. n1 >= threshold)
then
6742 grandchild1 = parent + (procs(1) / 2)
6744 grandchild1 = parent
6746 if(procs(2) > 1 .and. n2 >= threshold)
then
6747 grandchild2 = child + (procs(2) / 2)
6752 if(rank >= parent .and. rank < child)
then
6753 call mpi_comm_split(comm, 0, rank, newcomm, ierr)
6755 if(rank == parent)
then
6756 if(rank == grandchild1)
then
6757 call mpi_sendrecv(ary1((k + 1):n1), int(n1 - k), mpi_integer8, grandchild2, tag_p2h, &
6758 & ary2(1:m), int(m), mpi_integer8, child, tag_c2g, &
6759 & comm, mpi_status_ignore, ierr)
6761 call mpi_send(ary1((k + 1):n1), int(n1 - k,4), mpi_integer8, grandchild2, &
6762 & tag_p2h, comm, ierr)
6764 else if(rank == grandchild1)
then
6765 call mpi_recv(ary2(1:m), int(m), mpi_integer8, child, &
6766 & tag_c2g, comm, mpi_status_ignore, ierr)
6769 if(procs(1) > 1 .and. n1 >= threshold)
then
6770 call merge_mpi(k, ary1(1:k), m, ary2(1:m), k + m, ary3(1:(k + m)), 0, grandchild1 - parent, procs(1), newcomm)
6772 call merge(k, ary1(1:k), m, ary2(1:m), k + m, ary3(1:(k + m)))
6774 if(rank == parent)
then
6775 call mpi_recv(ary3((k + m + 1):n3), int(n3 - (k + m)), mpi_integer8, child, &
6776 & tag_p2c, comm, mpi_status_ignore, ierr)
6778 else if(rank >= child .and. rank < parent + nprocs)
then
6779 call mpi_comm_split(comm, 1, rank, newcomm, ierr)
6781 if(rank == child)
then
6782 if(rank == grandchild2)
then
6783 call mpi_sendrecv(ary2(1:m), int(m), mpi_integer8, grandchild1, tag_c2g, &
6784 & ary1((k + 1):n1), int(n1 - k), mpi_integer8, parent, tag_p2h, &
6785 & comm, mpi_status_ignore, ierr)
6787 call mpi_send(ary2(1:m), int(m,4), mpi_integer8, grandchild1, &
6788 & tag_c2g, comm, ierr)
6790 else if(rank == grandchild2)
then
6791 call mpi_recv(ary1((k + 1):n1), int(n1 - k), mpi_integer8, parent, &
6792 & tag_p2h, comm, mpi_status_ignore, ierr)
6795 if(procs(2) > 1 .and. n2 >= threshold)
then
6796 call merge_mpi(n2 - m, ary2((m + 1):n2), n1 - k, ary1((k + 1):n1), &
6797 & n3 - (k + m), ary3((k + m + 1):n3), 0, grandchild2 - child, procs(2), newcomm)
6799 call merge(n1 - k, ary1((k + 1):n1), n2 - m, ary2((m + 1):n2), &
6800 & n3 - (k + m), ary3((k + m + 1):n3))
6802 if(rank == child)
then
6803 call mpi_send(ary3((k + m + 1):n3), int(n3 - (k + m)), mpi_integer8, parent, &
6804 & tag_p2c, comm, ierr)
6807 call mpi_comm_split(comm, 3, rank, newcomm, ierr)
6808 write(*, *)
"something wrong in merge_mpi"
6810 end subroutine merge_mpi
6812 recursive subroutine merge_mpi_no_nest(n1, ary1, n2, ary2, n3, ary3, parent, child, comm)
6815 integer(8),
intent(in) :: n1, ary1(n1), n2, ary2(n2), n3
6816 integer,
intent(in) :: parent, child, comm
6817 integer(8),
intent(inout) :: ary3(n3)
6818 integer(8) k, m, pivot
6819 integer,
parameter :: threshold = 10
6821 integer,
parameter :: tag_p2c = 12345
6823 integer(8) pivot_dummy
6829 call mpi_comm_rank(comm, rank, ierr)
6832 if(rank == parent)
then
6835 call mpi_send(pivot, 1, mpi_integer8, child, tag_p2c, comm, ierr)
6837 call mpi_recv(m, 1, mpi_integer8, child, tag_p2c, comm, mpi_status_ignore, ierr)
6838 call mpi_sendrecv(ary1((k + 1):n1), int(n1 - k), mpi_integer8, child, tag_p2c, &
6839 & ary2(1:m), int(m), mpi_integer8, child, tag_p2c, &
6840 & comm, mpi_status_ignore, ierr)
6841 call merge(k, ary1(1:k), m, ary2(1:m), k + m, ary3(1:(k + m)))
6842 call mpi_recv(ary3((k + m + 1):n3), int(n3 - (k + m)), mpi_integer8, child, &
6843 & tag_p2c, comm, mpi_status_ignore, ierr)
6844 else if(rank == child)
then
6845 call mpi_recv(pivot, 1, mpi_integer8, parent, tag_p2c, comm, mpi_status_ignore, ierr)
6846 m = binary_search_i8(n2, ary2, pivot)
6847 call mpi_send(m, 1, mpi_integer8, parent, tag_p2c, comm, ierr)
6848 call mpi_sendrecv(ary2(1:m), int(m), mpi_integer8, parent, tag_p2c, &
6849 & ary1((k + 1):n1), int(n1 - k), mpi_integer8, parent, tag_p2c, &
6850 & comm, mpi_status_ignore, ierr)
6851 call merge(n1 - k, ary1((k + 1):n1), n2 - m, ary2((m + 1):n2), &
6852 & n3 - (k + m), ary3((k + m + 1):n3))
6853 call mpi_send(ary3((k + m + 1):n3), int(n3 - (k + m)), mpi_integer8, parent, &
6854 & tag_p2c, comm, ierr)
6856 end subroutine merge_mpi_no_nest
6861 subroutine scatter_grd_mpi(nrank,v3dg,v2dg,v3d,v2d)
6862 integer,
intent(in) :: nrank
6863 real(RP),
intent(in) :: v3dg(nlev,nlon,nlat,nv3d)
6864 real(RP),
intent(in) :: v2dg(nlon,nlat,nv2d)
6865 real(RP),
intent(out) :: v3d(nij1,nlev,nv3d)
6866 real(RP),
intent(out) :: v2d(nij1,nv2d)
6867 real(RP) :: bufs(nij1max,nlevall,NPRC_ENS)
6868 real(RP) :: bufr(nij1max,nlevall)
6869 integer :: j,k,n,ierr,ns,nr
6871 ns = nij1max * nlevall
6873 if( rank_ens == nrank )
then
6878 call grd_to_buf( nprc_ens, v3dg(k,:,:,n), bufs(:,j,:) )
6884 call grd_to_buf( nprc_ens, v2dg(:,:,n),bufs(:,j,:) )
6888 CALL mpi_scatter( bufs, ns, datatype, bufr, nr, datatype, nrank, comm_ens, ierr )
6894 v3d(:,k,n) = real( bufr(1:nij1,j), kind=rp )
6900 v2d(:,n) = real( bufr(1:nij1,j), kind=rp )
6904 end subroutine scatter_grd_mpi
6909 subroutine scatter_grd_mpi_all2all(mstart,mend,v3dg,v2dg,v3d,v2d)
6910 integer,
intent(in) :: mstart,mend
6911 real(RP),
intent(in) :: v3dg(nlev,nlon,nlat,nv3d)
6912 real(RP),
intent(in) :: v2dg(nlon,nlat,nv2d)
6913 real(RP),
intent(inout) :: v3d(nij1,nlev,nens,nv3d)
6914 real(RP),
intent(inout) :: v2d(nij1,nens,nv2d)
6915 real(RP) :: bufs(nij1max,nlevall,NPRC_ENS)
6916 real(RP) :: bufr(nij1max,nlevall,NPRC_ENS)
6917 integer :: k,n,j,m,mcount,ierr
6918 integer :: ns(NPRC_ENS),nst(NPRC_ENS),nr(NPRC_ENS),nrt(NPRC_ENS)
6920 mcount = mend - mstart + 1
6922 if(rank_ens < mcount)
then
6927 call grd_to_buf(nprc_ens,v3dg(k,:,:,n),bufs(:,j,:))
6932 call grd_to_buf(nprc_ens,v2dg(:,:,n),bufs(:,j,:))
6936 if(mcount == nprc_ens)
then
6937 call mpi_alltoall(bufs, nij1max*nlevall, datatype, &
6938 bufr, nij1max*nlevall, datatype, comm_ens, ierr)
6940 call set_all2allv_counts(mcount,nij1max*nlevall,nprc_ens,nr,nrt,ns,nst)
6941 call mpi_alltoallv(bufs, ns, nst, datatype, &
6942 bufr, nr, nrt, datatype, comm_ens, ierr)
6950 v3d(:,k,m,n) = real(bufr(1:nij1,j,m-mstart+1),kind=rp)
6955 v2d(:,m,n) = real(bufr(1:nij1,j,m-mstart+1),kind=rp)
6960 end subroutine scatter_grd_mpi_all2all
6965 subroutine gather_grd_mpi_all2all(mstart,mend,v3d,v2d,v3dg,v2dg)
6966 integer,
intent(in) :: mstart,mend
6967 real(RP),
intent(in) :: v3d(nij1,nlev,nens,nv3d)
6968 real(RP),
intent(in) :: v2d(nij1,nens,nv2d)
6969 real(RP),
intent(out) :: v3dg(nlev,nlon,nlat,nv3d)
6970 real(RP),
intent(out) :: v2dg(nlon,nlat,nv2d)
6971 real(RP) :: bufs(nij1max,nlevall,NPRC_ENS)
6972 real(RP) :: bufr(nij1max,nlevall,NPRC_ENS)
6973 integer :: k,n,j,m,mcount,ierr
6974 integer :: ns(NPRC_ENS),nst(NPRC_ENS),nr(NPRC_ENS),nrt(NPRC_ENS)
6976 mcount = mend - mstart + 1
6983 bufs(1:nij1,j,m-mstart+1) = real(v3d(:,k,m,n),kind=rp)
6988 bufs(1:nij1,j,m-mstart+1) = real(v2d(:,m,n),kind=rp)
6992 if(mcount == nprc_ens)
then
6993 call mpi_alltoall(bufs, nij1max*nlevall, datatype, &
6994 bufr, nij1max*nlevall, datatype, comm_ens, ierr)
6996 call set_all2allv_counts(mcount,nij1max*nlevall,nprc_ens,ns,nst,nr,nrt)
6997 call mpi_alltoallv(bufs, ns, nst, datatype, &
6998 bufr, nr, nrt, datatype, comm_ens, ierr)
7001 if(rank_ens < mcount)
then
7006 call buf_to_grd(nprc_ens,bufr(:,j,:),v3dg(k,:,:,n))
7011 call buf_to_grd(nprc_ens,bufr(:,j,:),v2dg(:,:,n))
7016 end subroutine gather_grd_mpi_all2all
7021 subroutine set_all2allv_counts(mcount,ngpblock,np,n_ens,nt_ens,n_mem,nt_mem)
7022 integer,
intent(in) :: mcount,ngpblock
7023 integer,
intent(in) :: np
7024 integer,
intent(out) :: n_ens(np),nt_ens(np),n_mem(np),nt_mem(np)
7033 if(rank_ens+1 == p)
then
7038 nt_ens(p) = nt_ens(p-1) + n_ens(p-1)
7039 nt_mem(p) = nt_mem(p-1) + n_mem(p-1)
7043 end subroutine set_all2allv_counts
7048 subroutine grd_to_buf(np,grd,buf)
7051 integer,
intent(in) :: np
7052 real(RP),
intent(in) :: grd(nlon,nlat)
7053 real(RP),
intent(out) :: buf(nij1max,np)
7054 integer :: i,j,m,ilon,ilat
7057 do i = 1, nij1node(m)
7058 j = m-1 + np * (i-1)
7059 ilon = mod(j,nlon) + 1
7060 ilat = (j-ilon+1) / nlon + 1
7062 buf(i,m) = grd(ilon,ilat)
7067 if( nij1node(m) < nij1max ) buf(nij1max,m) = undef
7071 end subroutine grd_to_buf
7076 subroutine buf_to_grd(np,buf,grd)
7077 integer,
intent(in) :: np
7078 real(RP),
intent(in) :: buf(nij1max,np)
7079 real(RP),
intent(out) :: grd(nlon,nlat)
7080 integer :: i,j,m,ilon,ilat
7083 do i = 1,nij1node(m)
7084 j = m-1 + np * (i-1)
7085 ilon = mod(j,nlon) + 1
7086 ilat = (j-ilon+1) / nlon + 1
7087 grd(ilon,ilat) = buf(i,m)
7092 end subroutine buf_to_grd
7103 subroutine calc_z_grd(nij, topo, z)
7111 integer,
intent(in) :: nij
7112 real(RP),
intent(in) :: topo(nij)
7113 real(RP),
intent(out) :: z(nij,nlev)
7125 end subroutine calc_z_grd
7127 subroutine qc_indexing_and_packing(&
7128 & na, nr, ne, ze, vr, radlon, radlat, radz, & ! input spherical
7129 & qcflag, input_is_dbz, attenuation, & ! input spherical
7130 & nlon, nlat, nlev, lon, lat, z, dlon, dlat, dz, & ! input cartesian
7131 & missing, & ! input param
7132 & lon0, lat0, & ! input param
7133 & nobs_each_elev, packed_grid, packed_data, packed_attn)
7139 integer,
intent(in) :: na, nr, ne
7140 real(RP),
intent(in) :: ze(na * nr, ne), vr(na * nr, ne)
7141 real(RP),
intent(in),
dimension(na * nr, ne) :: radlon, radlat, radz
7142 real(RP),
intent(in) :: qcflag(na * nr, ne)
7143 logical,
intent(in) :: input_is_dbz
7144 real(RP),
intent(in) :: lon0, lat0
7145 real(RP),
intent(in),
dimension(na * nr, ne) :: attenuation
7146 integer,
intent(in) :: nlon, nlat, nlev
7147 real(RP),
intent(in) :: lon(nlon), lat(nlat), z(nlev)
7148 real(RP),
intent(in) :: dlon, dlat, dz, missing
7149 integer(8),
intent(out) :: nobs_each_elev(ne)
7150 integer(8),
intent(out) :: packed_grid(na * nr * ne)
7151 real(RP),
intent(out) :: packed_data(5, na * nr * ne)
7152 logical,
intent(out) :: packed_attn(na * nr * ne)
7153 REAL(RP) :: ri, rj, rk, dlon_inv, dlat_inv, dz_inv, qced_ze, qced_vr
7157 real(RP) :: lon_coef, lat_coef, dist_i, dist_j
7158 real(RP) :: vr_min_dist_square
7160 lon_coef = d2r * (cos(lat0 * d2r) * radius)
7161 lat_coef = d2r * radius
7162 vr_min_dist_square = vr_min_dist * vr_min_dist
7164 dlon_inv = 1.0d0 / dlon
7165 dlat_inv = 1.0d0 / dlat
7178 IF(input_is_dbz .and. (qced_ze .ne. missing)) qced_ze = 10.0d0 ** (qced_ze * 0.1d0)
7182 if(qced_ze .LE. missing)
then
7186 if(qced_ze .LT. minz)
then
7189 if(qced_ze .LT. min_radar_ref_vr)
then
7199 if(use_qcflag .and. (qcflag(i, ie) .GE. 900.0d0))
then
7204 if(qced_ze == missing) cycle
7208 ri = (radlon(i, ie) - lon(1)) * dlon_inv
7209 rj = (radlat(i, ie) - lat(1)) * dlat_inv
7210 rk = (radz(i, ie) - z(1) ) * dz_inv
7212 if(ri < 0 .or. ri >= nlon .or. rj < 0 .or. rj >= nlat .or. rk < 0 .or. rk >= nlev) cycle
7214 if (qced_vr > missing)
then
7215 if (radz(i, ie) < radar_zmin)
then
7218 dist_i = (radlon(i, ie) - lon0) * lon_coef
7219 dist_j = (radlat(i, ie) - lat0) * lat_coef
7220 if (dist_i * dist_i + dist_j * dist_j < vr_min_dist_square)
then
7226 packed_grid(idx) = aint(ri) + (aint(rj) + aint(rk) * nlat) * nlon + 1
7227 packed_data(1, idx) = qced_ze
7228 packed_data(2, idx) = qced_vr
7229 packed_data(3, idx) = radlon(i, ie)
7230 packed_data(4, idx) = radlat(i, ie)
7231 packed_data(5, idx) = radz(i, ie)
7232 packed_attn(idx) = ((.not. use_attenuation) .or. &
7233 & (attenuation(i, ie) > attenuation_threshold))
7235 nobs_each_elev(ie) = nobs_each_elev(ie) + 1
7239 end subroutine qc_indexing_and_packing
7244 integer function uid_obs(id_obs)
7246 integer,
intent(in) :: id_obs
7266 case(id_radar_ref_obs)
7268 case(id_radar_ref_zero_obs)
7270 case(id_radar_vr_obs)
7272 case(id_radar_prh_obs)
7287 end function uid_obs
7292 integer function uid_obs_varlocal(id_obs)
7294 integer,
intent(in) :: id_obs
7298 case(id_u_obs, id_v_obs)
7299 uid_obs_varlocal = 1
7300 case(id_t_obs, id_tv_obs)
7301 uid_obs_varlocal = 2
7302 case(id_q_obs, id_rh_obs)
7303 uid_obs_varlocal = 3
7305 uid_obs_varlocal = 4
7307 uid_obs_varlocal = 5
7308 case(id_tclon_obs, id_tclat_obs, id_tcmip_obs)
7309 uid_obs_varlocal = 6
7310 case(id_radar_ref_obs, id_radar_ref_zero_obs, id_radar_prh_obs)
7311 uid_obs_varlocal = 7
7312 case(id_radar_vr_obs)
7313 uid_obs_varlocal = 8
7315 uid_obs_varlocal = 9
7317 uid_obs_varlocal = -1
7319 end function uid_obs_varlocal
7321 function binary_search_i8(n, ary, val)
7322 integer(8) binary_search_i8
7323 integer(8),
intent(in) :: n
7324 integer(8),
intent(in) :: ary(n), val
7325 integer(8) pivot, nmax, nmin
7329 if(ary(1) < ary(n))
then
7330 do while(nmax > nmin + 1)
7331 pivot = nmin + (nmax - nmin) / 2
7332 if(val == ary(pivot))
then
7335 else if(val < ary(pivot))
then
7342 do while(nmax > nmin + 1)
7343 pivot = nmin + (nmax - nmin) / 2
7344 if(val == ary(pivot))
then
7347 else if(val > ary(pivot))
then
7354 binary_search_i8 = nmin
7355 end function binary_search_i8
7360 subroutine rank_1d_2d(rank, rank_i, rank_j)
7363 integer,
intent(in) :: rank
7364 integer,
intent(out) :: rank_i, rank_j
7371 end subroutine rank_1d_2d
7376 subroutine rank_2d_1d(rank_i, rank_j, rank)
7378 integer,
intent(in) :: rank_i, rank_j
7379 integer,
intent(out) :: rank
7381 rank = rank_j * nproc_x + rank_i
7384 end subroutine rank_2d_1d
7396 subroutine rij_rank(ig, jg, rank)
7398 real(RP),
intent(in) :: ig
7399 real(RP),
intent(in) :: jg
7400 integer,
intent(out) :: rank
7401 integer :: rank_i, rank_j
7403 if (ig < real(start_x,kind=rp) .or. ig > real(nlon*nproc_x+xhalo,kind=rp) .or. &
7404 jg < real(start_y,kind=rp) .or. jg > real(nlat*nproc_y+yhalo,kind=rp))
then
7406 if (.not.(start_x == end_x .and. (ig < real(start_x,kind=rp) .or. ig > real(nlong+xhalo,kind=rp) )))
then
7414 rank_i = ceiling((ig-real(xhalo,kind=rp)-0.5_rp) / real(nlon,kind=rp)) - 1
7415 rank_j = ceiling((jg-real(yhalo,kind=rp)-0.5_rp) / real(nlat,kind=rp)) - 1
7416 call rank_2d_1d(rank_i, rank_j, rank)
7419 end subroutine rij_rank
7424 subroutine rij_g2l(rank, ig, jg, il, jl)
7426 integer,
intent(in) :: rank
7427 real(RP),
intent(in) :: ig
7428 real(RP),
intent(in) :: jg
7429 real(RP),
intent(out) :: il
7430 real(RP),
intent(out) :: jl
7432 integer :: rank_i, rank_j
7435 call rank_1d_2d( rank, rank_i, rank_j )
7436 il = ig - real(rank_i * nlon, kind=rp)
7437 jl = jg - real(rank_j * nlat, kind=rp)
7440 end subroutine rij_g2l
7445 subroutine ij_obsgrd( ctype, ri, rj, ogi, ogj )
7447 integer,
intent(in) :: ctype
7448 real(RP),
intent(in) :: ri, rj
7449 integer,
intent(out) :: ogi, ogj
7450 real(RP) :: ril, rjl
7452 call rij_g2l( rank_lcl, ri, rj, ril, rjl )
7453 ogi = ceiling( ( ril - real( xhalo, kind=rp ) - 0.5 ) * real( obsgrd(ctype)%ngrd_i, kind=rp ) / real( nlon, kind=rp ) )
7454 ogj = ceiling( ( rjl - real( yhalo, kind=rp ) - 0.5 ) * real( obsgrd(ctype)%ngrd_j, kind=rp ) / real( nlat, kind=rp ) )
7457 end subroutine ij_obsgrd
7462 subroutine ij_obsgrd_ext( ctype, ri, rj, ogi, ogj )
7464 integer,
intent(in) :: ctype
7465 real(RP),
intent(in) :: ri, rj
7466 integer,
intent(out) :: ogi, ogj
7468 real(RP) :: ril, rjl
7471 call rij_g2l( rank_lcl, ri, rj, ril, rjl )
7472 ogi = ceiling( ( ril - real( xhalo, kind=rp ) - 0.5 ) * real( obsgrd(ctype)%ngrd_i, kind=rp ) / real( nlon, kind=rp ) ) &
7473 + obsgrd(ctype)%ngrdsch_i
7474 ogj = ceiling( ( rjl - real( yhalo, kind=rp ) - 0.5 ) * real( obsgrd(ctype)%ngrd_j, kind=rp ) / real( nlat, kind=rp ) ) &
7475 + obsgrd(ctype)%ngrdsch_j
7478 end subroutine ij_obsgrd_ext
7483 subroutine obs_choose(ctype, proc, imin, imax, jmin, jmax, nn, nobs_use)
7485 integer,
intent(in) :: ctype
7486 integer,
intent(in) :: proc
7487 integer,
intent(in) :: imin, imax, jmin, jmax
7488 integer,
intent(inout) :: nn
7489 integer,
intent(inout),
optional :: nobs_use(:)
7492 if (imin > imax .or. jmin > jmax)
return
7493 if (obsgrd(ctype)%tot(proc) == 0)
return
7496 if (
present(nobs_use))
then
7497 do n = obsgrd(ctype)%ac(imin-1,j,proc)+1, obsgrd(ctype)%ac(imax,j,proc)
7502 nn = nn + obsgrd(ctype)%ac(imax,j,proc) - obsgrd(ctype)%ac(imin-1,j,proc)
7507 end subroutine obs_choose
7512 subroutine obs_choose_ext(ctype, imin, imax, jmin, jmax, nn, nobs_use)
7514 integer,
intent(in) :: ctype
7515 integer,
intent(in) :: imin, imax, jmin, jmax
7516 integer,
intent(inout) :: nn
7517 integer,
intent(out),
optional :: nobs_use(:)
7522 if (imin > imax .or. jmin > jmax)
return
7523 if (obsgrd(ctype)%tot_ext == 0)
return
7526 if (
present(nobs_use))
then
7527 do n = obsgrd(ctype)%ac_ext(imin-1,j)+1, obsgrd(ctype)%ac_ext(imax,j)
7532 nn = nn + obsgrd(ctype)%ac_ext(imax,j) - obsgrd(ctype)%ac_ext(imin-1,j)
7537 end subroutine obs_choose_ext
7539 subroutine obs_info_allocate(obs, extended)
7541 type(
obs_info),
intent(inout) :: obs
7542 logical,
optional,
intent(in) :: extended
7544 call obs_info_deallocate(obs)
7546 allocate( obs%elm (obs%nobs) )
7547 allocate( obs%lon (obs%nobs) )
7548 allocate( obs%lat (obs%nobs) )
7549 allocate( obs%lev (obs%nobs) )
7550 allocate( obs%dat (obs%nobs) )
7551 allocate( obs%err (obs%nobs) )
7552 allocate( obs%typ (obs%nobs) )
7553 allocate( obs%dif (obs%nobs) )
7564 if (
present(extended))
then
7566 allocate( obs%ri (obs%nobs) )
7567 allocate( obs%rj (obs%nobs) )
7568 allocate( obs%rank (obs%nobs) )
7577 end subroutine obs_info_allocate
7579 subroutine obs_info_deallocate(obs)
7581 type(
obs_info),
intent(inout) :: obs
7583 if(
allocated(obs%elm) )
deallocate(obs%elm)
7584 if(
allocated(obs%lon) )
deallocate(obs%lon)
7585 if(
allocated(obs%lat) )
deallocate(obs%lat)
7586 if(
allocated(obs%lev) )
deallocate(obs%lev)
7587 if(
allocated(obs%dat) )
deallocate(obs%dat)
7588 if(
allocated(obs%err) )
deallocate(obs%err)
7589 if(
allocated(obs%typ) )
deallocate(obs%typ)
7590 if(
allocated(obs%dif) )
deallocate(obs%dif)
7591 if(
allocated(obs%ri ) )
deallocate(obs%ri)
7592 if(
allocated(obs%rj ) )
deallocate(obs%rj)
7593 if(
allocated(obs%rank))
deallocate(obs%rank)
7596 end subroutine obs_info_deallocate
7598 subroutine obs_da_value_allocate( obsda, member )
7601 integer,
intent(in) :: member
7603 call obs_da_value_deallocate( obsda )
7605 allocate( obsda%set (obsda%nobs) )
7606 allocate( obsda%idx (obsda%nobs) )
7607 allocate( obsda%key (obsda%nobs) )
7608 allocate( obsda%val (obsda%nobs) )
7609 allocate( obsda%qc (obsda%nobs) )
7611 allocate( obsda%tm (obsda%nobs) )
7612 allocate( obsda%pm (obsda%nobs) )
7613 allocate( obsda%qv (obsda%nobs) )
7615 obsda%nobs_in_key = 0
7625 if (member > 0)
then
7626 allocate( obsda%ensval(member,obsda%nobs) )
7627 obsda%ensval = 0.0d0
7629 allocate( obsda%eqv(member,obsda%nobs) )
7634 end subroutine obs_da_value_allocate
7636 subroutine obs_da_value_deallocate( obsda )
7640 obsda%nobs_in_key = 0
7642 if(
allocated(obsda%set ) )
deallocate(obsda%set )
7643 if(
allocated(obsda%idx ) )
deallocate(obsda%idx )
7644 if(
allocated(obsda%key ) )
deallocate(obsda%key )
7645 if(
allocated(obsda%val ) )
deallocate(obsda%val )
7646 if(
allocated(obsda%ensval) )
deallocate(obsda%ensval)
7647 if(
allocated(obsda%qc ) )
deallocate(obsda%qc )
7649 if(
allocated(obsda%tm ) )
deallocate(obsda%tm)
7650 if(
allocated(obsda%pm ) )
deallocate(obsda%pm)
7651 if(
allocated(obsda%eqv ) )
deallocate(obsda%eqv)
7652 if(
allocated(obsda%qv ) )
deallocate(obsda%qv)
7655 end subroutine obs_da_value_deallocate
7657 subroutine obs_da_value_allreduce( obsda )
7661 real(RP),
allocatable :: ensval_bufs(:,:)
7662 real(RP),
allocatable :: ensval_bufr(:,:)
7663 real(RP),
allocatable :: ensval_bufs2(:,:)
7664 real(RP),
allocatable :: ensval_bufr2(:,:)
7666 integer :: cntr(NPRC_ENS)
7667 integer :: dspr(NPRC_ENS)
7668 integer :: current_shape(2)
7669 integer :: ie, it, im, imb, ierr
7671 if( obsda%nobs <= 0 )
then
7678 cntr(ie) = cntr(ie) + 1
7680 allocate( ensval_bufs(obsda%nobs, cntr(rank_ens+1)) )
7681 allocate( ensval_bufr(obsda%nobs, nprc_ens) )
7682 allocate( ensval_bufs2(obsda%nobs, cntr(rank_ens+1)) )
7683 allocate( ensval_bufr2(obsda%nobs, nprc_ens) )
7685 do im = 1, cntr(rank_ens+1)
7686 ensval_bufs(:,im) = obsda%ensval(im,:)
7687 ensval_bufs2(:,im) = obsda%eqv(im,:)
7690 cntr(:) = cntr(:) * obsda%nobs
7691 cnts = cntr(rank_ens+1)
7694 dspr(ie) = dspr(ie-1) + cntr(ie-1)
7697 call mpi_allgatherv( ensval_bufs, cnts, datatype, ensval_bufr, cntr, dspr, datatype, comm_ens, ierr )
7698 call mpi_allgatherv( ensval_bufs2, cnts, datatype, ensval_bufr2, cntr, dspr, datatype, comm_ens, ierr )
7700 current_shape = shape(obsda%ensval)
7701 if (current_shape(1) < nprc_ens)
then
7702 deallocate (obsda%ensval)
7703 allocate (obsda%ensval(nprc_ens, obsda%nobs))
7704 deallocate (obsda%eqv)
7705 allocate (obsda%eqv(nprc_ens, obsda%nobs))
7709 obsda%ensval(ie,:) = ensval_bufr(:,ie)
7710 obsda%eqv(ie,:) = ensval_bufr2(:,ie)
7712 deallocate(ensval_bufs, ensval_bufr)
7713 deallocate(ensval_bufs2, ensval_bufr2)
7716 if( nprc_ens > 1 )
then
7717 call mpi_allreduce( mpi_in_place, obsda%qc(:), obsda%nobs, mpi_integer, mpi_max, comm_ens, ierr )
7718 call mpi_allreduce( mpi_in_place, obsda%tm(:), obsda%nobs, datatype, mpi_sum, comm_ens, ierr )
7719 call mpi_allreduce( mpi_in_place, obsda%pm(:), obsda%nobs, datatype, mpi_sum, comm_ens, ierr )
7721 obsda%tm = obsda%tm / real( nprc_ens, kind=rp )
7722 obsda%pm = obsda%pm / real( nprc_ens, kind=rp )
7725 end subroutine obs_da_value_allreduce
7727 subroutine obs_da_value_partial_reduce_iter(obsda, iter, nstart, nobs, ensval, qc, eqv, tm, pm )
7730 integer,
intent(in) :: iter
7731 integer,
intent(in) :: nstart
7732 integer,
intent(in) :: nobs
7733 real(RP),
intent(in) :: ensval(nobs)
7734 integer,
intent(in) :: qc(nobs)
7735 real(RP),
intent(in) :: eqv(nobs)
7736 real(RP),
intent(in) :: tm(nobs)
7737 real(RP),
intent(in) :: pm(nobs)
7744 nend = nstart + nobs - 1
7747 obsda%ensval(iter,nstart:nend) = ensval
7748 obsda%eqv(iter,nstart:nend) = eqv
7751 obsda%qc(nstart:nend) = max(obsda%qc(nstart:nend), qc)
7753 obsda%tm(nstart:nend) = obsda%tm(nstart:nend) + tm
7754 obsda%pm(nstart:nend) = obsda%pm(nstart:nend) + pm
7757 end subroutine obs_da_value_partial_reduce_iter
7762 subroutine read_ens_mpi_addiinfl(v3d, v2d)
7764 real(RP),
intent(out) :: v3d(:,:,:,:)
7765 real(RP),
intent(out) :: v2d(:,:,:)
7767 character(len=H_LONG) :: filename
7768 real(RP) :: v3dg(nlev,nlon,nlat,nv3d)
7769 real(RP) :: v2dg(nlon,nlat,nv2d)
7770 integer :: it, im, mstart, mend
7795 end subroutine read_ens_mpi_addiinfl
7800 subroutine monit_obs( OBS_IN_NUM, OBS_IN_FORMAT, v3dg, v2dg, nobs, bias, rmse, monit_type, use_key, step )
7807 integer,
intent(in) :: OBS_IN_NUM
7808 character(len=H_LONG),
intent(in) :: OBS_IN_FORMAT(:)
7810 real(RP),
intent(in) :: v3dg(nlev,nlon,nlat,nv3d)
7811 real(RP),
intent(in) :: v2dg(nlon,nlat,nv2d)
7812 integer,
intent(out) :: nobs(nid_obs)
7813 real(RP),
intent(out) :: bias(nid_obs)
7814 real(RP),
intent(out) :: rmse(nid_obs)
7815 logical,
intent(out) :: monit_type(nid_obs)
7816 logical,
intent(in) :: use_key
7817 integer,
intent(in) :: step
7819 real(RP) :: v3dgh(nlevh,nlonh,nlath,nv3dd)
7820 real(RP) :: v2dgh(nlonh,nlath,nv2dd)
7823 integer :: iset,iidx
7824 real(RP) :: ril,rjl,rk,rkz
7826 real(RP),
allocatable :: oelm(:)
7827 real(RP),
allocatable :: ohx(:)
7828 integer,
allocatable :: oqc(:)
7830 call state_to_history(v3dg, v2dg, v3dgh, v2dgh)
7833 nnobs = obsda_sort%nobs_in_key
7835 nnobs = obsda_sort%nobs
7838 allocate (oelm(nnobs))
7839 allocate (ohx(nnobs))
7840 allocate (oqc(nnobs))
7844 allocate (obsdep_set(obsdep_nobs))
7845 allocate (obsdep_idx(obsdep_nobs))
7846 allocate (obsdep_qc(obsdep_nobs))
7847 allocate (obsdep_omb(obsdep_nobs))
7848 allocate (obsdep_oma(obsdep_nobs))
7856 nn = obsda_sort%key(n)
7861 iset = obsda_sort%set(nn)
7862 iidx = obsda_sort%idx(nn)
7865 obsdep_set(n) = iset
7866 obsdep_idx(n) = iidx
7869 oelm(n) = obs(iset)%elm(iidx)
7870 call rij_g2l(
prc_myrank, obs(iset)%ri(iidx), obs(iset)%rj(iidx), ril, rjl)
7872 if (departure_stat_t_range <= 0.0d0 .or. &
7873 abs(obs(iset)%dif(iidx)) <= departure_stat_t_range)
then
7877 select case (obs_in_format(iset))
7881 call phys2ijk(v3dgh(:,:,:,iv3dd_p),obs(iset)%elm(iidx), &
7882 ril,rjl,obs(iset)%lev(iidx),rk,oqc(n))
7883 if (oqc(n) == iqc_good)
then
7884 call trans_xtoy(obs(iset)%elm(iidx),ril,rjl,rk, &
7885 obs(iset)%lon(iidx),obs(iset)%lat(iidx), &
7886 v3dgh,v2dgh,ohx(n),oqc(n),stggrd=1)
7889 case (
'RADAR',
'PAWR_TOSHIBA',
'MP_PAWR_TOSHIBA',
'PAWR_JRC',
'HIMAWARI8' )
7891 if (departure_stat_radar)
then
7892 call phys2ijkz(v3dgh(:,:,:,iv3dd_hgt),ril,rjl,obs(iset)%lev(iidx),rkz,oqc(n))
7893 if (oqc(n) == iqc_good)
then
7894 call trans_xtoy_radar(obs(iset)%elm(iidx),obs(iset)%meta(1), &
7895 obs(iset)%meta(2),obs(iset)%meta(3),ril,rjl,rkz, &
7896 obs(iset)%lon(iidx),obs(iset)%lat(iidx), &
7897 obs(iset)%lev(iidx),v3dgh,v2dgh,ohx(n),oqc(n),stggrd=1)
7898 if (oqc(n) == iqc_ref_low) oqc(n) = iqc_good
7904 if (oqc(n) == iqc_good)
then
7905 ohx(n) = obs(iset)%dat(iidx) - ohx(n)
7911 obsdep_qc(n) = oqc(n)
7912 obsdep_omb(n) = ohx(n)
7913 else if (step == 2)
then
7914 if (obsdep_qc(n) == iqc_good)
then
7915 obsdep_qc(n) = oqc(n)
7917 obsdep_oma(n) = ohx(n)
7926 call monit_dep(nnobs,oelm,ohx,oqc,nobs,bias,rmse)
7928 monit_type = .false.
7929 monit_type(uid_obs(id_u_obs)) = .true.
7930 monit_type(uid_obs(id_v_obs)) = .true.
7931 monit_type(uid_obs(id_t_obs)) = .true.
7932 monit_type(uid_obs(id_tv_obs)) = .true.
7933 monit_type(uid_obs(id_q_obs)) = .true.
7934 monit_type(uid_obs(id_rh_obs)) = .true.
7935 monit_type(uid_obs(id_ps_obs)) = .true.
7936 if (departure_stat_radar)
then
7937 monit_type(uid_obs(id_radar_ref_obs)) = .true.
7938 monit_type(uid_obs(id_radar_ref_zero_obs)) = .true.
7939 monit_type(uid_obs(id_radar_vr_obs)) = .true.
7947 end subroutine monit_obs
7952 subroutine monit_obs_mpi( OBS_IN_NUM, OBS_IN_FORMAT, v3dg, v2dg, monit_step, timelabel )
7957 integer,
intent(in) :: OBS_IN_NUM
7958 character(len=H_LONG),
intent(in) :: OBS_IN_FORMAT(:)
7960 real(RP),
intent(in) :: v3dg(nlev,nlon,nlat,nv3d)
7961 real(RP),
intent(in) :: v2dg(nlon,nlat,nv2d)
7962 integer,
intent(in) :: monit_step
7963 character(15),
intent(in),
optional :: timelabel
7965 integer :: nobs(nid_obs)
7966 integer :: nobs_g(nid_obs)
7967 real(RP) :: bias(nid_obs)
7968 real(RP) :: bias_g(nid_obs)
7969 real(RP) :: rmse(nid_obs)
7970 real(RP) :: rmse_g(nid_obs)
7971 logical :: monit_type(nid_obs)
7972 integer :: obsdep_g_nobs
7973 integer,
allocatable :: obsdep_g_set(:)
7974 integer,
allocatable :: obsdep_g_idx(:)
7975 integer,
allocatable :: obsdep_g_qc(:)
7976 real(RP),
allocatable :: obsdep_g_omb(:)
7977 real(RP),
allocatable :: obsdep_g_oma(:)
7979 integer :: cntr(NPRC_LCL)
7980 integer :: dspr(NPRC_LCL)
7981 integer :: i, ip, ierr
7983 if( rank_ens == 0 )
then
7984 call monit_obs( obs_in_num, obs_in_format, v3dg, v2dg, nobs, bias, rmse, monit_type, .true., monit_step )
7987 if (monit_type(i))
then
7989 if (nobs(i) == 0)
then
7993 bias_g(i) = bias(i) * real( nobs(i), kind=rp )
7994 rmse_g(i) = rmse(i) * rmse(i) * real( nobs(i), kind=rp )
7999 if( nprc_lcl > 1 )
then
8000 call mpi_allreduce(mpi_in_place, nobs_g, nid_obs, mpi_integer, mpi_sum, comm_lcl, ierr)
8001 call mpi_allreduce(mpi_in_place, bias_g, nid_obs, datatype, mpi_sum, comm_lcl, ierr)
8002 call mpi_allreduce(mpi_in_place, rmse_g, nid_obs, datatype, mpi_sum, comm_lcl, ierr)
8006 if (monit_type(i))
then
8007 if (nobs_g(i) == 0)
then
8011 bias_g(i) = bias_g(i) / real(nobs_g(i),kind=rp)
8012 rmse_g(i) = sqrt(rmse_g(i) / real(nobs_g(i),kind=rp))
8021 if (monit_step == 2)
then
8024 cntr(rank_lcl+1) = cnts
8025 call mpi_allreduce( mpi_in_place, cntr, nprc_lcl, mpi_integer, mpi_sum, comm_lcl, ierr )
8027 do ip = 1, nprc_lcl-1
8028 dspr(ip+1) = dspr(ip) + cntr(ip)
8031 obsdep_g_nobs = dspr(nprc_lcl) + cntr(nprc_lcl)
8032 allocate (obsdep_g_set(obsdep_g_nobs))
8033 allocate (obsdep_g_idx(obsdep_g_nobs))
8034 allocate (obsdep_g_qc(obsdep_g_nobs))
8035 allocate (obsdep_g_omb(obsdep_g_nobs))
8036 allocate (obsdep_g_oma(obsdep_g_nobs))
8038 if (obsdep_g_nobs > 0)
then
8039 call mpi_gatherv(obsdep_set, cnts, mpi_integer, obsdep_g_set, cntr, dspr, mpi_integer, 0, comm_lcl, ierr)
8040 call mpi_gatherv(obsdep_idx, cnts, mpi_integer, obsdep_g_idx, cntr, dspr, mpi_integer, 0, comm_lcl, ierr)
8041 call mpi_gatherv(obsdep_qc, cnts, mpi_integer, obsdep_g_qc, cntr, dspr, mpi_integer, 0, comm_lcl, ierr)
8042 call mpi_gatherv(obsdep_omb, cnts, datatype, obsdep_g_omb, cntr, dspr, datatype, 0, comm_lcl, ierr)
8043 call mpi_gatherv(obsdep_oma, cnts, datatype, obsdep_g_oma, cntr, dspr, datatype, 0, comm_lcl, ierr)
8046 deallocate (obsdep_g_set)
8047 deallocate (obsdep_g_idx)
8048 deallocate (obsdep_g_qc )
8049 deallocate (obsdep_g_omb)
8050 deallocate (obsdep_g_oma)
8054 if (monit_step == 2)
then
8055 deallocate (obsdep_set)
8056 deallocate (obsdep_idx)
8057 deallocate (obsdep_qc )
8058 deallocate (obsdep_omb)
8059 deallocate (obsdep_oma)
8063 call mpi_bcast(nobs, nid_obs, mpi_integer, 0, comm_ens, ierr)
8064 call mpi_bcast(bias, nid_obs, datatype, 0, comm_ens, ierr)
8065 call mpi_bcast(rmse, nid_obs, datatype, 0, comm_ens, ierr)
8066 call mpi_bcast(nobs_g, nid_obs, mpi_integer, 0, comm_ens, ierr)
8067 call mpi_bcast(bias_g, nid_obs, datatype, 0, comm_ens, ierr)
8068 call mpi_bcast(rmse_g, nid_obs, datatype, 0, comm_ens, ierr)
8069 call mpi_bcast(monit_type, nid_obs, mpi_logical, 0, comm_ens, ierr)
8071 if( monit_step == 1 )
then
8072 log_info(
"LETKF_debug",
'(1x,A)')
'OBSERVATIONAL DEPARTURE STATISTICS [GUESS] (IN THIS SUBDOMAIN):'
8073 else if (monit_step == 2)
then
8074 log_info(
"LETKF_debug",
'(1x,A)')
'OBSERVATIONAL DEPARTURE STATISTICS [ANALYSIS] (IN THIS SUBDOMAIN):'
8076 call monit_print(nobs, bias, rmse, monit_type)
8078 if (monit_step == 1)
then
8079 log_info(
"LETKF_debug",
'(1x,A)')
'OBSERVATIONAL DEPARTURE STATISTICS [GUESS] (GLOBAL):'
8080 else if (monit_step == 2)
then
8081 log_info(
"LETKF_debug",
'(1x,A)')
'OBSERVATIONAL DEPARTURE STATISTICS [ANALYSIS] (GLOBAL):'
8083 call monit_print(nobs_g, bias_g, rmse_g, monit_type)
8086 end subroutine monit_obs_mpi
8091 SUBROUTINE monit_dep(nn,elm,dep,qc,nobs,bias,rmse)
8095 INTEGER,
INTENT(IN) :: nn
8096 REAL(RP),
INTENT(IN) :: elm(nn)
8097 REAL(RP),
INTENT(IN) :: dep(nn)
8098 INTEGER,
INTENT(IN) :: qc(nn)
8099 INTEGER,
INTENT(OUT) :: nobs(nid_obs)
8100 REAL(RP),
INTENT(OUT) :: bias(nid_obs)
8101 REAL(RP),
INTENT(OUT) :: rmse(nid_obs)
8109 IF(qc(n) /= iqc_good) cycle
8112 if (ielm == id_tv_obs)
then
8115 if (ielm == id_radar_ref_zero_obs)
then
8116 ielm = id_radar_ref_obs
8120 nobs(i) = nobs(i) + 1
8121 bias(i) = bias(i) + dep(n)
8122 rmse(i) = rmse(i) + dep(n)**2
8126 IF(nobs(i) == 0)
THEN
8130 bias(i) = bias(i) / real(nobs(i),kind=rp)
8131 rmse(i) = sqrt(rmse(i) / real(nobs(i),kind=rp))
8136 END SUBROUTINE monit_dep
8141 SUBROUTINE monit_print(nobs,bias,rmse,monit_type)
8143 INTEGER,
INTENT(IN) :: nobs(nid_obs)
8144 REAL(RP),
INTENT(IN) :: bias(nid_obs)
8145 REAL(RP),
INTENT(IN) :: rmse(nid_obs)
8146 LOGICAL,
INTENT(IN),
OPTIONAL :: monit_type(nid_obs)
8148 character(12) :: var_show(nid_obs)
8149 character(12) :: nobs_show(nid_obs)
8150 character(12) :: bias_show(nid_obs)
8151 character(12) :: rmse_show(nid_obs)
8154 character(4) :: nstr
8156 logical :: monit_type_(nid_obs)
8158 monit_type_ = .true.
8159 if (
present(monit_type)) monit_type_ = monit_type
8163 if (monit_type_(i) .and. i /= uid_obs(id_tv_obs) .and. i /= uid_obs(id_radar_ref_zero_obs))
then
8165 write(var_show(n),
'(A12)') obelmlist(i)
8166 write(nobs_show(n),
'(I12)') nobs(i)
8167 if (nobs(i) > 0)
then
8168 write(bias_show(n),
'(ES12.3)') bias(i)
8169 write(rmse_show(n),
'(ES12.3)') rmse(i)
8171 write(bias_show(n),
'(A12)')
'N/A'
8172 write(rmse_show(n),
'(A12)')
'N/A'
8176 write(nstr,
'(I4)') n
8178 log_info(
"LETKF_debug",
'(1x,A,' // trim(nstr) //
"('============'))")
'======'
8179 log_info(
"LETKF_debug",
'(1x,A,' // trim(nstr) //
'A)' )
' ', var_show(1:n)
8180 log_info(
"LETKF_debug",
'(1x,A,' // trim(nstr) //
"('------------'))")
'------'
8181 log_info(
"LETKF_debug",
'(1x,A,' // trim(nstr) //
'A)' )
'BIAS ', bias_show(1:n)
8182 log_info(
"LETKF_debug",
'(1x,A,' // trim(nstr) //
'A)' )
'RMSE ', rmse_show(1:n)
8183 log_info(
"LETKF_debug",
'(1x,A,' // trim(nstr) //
'A)' )
'NUMBER', nobs_show(1:n)
8184 log_info(
"LETKF_debug",
'(1x,A,' // trim(nstr) //
"('============'))")
'======'
8187 END SUBROUTINE monit_print
8199 subroutine state_to_history(v3dg, v2dg, v3dgh, v2dgh)
8212 atmos_saturation_psat_all
8215 real(RP),
intent(in) :: v3dg(nlev,nlon,nlat,nv3d)
8216 real(RP),
intent(in) :: v2dg(nlon,nlat,nv2d)
8217 real(RP),
intent(out) :: v3dgh(nlevh,nlonh,nlath,nv3dd)
8218 real(RP),
intent(out) :: v2dgh(nlonh,nlath,nv2dd)
8219 real(RP) :: v3dgh_RP(nlevh,nlonh,nlath,nv3dd)
8220 real(RP) :: v2dgh_RP(nlonh,nlath,nv2dd)
8222 integer :: i, j, k, iv3d, iv2d
8224 real(RP) :: utmp, vtmp
8225 real(RP) :: qdry, Rtot
8226 real(RP) :: psat(nlevh,nlonh,nlath)
8230 v3dgh_rp(:,:,:,:) = undef
8231 v2dgh_rp(:,:,:) = undef
8233 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_u) = v3dg(:,:,:,iv3d_u)
8234 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_v) = v3dg(:,:,:,iv3d_v)
8235 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_w) = v3dg(:,:,:,iv3d_w)
8236 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_t) = v3dg(:,:,:,iv3d_t)
8237 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_p) = v3dg(:,:,:,iv3d_p)
8238 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_q) = v3dg(:,:,:,iv3d_q)
8239 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_qc) = v3dg(:,:,:,iv3d_qc)
8240 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_qr) = v3dg(:,:,:,iv3d_qr)
8241 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_qi) = v3dg(:,:,:,iv3d_qi)
8242 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_qs) = v3dg(:,:,:,iv3d_qs)
8243 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_qg) = v3dg(:,:,:,iv3d_qg)
8247 do j = start_y, end_y
8248 do i = start_x, end_x
8249 do k = start_z, end_z
8250 utmp = v3dgh_rp(k,i,j,iv3d_u)
8251 vtmp = v3dgh_rp(k,i,j,iv3d_v)
8253 v3dgh_rp(k,i,j,iv3d_u) = utmp * rotc(i,j,1) - vtmp * rotc(i,j,2)
8254 v3dgh_rp(k,i,j,iv3d_v) = utmp * rotc(i,j,2) + vtmp * rotc(i,j,1)
8262 call atmos_saturation_psat_all( nlevh, start_z, end_z, &
8263 nlonh, start_x, end_x, &
8264 nlath, start_y, end_y, &
8265 v3dgh_rp(:,:,:,iv3dd_t), &
8269 do j = start_y, end_y
8270 do i = start_x, end_x
8271 do k = start_z, end_z
8273 do iv3d = iv3dd_q, iv3dd_qg
8274 qdry = qdry - v3dgh_rp(k,i,j,iv3d)
8276 rtot = rdry * qdry + rvap * v3dgh_rp(k,i,j,iv3dd_q)
8278 v3dgh_rp(k,i,j,iv3dd_rh) = v3dgh_rp(k,i,j,iv3dd_q) * v3dgh_rp(k,i,j,iv3dd_p) / psat(k,i,j) * rvap / rtot
8286 v3dgh_rp(start_z:end_z,start_x:end_x,start_y:end_y,iv3dd_hgt) = real_cz(start_z:end_z,start_x:end_x,start_y:end_y)
8291 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_topo) = real_cz(start_z,start_x:end_x,start_y:end_y)
8293 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_ps) = v3dg(1,1:nlon,1:nlat,iv3d_p)
8294 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_u10m) = v3dg(1,1:nlon,1:nlat,iv3d_u)
8295 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_v10m) = v3dg(1,1:nlon,1:nlat,iv3d_v)
8296 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_t2m) = v3dg(1,1:nlon,1:nlat,iv3d_t)
8297 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_q2m) = v3dg(1,1:nlon,1:nlat,iv3d_q)
8301 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_ps) = v3dg(1,:,:,iv3d_p)
8302 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_u10m) = v3dg(1,:,:,iv3d_u)
8303 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_v10m) = v3dg(1,:,:,iv3d_v)
8304 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_t2m) = v3dg(1,:,:,iv3d_t)
8305 v2dgh_rp(start_x:end_x,start_y:end_y,iv2dd_q2m) = v3dg(1,:,:,iv3d_q)
8310 do j = start_y, end_y
8311 do i = start_x, end_x
8312 v3dgh_rp( 1:start_z-1,i,j,iv3d) = v3dgh_rp(start_z,i,j,iv3d)
8313 v3dgh_rp(end_z+1:nlevh ,i,j,iv3d) = v3dgh_rp(end_z ,i,j,iv3d)
8322 call comm_vars8( v3dgh_rp(:,:,:,iv3d), iv3d )
8325 call comm_wait ( v3dgh_rp(:,:,:,iv3d), iv3d )
8329 call comm_vars8( v2dgh_rp(:,:,iv2d), iv2d )
8332 call comm_wait ( v2dgh_rp(:,:,iv2d), iv2d )
8335 v3dgh = real(v3dgh_rp, kind=rp)
8336 v2dgh = real(v2dgh_rp, kind=rp)
8339 end subroutine state_to_history