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