da_scan_obs_ssmi.inc
References to this file elsewhere.
1 subroutine da_scan_obs_ssmi (iv, filename)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Read the header of a SSMI 2.0 GTS observation file
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type(iv_type), intent(inout) :: iv
10 character(len=*), intent(in) :: filename
11
12 character (len = 10) :: fmt_name
13 character (len = 160) :: fmt_info, fmt_loc
14 character (len = 120) :: char_ned
15
16 integer :: iost, fm,iunit
17
18 type (model_loc_type) :: loc
19 type (info_type) :: info
20 type (field_type) :: speed, tpw
21
22 type (field_type) :: tb19v, tb19h, tb22v
23 type (field_type) :: tb37v, tb37h, tb85v, tb85h
24
25 logical :: isfilter,ipass
26 logical :: outside
27 integer :: irain, iprecip
28
29 if (trace_use) call da_trace_entry("da_scan_obs_ssmi")
30
31 isfilter = .true. ! filter out rain points
32 irain = 0
33
34 ! open FILE
35 call da_get_unit(iunit)
36 open(unit = iunit, &
37 FILE = trim(filename), &
38 FORM = 'FORMATTED', &
39 ACCESS = 'SEQUENTIAL', &
40 iostat = iost, &
41 STATUS = 'OLD')
42
43 if (iost /= 0) then
44 call da_warning(__FILE__,__LINE__, (/"Cannot open SSMI file "//trim(filename)/))
45 call da_free_unit(iunit)
46 return
47 end if
48
49 rewind (unit = iunit)
50
51 ! 2. read header
52 ! ===============
53
54 ! 2.1 read first line
55 ! ---------------
56
57 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
58 if (iost /= 0) then
59 call da_error(__FILE__,__LINE__, (/"Cannot read SSMI file "//trim(filename)/))
60 end if
61
62 ! 2.3 read number of reports
63 do
64 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
65 if (iost /= 0) then
66 call da_error(__FILE__,__LINE__, (/"Cannot read SSMI file "//trim(filename)/))
67 end if
68 if (char_ned(1:6) == 'NESTIX') exit
69 end do
70
71 do
72 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
73 if (char_ned(1:6) == 'INFO ') exit
74 end do
75
76 read (unit = iunit, fmt = '(A)', iostat = iost) char_ned
77
78 ! read formats
79 read (unit=iunit, fmt = '(A,1X,A)') &
80 fmt_name, fmt_info, &
81 fmt_name, fmt_loc
82
83 ! skip 1 line
84 read (unit=iunit, fmt = '(A)') fmt_name
85
86 ! loop over records
87 reports: do
88 ! read station general info
89 read (unit=iunit, fmt = fmt_info, iostat = iost) &
90 info % platform, &
91 info % date_char, &
92 info % name, &
93 info % levels, &
94 info % lat, &
95 info % lon, &
96 info % elv, &
97 info % id
98
99 read(unit=info % platform (4:6),fmt='(I3)') fm
100 if (iost /= 0) exit reports
101
102 select case(fm)
103 case (125) ;
104 ! read surface wind speed and precipitable water
105 read (unit=iunit, fmt = fmt_loc) speed, tpw
106 case (126) ;
107 read (unit=iunit, fmt = fmt_loc) tb19v, tb19h, tb22v, tb37v, tb37h, tb85v, tb85h
108 case default;
109 write(unit=message(1), fmt='(a, i6)') 'unsaved ssmi obs found, fm=', fm
110 write(unit=message(2), fmt='(a, 2f12.6)') 'info%(lon,lat)=', info%lon, info%lat
111 call da_warning(__FILE__,__LINE__,message(1:2))
112 end select
113
114 ! check if obs is in horizontal domain
115
116 ! Compute the model horizontal coordinate x, y
117 ! Check if obs is wihin horizontal domain
118
119 call da_llxy (info, loc, outside)
120
121 select case(fm)
122 case (125) ;
123 if (.not. use_ssmiretrievalobs .or. iv%info(ssmi_rv)%ntotal == max_ssmi_rv_input) cycle reports
124
125 ! Check if at least one field is present
126 if ((tpw % qc == missing) .AND. (speed % qc == missing)) then
127 cycle reports
128 end if
129 iv%info(ssmi_rv)%ntotal = iv%info(ssmi_rv)%ntotal + 1
130 if (outside) cycle reports
131
132 iv%info(ssmi_rv)%nlocal = iv%info(ssmi_rv)%nlocal + 1
133
134 case (126) ;
135 if (.not. use_ssmitbobs .or. iv%info(ssmi_tb)%ntotal == max_ssmi_tb_input) cycle reports
136
137 ! Check if at least one field is present
138
139 if ((tb19v % qc == missing) .AND. (tb19h % qc == missing) .AND. &
140 (tb22v % qc == missing) .AND. &
141 (tb37v % qc == missing) .AND. (tb37h % qc == missing) .AND. &
142 (tb85v % qc == missing) .AND. (tb85h % qc == missing)) then
143 cycle reports
144 end if
145
146 ! filter rain pixels
147 ! ====================================
148
149 if (isfilter) then
150 ipass = .false.
151 iprecip = 0
152 call filter(tb19v%inv, tb19h%inv, tb22v%inv, tb37v%inv, &
153 tb37h%inv, tb85v%inv, tb85h%inv, ipass, iprecip, &
154 info%lat)
155 if (iprecip == 1) then
156 irain = irain + 1
157 cycle reports
158 end if
159 end if
160
161 iv%info(ssmi_tb)%ntotal = iv%info(ssmi_tb)%ntotal + 1
162 if (outside) cycle reports
163 iv%info(ssmi_tb)%nlocal = iv%info(ssmi_tb)%nlocal + 1
164
165 case default;
166 ! Do nothing.
167 end select
168
169 end do reports
170
171 iv%info(ssmt1)%max_lev = 1
172 iv%info(ssmt2)%max_lev = 1
173 iv%info(ssmi_tb)%max_lev = 1
174 iv%info(ssmi_rv)%max_lev = 1
175
176 close(unit=iunit)
177 call da_free_unit(iunit)
178
179 if (irain > 0) then
180 write(unit=stdout, fmt='(/,5x,a,i6/)') 'Rejected rain contaminated ssmi_tb =', irain
181 write(unit=stdout, fmt='(A)') ' '
182 end if
183
184 if (trace_use) call da_trace_exit("da_scan_obs_ssmi")
185
186 end subroutine da_scan_obs_ssmi
187
188