da_read_kma1dvar.inc
References to this file elsewhere.
1 subroutine da_read_kma1dvar (inst,iv, ob, xp, infile)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: read in kma 1dvar innovation output to innovation and obs structure
5 !
6 ! METHOD: use F90 sequantial data structure to avoid read file twice
7 ! so that da_scan_bufrtovs is not necessary any more.
8 ! 1. read radiance data in sequential data structure
9 ! 2. assign sequential data structure to innovation structure
10 ! and deallocate sequential data structure
11 !---------------------------------------------------------------------------
12
13 implicit none
14
15 ! subroutine argument
16 integer , intent (in) :: inst
17 character(20) , intent (in) :: infile
18 type (xpose_type) , intent (in) :: xp
19 type (y_type) , intent (inout) :: ob
20 type (ob_type) , intent (inout) :: iv
21
22 ! local variables
23 integer :: iost, n, i, j,inunit
24
25 ! Declare local variables
26 logical outside,inside_halo
27 character(20) :: instrument
28
29 ! pixel information
30 integer :: idate, npass, scanpos
31 real :: rlat, rlon ! lat/lon in degrees for Anfovs
32 real :: satzen ! scan angles for Anfovs
33 integer :: landsea_mask
34 real :: srf_height
35 integer :: nchanl,knchan ! number of channels
36
37 real :: fastem(5)
38 real , allocatable :: & ! bright temperatures
39 otb(:),oerr(:),inov(:), emiss(:), &
40 tb(:), inv(:),err(:)
41 integer, allocatable :: kchan(:),qc(:)
42
43 type (datalink_type), pointer :: head, p, current
44
45 integer :: size, error
46 type(info_type) :: info
47 type(model_loc_type) :: loc
48
49 if (trace_use) call da_trace_entry("da_read_kma1dvar")
50
51 !**************************************************************************
52 ! Initialize variables
53
54 nchanl = iv%instid(inst)%nchan
55
56 allocate ( kchan(nchanl) )
57 allocate ( otb(nchanl) )
58 allocate ( oerr(nchanl) )
59 allocate ( inov(nchanl) )
60 allocate ( emiss(nchanl) )
61
62 allocate ( tb(nchanl) )
63 allocate ( inv(nchanl) )
64 allocate ( err(nchanl) )
65 allocate ( qc(nchanl) )
66
67 ! 0.0 Open unit for kma 1dvar file and read file header
68 !--------------------------------------------------------------
69
70 call da_get_unit(inunit)
71 open(unit=inunit,file=trim(infile),form='formatted', &
72 iostat = iost, status = 'old')
73 if (iost /= 0) then
74 call da_error(__FILE__,__LINE__,&
75 (/"Cannot open file "//trim(infile)/))
76 end if
77
78 read(unit=inunit,fmt='(A,I10,I12)') instrument,idate,npass
79
80 ! Loop to read pixel and assign information to a sequential structure
81 !-------------------------------------------------------------------------
82
83 allocate ( head )
84 nullify ( head % next )
85 p => head
86 size = 0
87
88 do j=1,npass
89 read(inunit,'(2F12.4)') rlat,rlon
90 read(inunit,'(I5,20I5)') knchan,(kchan(i),i=1,knchan)
91 ! landsea_mask: 0:land ; 1:sea (same as RTTOV)
92 read(inunit,'(I5,F12.4,I5,F12.4)') scanpos,satzen,landsea_mask,srf_height
93 read(inunit,'(20F12.4)') (oerr(i),i=1,knchan)
94 read(inunit,'(20F12.4)') (emiss(i),i=1,knchan)
95 read(inunit,'(5F12.4)') (fastem(i),i=1,5)
96 read(inunit,'(20F12.4)') (otb(i),i=1,knchan)
97 read(inunit,'(20F12.4)') (inov(i),i=1,knchan)
98
99 ! 1.0 Extract observation location and other required information
100 ! judge if data is in the domain, read next record if not
101 !-------------------------------------------------------------------
102
103 info%lat = rlat
104 info%lon = rlon
105 call da_ll_to_xy (info, loc, xp, outside, inside_halo )
106 if (outside) cycle
107
108 info%elv = srf_height
109 write(unit=info%date_char, fmt='(i10)') idate
110
111 tb(1:nchanl) = missing_r
112 inv(1:nchanl) = missing_r
113 err(1:nchanl) = missing_r
114 qc(1:nchanl) = qc_bad
115
116 do i=1,knchan
117 tb(kchan(i)) = otb(i)
118 inv(kchan(i)) = inov(i)
119 err(kchan(i)) = oerr(i)
120 qc(kchan(i)) = qc_good
121 end do
122
123 ! 2.0 assign information to sequential radiance structure
124 !--------------------------------------------------------------------
125
126 allocate (p % tb_inv (1:nchanl))
127 allocate (p % tb_error (1:nchanl))
128 allocate (p % tb_qc (1:nchanl))
129 allocate (p % tb_ob (1:nchanl))
130 p%info = info
131 p%loc = loc
132 p%landsea_mask = landsea_mask
133 p%scanpos = scanpos
134 p%satzen = satzen
135 p%satazi = missing_r
136 p%solzen = missing_r
137 p%tb_inv(1:nchanl) = inv(1:nchanl)
138 p%tb_error(1:nchanl) = err(1:nchanl)
139 p%tb_qc(1:nchanl) = qc(1:nchanl)
140 p%tb_ob(1:nchanl) = tb(1:nchanl)
141 p%sensor_index = inst
142
143 size = size + 1
144 allocate ( p%next, stat=error) ! add next data
145 if (error /= 0 ) then
146 call da_error(__FILE__,__LINE__, &
147 (/"Cannot allocate radiance structure"/))
148 end if
149
150 p => p%next
151 nullify (p%next)
152 end do
153
154 iv%total_obs = iv%total_obs + size
155 iv%total_rad_pixel = iv%total_rad_pixel + size
156 iv%total_rad_channel = iv%total_rad_channel + size*nchanl
157
158 deallocate (kchan)
159 deallocate (otb)
160 deallocate (oerr)
161 deallocate (inov)
162 deallocate (emiss)
163
164 deallocate (tb)
165 deallocate (inv)
166 deallocate (err)
167 deallocate (qc)
168
169 close(inunit)
170 call da_free_unit(inunit)
171
172 ! 3.0 allocate innovation and obs radiance structure
173 !----------------------------------------------------------------
174 iv%instid(inst)%num_rad = size
175 ob%instid(inst)%num_rad = size
176
177 write(unit=stdout,fmt='(a,i3,2x,a,3x,i10)') ' allocating space for radiance innov structure', &
178 i, iv%instid(inst)%rttovid_string, iv%instid(inst)%num_rad
179
180 ! 4.0 assign sequential structure to innovation structure
181 !-------------------------------------------------------------
182 p => head
183 allocate (iv%instid(inst)%tb_inv(1:nchanl,size))
184 allocate (iv%instid(inst)%tb_error(1:nchanl,size))
185 allocate (iv%instid(inst)%tb_qc(1:nchanl,size))
186 allocate (ob%instid(inst)%tb(1:nchanl,size))
187 do n = 1, size
188 iv%instid(inst)%info(n) = p%info
189 iv%instid(inst)%loc_i(n) = p%loc%i
190 iv%instid(inst)%loc_j(n) = p%loc%j
191 iv%instid(inst)%loc_dx(n) = p%loc%dx
192 iv%instid(inst)%loc_dy(n) = p%loc%dy
193 iv%instid(inst)%loc_dxm(n) = p%loc%dxm
194 iv%instid(inst)%loc_dym(n) = p%loc%dym
195 iv%instid(inst)%landsea_mask(n) = p%landsea_mask
196 iv%instid(inst)%scanpos(n) = p%scanpos
197 iv%instid(inst)%satzen(n) = p%satzen
198 iv%instid(inst)%satazi(n) = p%satazi
199 iv%instid(inst)%solzen(n) = p%solzen
200 iv%instid(inst)%tb_inv(1:nchanl,n) = p%tb_inv(1:nchanl)
201 iv%instid(inst)%tb_error(1:nchanl,n) = p%tb_error(1:nchanl)
202 iv%instid(inst)%tb_qc(1:nchanl,n) = p%tb_qc(1:nchanl)
203 ob%instid(inst)%tb(1:nchanl,n) = p%tb_ob(1:nchanl)
204
205 ! write(unit=stdout,*) inst, nread, iv%instid(inst)%tb_inv(1:nchanl,n)
206 current => p
207 p => p%next
208
209 ! free current data
210 deallocate (current % tb_ob)
211 deallocate (current % tb_inv)
212 deallocate (current % tb_error)
213 deallocate (current % tb_qc)
214 deallocate (current)
215 end do
216
217 ! check if sequential structure has been freed
218 !
219 ! p => head
220 ! do i = 1, size
221 ! write (unit=stdout,fmt=*) i, p%tb(1:nchanl)%inv
222 ! p => p%next
223 ! end do
224
225 if (trace_use) call da_trace_exit("da_readkma1dvar")
226
227 end subroutine da_read_kma1dvar
228
229