FV3 Bundle
m_diag_raob.f90
Go to the documentation of this file.
1 module m_diag_raob
2  use fckit_log_module, only: fckit_log
3  use ncd_kinds, only: i_kind,r_single,r_kind,r_double
4  use nc_diag_write_mod, only: nc_diag_init, nc_diag_metadata, nc_diag_write
5 
6  use nc_diag_read_mod, only: nc_diag_read_get_dim
7  use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_close
8  use nc_diag_read_mod, only: nc_diag_read_get_var
9  use nc_diag_read_mod, only: nc_diag_read_get_global_attr
10 
11  use type_distribution, only: random_distribution
12 
13  implicit none
14 
15  private
16  save
17 
18  public :: diag_raob_header
19  public :: diag_raob_mass !generic name for non-wind obs - obs w/ single obs associated
20  public :: diag_raob_wind ! name for wind obs - obs w/ two obs (u,v) associated
21 
22  public :: write_split_raob_diag_nc
23 
24  public :: read_raob_diag_nc_header
25  public :: read_raob_diag_nc_mass
26 
27  type diag_raob_header
28  character(3), dimension(:), allocatable :: obstype
29  integer(i_kind) :: n_obstype
30  integer(i_kind), dimension(:), allocatable :: n_observations
31  integer(i_kind) :: n_observations_mass
32  integer(i_kind) :: n_observations_wind
33  integer(i_kind) :: n_observations_total
34  integer(i_kind) :: date
35  end type diag_raob_header
36 
37  type diag_raob_mass
38  character(8) :: station_id
39  character(3) :: observation_class
40  real(r_kind) :: observation_type
41  real(r_kind) :: observation_subtype
42  real(r_kind) :: latitude
43  real(r_kind) :: longitude
44  real(r_kind) :: station_elevation
45  real(r_kind) :: pressure
46  real(r_kind) :: height
47  real(r_kind) :: time
48  real(r_kind) :: prep_qc_mark
49  real(r_kind) :: setup_qc_mark
50  real(r_kind) :: prep_use_flag
51  real(r_kind) :: analysis_use_flag
52  real(r_kind) :: nonlinear_qc_rel_wgt
53  real(r_kind) :: errinv_input
54  real(r_kind) :: errinv_adjust
55  real(r_kind) :: errinv_final
56  real(r_kind) :: observation
57  real(r_kind) :: obs_minus_forecast_adjusted
58  real(r_kind) :: obs_minus_forecast_unadjusted
59 
60  end type diag_raob_mass
61 
62  type diag_raob_wind
63  character(8) :: station_id
64  character(3) :: observation_class
65  real(r_kind) :: observation_type
66  real(r_kind) :: observation_subtype
67  real(r_kind) :: latitude
68  real(r_kind) :: longitude
69  real(r_kind) :: station_elevation
70  real(r_kind) :: pressure
71  real(r_kind) :: height
72  real(r_kind) :: time
73  real(r_kind) :: prep_qc_mark
74  real(r_kind) :: setup_qc_mark
75  real(r_kind) :: prep_use_flag
76  real(r_kind) :: analysis_use_flag
77  real(r_kind) :: nonlinear_qc_rel_wgt
78  real(r_kind) :: errinv_input
79  real(r_kind) :: errinv_adjust
80  real(r_kind) :: errinv_final
81  real(r_kind) :: u_observation
82  real(r_kind) :: u_obs_minus_forecast_adjusted
83  real(r_kind) :: u_obs_minus_forecast_unadjusted
84  real(r_kind) :: v_observation
85  real(r_kind) :: v_obs_minus_forecast_adjusted
86  real(r_kind) :: v_obs_minus_forecast_unadjusted
87  real(r_kind) :: wind_reduction_factor_at_10m
88  end type diag_raob_wind
89 
90  integer,parameter :: maxobstype=30
91  integer :: nobstype
92  integer :: raob_mass_type = 120 ! type for RAOB T and Q
93  integer :: raob_wind_type = 220 ! type for RAOB U and V
94  integer :: t_qcmark = 0 ! 0=tv; 1=tdry
95 
96 contains
97 
98  integer(i_kind) function get_obstype_index(obstype, obstypearr)
99 ! integer(i_kind) :: get_obstype_index
100  character(3),intent(in) :: obstype
101  character(3),intent(inout),dimension(*) :: obstypearr
102 
103  integer :: i, idx
104  logical :: matched
105 
106  matched = .false.
107 
108  if (nobstype .eq. 0) then
109  nobstype = 1
110  obstypearr(1) = obstype
111  idx = nobstype
112  print *,'obstype=',obstype,' set to index',idx
113  else
114  do i=1,nobstype
115  if (obstype .eq. obstypearr(i)) then
116  idx = i
117  matched = .true.
118  endif
119  enddo
120  if (.not. matched) then
121  nobstype = nobstype + 1
122  obstypearr(nobstype) = obstype
123  idx = nobstype
124  print *,'obstype=',obstype,' set to index',idx
125  endif
126  endif
127 
128 
129  get_obstype_index = idx
130 
131  end function get_obstype_index
132 
133 
134  subroutine write_split_raob_diag_nc(infn,raob_header, raob_mass, raob_wind, append_suffix)
135  character(120), intent(in) :: infn
136  type(diag_raob_header), intent(in) :: raob_header
137  type(diag_raob_mass),dimension(raob_header%n_Observations_Mass),intent(in) :: raob_mass
138  type(diag_raob_wind),dimension(raob_header%n_Observations_Wind),intent(in) :: raob_wind
139  logical, intent(in) :: append_suffix
140 
141  character(120) :: outfn
142  character(20) :: str, str2
143  integer :: strlen
144  integer :: i, itype
145 
146  do itype=1, raob_header%n_ObsType
147  str = raob_header%ObsType(itype)
148  if (.not. append_suffix) then
149  str2 = 'diag_raob_' // trim(adjustl(str))
150  outfn = replace_text(trim(infn),'diag_raob',str2)
151  strlen = len(trim(outfn))
152  outfn = outfn(1:strlen-3) // 'nc4'
153  else
154  outfn = trim(infn) // '.' // trim(adjustl(str)) // '.nc4'
155  endif
156 
157  print *,outfn
158 
159  call nc_diag_init(outfn)
160 
161  if (raob_header%ObsType(itype) .eq. ' uv') then
162  do i=1,raob_header%n_Observations_Wind
163  call nc_diag_metadata("Station_ID", raob_wind(i)%Station_ID )
164  call nc_diag_metadata("Observation_Class", raob_wind(i)%Observation_Class )
165  call nc_diag_metadata("Observation_Type", raob_wind(i)%Observation_Type )
166  call nc_diag_metadata("Observation_Subtype", raob_wind(i)%Observation_Subtype )
167  call nc_diag_metadata("Latitude", raob_wind(i)%Latitude )
168  call nc_diag_metadata("Longitude", raob_wind(i)%Longitude )
169  call nc_diag_metadata("Station_Elevation", raob_wind(i)%Station_Elevation )
170  call nc_diag_metadata("Pressure", raob_wind(i)%Pressure )
171  call nc_diag_metadata("Height", raob_wind(i)%Height )
172  call nc_diag_metadata("Time", raob_wind(i)%Time )
173  call nc_diag_metadata("Prep_QC_Mark", raob_wind(i)%Prep_QC_Mark )
174  call nc_diag_metadata("Setup_QC_Mark", raob_wind(i)%Setup_QC_Mark )
175  call nc_diag_metadata("Prep_Use_Flag", raob_wind(i)%Prep_Use_Flag )
176  call nc_diag_metadata("Analysis_Use_Flag", raob_wind(i)%Analysis_Use_Flag )
177  call nc_diag_metadata("Nonlinear_QC_Rel_Wgt", raob_wind(i)%Nonlinear_QC_Rel_Wgt )
178  call nc_diag_metadata("Errinv_Input", raob_wind(i)%Errinv_Input )
179  call nc_diag_metadata("Errinv_Adjust", raob_wind(i)%Errinv_Adjust )
180  call nc_diag_metadata("Errinv_Final", raob_wind(i)%Errinv_Final )
181  call nc_diag_metadata("u_Observation", raob_wind(i)%u_Observation )
182  call nc_diag_metadata("u_Obs_Minus_Forecast_adjusted", raob_wind(i)%u_Obs_Minus_Forecast_adjusted )
183  call nc_diag_metadata("u_Obs_Minus_Forecast_unadjusted", raob_wind(i)%u_Obs_Minus_Forecast_unadjusted )
184  call nc_diag_metadata("v_Observation", raob_wind(i)%v_Observation )
185  call nc_diag_metadata("v_Obs_Minus_Forecast_adjusted", raob_wind(i)%v_Obs_Minus_Forecast_adjusted )
186  call nc_diag_metadata("v_Obs_Minus_Forecast_unadjusted", raob_wind(i)%v_Obs_Minus_Forecast_unadjusted )
187  call nc_diag_metadata("Wind_Reduction_Factor_at_10m", raob_wind(i)%Wind_Reduction_Factor_at_10m )
188  enddo
189  else
190  do i=1,raob_header%n_Observations_Mass
191  if (raob_mass(i)%Observation_Class .eq. raob_header%ObsType(itype) ) then
192  call nc_diag_metadata("Station_ID", raob_mass(i)%Station_ID )
193  call nc_diag_metadata("Observation_Class", raob_mass(i)%Observation_Class )
194  call nc_diag_metadata("Observation_Type", raob_mass(i)%Observation_Type )
195  call nc_diag_metadata("Observation_Subtype", raob_mass(i)%Observation_Subtype )
196  call nc_diag_metadata("Latitude", raob_mass(i)%Latitude )
197  call nc_diag_metadata("Longitude", raob_mass(i)%Longitude )
198  call nc_diag_metadata("Station_Elevation", raob_mass(i)%Station_Elevation )
199  call nc_diag_metadata("Pressure", raob_mass(i)%Pressure )
200  call nc_diag_metadata("Height", raob_mass(i)%Height )
201  call nc_diag_metadata("Time", raob_mass(i)%Time )
202  call nc_diag_metadata("Prep_QC_Mark", raob_mass(i)%Prep_QC_Mark )
203  call nc_diag_metadata("Setup_QC_Mark", raob_mass(i)%Setup_QC_Mark )
204  call nc_diag_metadata("Prep_Use_Flag", raob_mass(i)%Prep_Use_Flag )
205  call nc_diag_metadata("Analysis_Use_Flag", raob_mass(i)%Analysis_Use_Flag )
206  call nc_diag_metadata("Nonlinear_QC_Rel_Wgt", raob_mass(i)%Nonlinear_QC_Rel_Wgt )
207  call nc_diag_metadata("Errinv_Input", raob_mass(i)%Errinv_Input )
208  call nc_diag_metadata("Errinv_Adjust", raob_mass(i)%Errinv_Adjust )
209  call nc_diag_metadata("Errinv_Final", raob_mass(i)%Errinv_Final )
210  call nc_diag_metadata("Observation", raob_mass(i)%Observation )
211  call nc_diag_metadata("Obs_Minus_Forecast_adjusted", raob_mass(i)%Obs_Minus_Forecast_adjusted )
212  call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", raob_mass(i)%Obs_Minus_Forecast_unadjusted )
213 
214  endif
215  enddo
216  endif
217 
218  call nc_diag_write
219 
220  enddo
221  end subroutine write_split_raob_diag_nc
222 
223  subroutine read_raob_diag_nc_header(infn,raob_header)
224  character(len=*), intent(in) :: infn
225  type(diag_raob_header), intent(inout) :: raob_header
226 
227  integer(i_kind) :: fid,nobs, gnobs
228  type(random_distribution) :: distribution
229 
230  nobs=0
231  call nc_diag_read_init(infn,fid)
232  gnobs = nc_diag_read_get_dim(fid,'nobs')
233  distribution=random_distribution(gnobs)
234  nobs=distribution%nobs_pe()
235  call nc_diag_read_get_global_attr(fid,"date_time", raob_header%date )
236  raob_header%n_Observations_Mass = nobs
237  call nc_diag_read_close(infn)
238 
239  end subroutine read_raob_diag_nc_header
240 
241  subroutine read_raob_diag_nc_mass(infn, raob_header, raob_mass, ierr)
242  character(len=*), intent(in) :: infn
243  type(diag_raob_header), intent(inout) :: raob_header
244  type(diag_raob_mass),pointer, intent(inout) :: raob_mass(:)
245  integer, intent(out) :: ierr
246 
247  character(20) :: str, str2
248  integer(i_kind) :: ii, ic, fid, nobs, nraob, ncount(1)
249  character(len=8),allocatable, dimension(:) :: c_var
250  integer(i_kind), allocatable, dimension(:) :: i_var
251  real(r_kind), allocatable, dimension(:) :: r_var
252 
253  type(diag_raob_mass), pointer :: rtmp_mass(:)
254 
255  type(random_distribution) :: distribution
256 
257  ierr=0
258  call nc_diag_read_init(infn,fid)
259 
260  nobs = nc_diag_read_get_dim(fid,'nobs')
261  allocate(rtmp_mass(nobs))
262  allocate(c_var(nobs))
263  allocate(i_var(nobs))
264  allocate(r_var(nobs))
265 
266 ! if (raob_mass(i)%Observation_Class .eq. raob_header%ObsType(itype) ) then
267 ! call nc_diag_read_get_var(fid,"Station_ID", c_var ); rtmp_mass(:)%Station_ID = c_var
268 ! call nc_diag_read_get_var(fid,"Observation_Class", c_var ); rtmp_mass(:)%Observation_Class = c_var
269  call nc_diag_read_get_var(fid,"Observation_Type", i_var ); rtmp_mass(:)%Observation_Type = i_var
270  call nc_diag_read_get_var(fid,"Observation_Subtype", i_var ); rtmp_mass(:)%Observation_Subtype = i_var
271  call nc_diag_read_get_var(fid,"Latitude", r_var ); rtmp_mass(:)%Latitude = r_var
272  call nc_diag_read_get_var(fid,"Longitude", r_var ); rtmp_mass(:)%Longitude = r_var
273  call nc_diag_read_get_var(fid,"Pressure", r_var ); rtmp_mass(:)%Pressure = r_var
274  call nc_diag_read_get_var(fid,"Height", r_var ); rtmp_mass(:)%Height = r_var
275  call nc_diag_read_get_var(fid,"Time", r_var ); rtmp_mass(:)%Time = r_var
276  call nc_diag_read_get_var(fid,"Prep_QC_Mark", r_var ); rtmp_mass(:)%Prep_QC_Mark = r_var
277  call nc_diag_read_get_var(fid,"Setup_QC_Mark", r_var ); rtmp_mass(:)%Setup_QC_Mark = r_var
278  call nc_diag_read_get_var(fid,"Prep_Use_Flag", r_var ); rtmp_mass(:)%Prep_Use_Flag = r_var
279  call nc_diag_read_get_var(fid,"Analysis_Use_Flag", r_var ); rtmp_mass(:)%Analysis_Use_Flag = r_var
280  call nc_diag_read_get_var(fid,"Nonlinear_QC_Rel_Wgt", r_var ); rtmp_mass(:)%Nonlinear_QC_Rel_Wgt = r_var
281  call nc_diag_read_get_var(fid,"Errinv_Input", r_var ); rtmp_mass(:)%Errinv_Input = r_var
282  call nc_diag_read_get_var(fid,"Errinv_Adjust", r_var ); rtmp_mass(:)%Errinv_Adjust = r_var
283  call nc_diag_read_get_var(fid,"Errinv_Final", r_var ); rtmp_mass(:)%Errinv_Final = r_var
284  call nc_diag_read_get_var(fid,"Observation", r_var ); rtmp_mass(:)%Observation = r_var
285  call nc_diag_read_get_var(fid,"Obs_Minus_Forecast_adjusted", r_var ); rtmp_mass(:)%Obs_Minus_Forecast_adjusted = r_var
286  call nc_diag_read_get_var(fid,"Obs_Minus_Forecast_unadjusted",r_var ); rtmp_mass(:)%Obs_Minus_Forecast_unadjusted = r_var
287 ! endif
288  deallocate(r_var)
289  deallocate(i_var)
290  deallocate(c_var)
291 
292  ic=0
293  do ii=1,nobs
294  if(rtmp_mass(ii)%Observation_Type==raob_mass_type.and.&
295  rtmp_mass(ii)%Setup_QC_Mark==t_qcmark) then
296  ic=ic+1
297  endif
298  enddo
299  !> randomly distribution among PEs
300  distribution=random_distribution(ic)
301  nraob=distribution%nobs_pe()
302  ic=0
303  call distribution%reset()
304  do ii=1,nobs
305  if(rtmp_mass(ii)%Observation_Type==raob_mass_type.and.&
306  rtmp_mass(ii)%Setup_QC_Mark==t_qcmark) then
307  ic=ic+1
308  if (distribution%received(ic)) &
309  distribution%indx(distribution%nobs_pe())=ii
310  endif
311  enddo
312  print *, 'found this many raob ', nraob, ' on PE ', distribution%mype()
313  if(distribution%nobs_pe() /= size(distribution%indx)) then
314  print *, 'error determining Raob, inconsistent nraob, nobs_pe=', nraob, distribution%nobs_pe()
315  deallocate(raob_mass)
316  ierr = 99
317  return
318  endif
319 
320  if(associated(raob_mass)) deallocate(raob_mass)
321  allocate (raob_mass(nraob))
322  raob_header%n_Observations_Mass = nraob
323 
324 ! raob_mass(:)%Station_ID = rtmp_mass(distribution%indx)
325 ! raob_mass(:)%Observation_Class = c_var
326  raob_mass(:)%Observation_Type = rtmp_mass(distribution%indx)%Observation_Type
327  raob_mass(:)%Observation_Subtype = rtmp_mass(distribution%indx)%Observation_Subtype
328  raob_mass(:)%Latitude = rtmp_mass(distribution%indx)%Latitude
329  raob_mass(:)%Longitude = rtmp_mass(distribution%indx)%Longitude
330  raob_mass(:)%Pressure = rtmp_mass(distribution%indx)%Pressure
331  raob_mass(:)%Height = rtmp_mass(distribution%indx)%Height
332  raob_mass(:)%Time = rtmp_mass(distribution%indx)%Time
333  raob_mass(:)%Prep_QC_Mark = rtmp_mass(distribution%indx)%Prep_QC_Mark
334  raob_mass(:)%Setup_QC_Mark = rtmp_mass(distribution%indx)%Setup_QC_Mark
335  raob_mass(:)%Prep_Use_Flag = rtmp_mass(distribution%indx)%Prep_Use_Flag
336  raob_mass(:)%Analysis_Use_Flag = rtmp_mass(distribution%indx)%Analysis_Use_Flag
337  raob_mass(:)%Nonlinear_QC_Rel_Wgt= rtmp_mass(distribution%indx)%Nonlinear_QC_Rel_Wgt
338  raob_mass(:)%Errinv_Input = rtmp_mass(distribution%indx)%Errinv_Input
339  raob_mass(:)%Errinv_Adjust = rtmp_mass(distribution%indx)%Errinv_Adjust
340  raob_mass(:)%Errinv_Final = rtmp_mass(distribution%indx)%Errinv_Final
341  raob_mass(:)%Observation = rtmp_mass(distribution%indx)%Observation
342  raob_mass(:)%Obs_Minus_Forecast_adjusted = rtmp_mass(distribution%indx)%Obs_Minus_Forecast_adjusted
343  raob_mass(:)%Obs_Minus_Forecast_unadjusted = rtmp_mass(distribution%indx)%Obs_Minus_Forecast_unadjusted
344 
345  deallocate(rtmp_mass)
346 
347  call nc_diag_read_close(infn)
348  end subroutine read_raob_diag_nc_mass
349 
350  function replace_text (s,text,rep) result(outs)
351  character(*) :: s,text,rep
352  character(len(s)+100) :: outs ! provide outs with extra 100 char len
353  integer :: i, nt, nr
354 
355  outs = s ; nt = len_trim(text) ; nr = len_trim(rep)
356  i = index(outs,text(:nt)) !; IF (i == 0) EXIT
357  outs = outs(:i-1) // rep(:nr) // outs(i+nt:)
358  end function replace_text
359 
360 
361 end module m_diag_raob
m_diag_raob::raob_mass_type
integer raob_mass_type
Definition: m_diag_raob.f90:92
m_diag_raob::write_split_raob_diag_nc
subroutine, public write_split_raob_diag_nc(infn, raob_header, raob_mass, raob_wind, append_suffix)
Definition: m_diag_raob.f90:135
crtm_parameters::strlen
integer, parameter, public strlen
Definition: CRTM_Parameters.f90:116
ncd_kinds::i_kind
integer, parameter, public i_kind
Definition: ncd_kinds.F90:71
nc_diag_read_mod
Definition: nc_diag_read_mod.F90:1
m_diag_raob::maxobstype
integer, parameter maxobstype
Definition: m_diag_raob.f90:90
m_diag_raob::diag_raob_mass
Definition: m_diag_raob.f90:37
nc_diag_write_mod::nc_diag_init
subroutine nc_diag_init(filename, append)
Definition: nc_diag_write_mod.F90:292
type_distribution
Definition: m_distribution.f90:1
m_diag_raob::read_raob_diag_nc_mass
subroutine, public read_raob_diag_nc_mass(infn, raob_header, raob_mass, ierr)
Definition: m_diag_raob.f90:242
m_diag_raob::diag_raob_header
Definition: m_diag_raob.f90:27
m_diag_raob::get_obstype_index
integer(i_kind) function get_obstype_index(obstype, obstypearr)
Definition: m_diag_raob.f90:99
m_diag_raob::diag_raob_wind
Definition: m_diag_raob.f90:62
m_diag_raob::nobstype
integer nobstype
Definition: m_diag_raob.f90:91
nc_diag_write_mod
Definition: nc_diag_write_mod.F90:19
m_diag_raob::raob_wind_type
integer raob_wind_type
Definition: m_diag_raob.f90:93
m_diag_raob
Definition: m_diag_raob.f90:1
nc_diag_read_mod::nc_diag_read_close
subroutine nc_diag_read_close(filename, file_ncdr_id, from_pop)
Definition: nc_diag_read_mod.F90:225
ncd_kinds
Definition: ncd_kinds.F90:1
type_distribution::random_distribution
Definition: m_distribution.f90:10
m_diag_raob::t_qcmark
integer t_qcmark
Definition: m_diag_raob.f90:94
ncd_kinds::r_double
integer, parameter, public r_double
Definition: ncd_kinds.F90:80
nc_diag_write_mod::nc_diag_write
subroutine nc_diag_write
Definition: nc_diag_write_mod.F90:505
ncd_kinds::r_single
integer, parameter, public r_single
Definition: ncd_kinds.F90:79
m_diag_raob::replace_text
character(len(s)+100) function replace_text(s, text, rep)
Definition: m_diag_raob.f90:351
ncd_kinds::r_kind
integer, parameter, public r_kind
Definition: ncd_kinds.F90:108
nc_diag_read_mod::nc_diag_read_init
subroutine nc_diag_read_init(filename, file_ncdr_id, from_push)
Definition: nc_diag_read_mod.F90:176
m_diag_raob::read_raob_diag_nc_header
subroutine, public read_raob_diag_nc_header(infn, raob_header)
Definition: m_diag_raob.f90:224