da_rad_diags.f90
References to this file elsewhere.
1 program da_rad_diags
2 !
3 ! Author: Hui-Chuan Lin
4 !
5 ! program to read multiple-time of radiance innovation files and write out in
6 ! netcdf format for ncl time-series plotting
7 !
8 ! input files: (1) namelist.rw_rad_diags
9 ! &record1
10 ! nproc = 16 (the proc numbers used when inv files were written out)
11 ! instid = 'dmsp-16-ssmis' (inst names, can be more than one instid)
12 ! file_prefix = "inv"
13 ! start_date = '2005082000'
14 ! end_date = '2005083018'
15 ! cycle_period = 6
16 ! /
17 ! (2) inv_* or oma_* from wrfvar
18
19 implicit none
20
21 #include "netcdf.inc"
22 !
23 ! namelist variables
24 !
25 ! nproc: number of processsors used when writing out inv files
26 ! instid, eg dmsp-16-ssmis
27 ! file_prefix, inv or oma
28 ! start_date, end_date, eg 2006100100, 2006102800
29 ! cycle_period (hours) between dates, eg 6 or 12
30 integer, parameter :: maxnum = 20, maxlvl = 100
31 integer :: nml_unit = 87
32 integer :: nproc, nlev, ilev, ich
33 integer :: cycle_period
34 character(len=20), dimension(maxnum) :: instid
35 character(len=3) :: file_prefix
36 character(len=10) :: start_date, end_date
37 namelist /record1/ nproc, instid, file_prefix, start_date, end_date, cycle_period
38 !
39 ! netcdf variables
40 !
41 character(len=80) :: ncname
42 integer :: ncid, dimid, varid
43 integer, dimension(3) :: ishape, istart, icount
44 !
45 logical :: isfile, prf_found, jac_found
46 integer, parameter :: datelen1 = 10
47 integer, parameter :: datelen2 = 19
48 real, parameter :: missing_r = -888888.00
49 character(len=250) :: buf, inst
50 character(len=7) :: numbuf
51 character(len=datelen1) :: valid_date
52 integer :: ninst, iinst, iproc, ipixel, ifirst
53 integer :: ios, i, n, ips, ipe, nerr, itime
54 integer :: ntime, nchan, total_npixel
55 integer, dimension(:), allocatable :: ichan, npixel, iunit, scanpos, isflg
56 integer, dimension(:), allocatable :: landsea_mask, soiltyp, vegtyp
57 real, dimension(:), allocatable :: lat, lon, elv, elev
58 real, dimension(:), allocatable :: satzen, satazi, t2m, mr2m, u10, v10, ps, ts
59 real, dimension(:), allocatable :: smois, tslb, snowh, vegfra, clwp
60 integer, dimension(:,:), allocatable :: tb_qc
61 real, dimension(:,:), allocatable :: tb_obs, tb_bak, tb_inv, tb_oma, tb_err, ems, ems_jac
62 real, dimension(:,:), allocatable :: prf_pfull, prf_phalf, prf_t, prf_q, prf_water
63 real, dimension(:,:), allocatable :: prf_ice, prf_rain, prf_snow, prf_grau, prf_hail
64 real, dimension(:,:), allocatable :: prf_water_reff, prf_ice_reff, prf_rain_reff
65 real, dimension(:,:), allocatable :: prf_snow_reff, prf_grau_reff, prf_hail_reff
66 real, dimension(:,:,:), allocatable :: prf_t_jac, prf_q_jac, prf_water_jac, prf_ice_jac
67 real, dimension(:,:,:), allocatable :: prf_rain_jac, prf_snow_jac, prf_grau_jac, prf_hail_jac
68 real, dimension(:,:,:), allocatable :: prf_water_reff_jac, prf_ice_reff_jac, prf_rain_reff_jac
69 real, dimension(:,:,:), allocatable :: prf_snow_reff_jac, prf_grau_reff_jac, prf_hail_reff_jac
70 character(len=80), dimension(:), allocatable :: inpname
71 character(len=datelen1), dimension(:), allocatable :: datestr1 ! date string
72 character(len=datelen2), dimension(:), allocatable :: datestr2 ! date string
73 !
74 ! initialize some variables
75
76 instid = ''
77 ninst = 0
78 file_prefix = 'inv'
79 ilev = 0
80 nlev = 0
81 prf_found = .false.
82 jac_found = .false.
83 !
84 ! read namelist
85 !
86 open(unit=nml_unit, file='namelist.da_rad_diags', status='old', form='formatted')
87 read(unit=nml_unit, nml=record1, iostat=ios)
88 write(0,nml=record1)
89 if ( ios /= 0 ) then
90 write(0,*) 'Error reading namelist record1'
91 stop
92 end if
93 !
94 ! find out how many instruments to process
95 !
96 do i = 1, maxnum
97 if ( len_trim(instid(i)) /= 0 ) then
98 ninst = ninst + 1
99 end if
100 end do
101 !
102 ! find out how many dates to process
103 !
104 valid_date = start_date
105 do while ( valid_date <= end_date )
106 ntime = ntime + 1
107 call advance_cymdh(valid_date, cycle_period, valid_date)
108 end do
109 write(0,*) 'ntime = ', ntime
110
111 allocate ( datestr1(ntime) )
112 valid_date = start_date
113 datestr1(1) = start_date
114 do i = 2, ntime
115 call advance_cymdh(datestr1(i-1), cycle_period, datestr1(i))
116 end do
117
118 ntime_loop: do itime = 1, ntime
119
120 write(0,*) '=================='
121 write(0,*) trim(datestr1(itime))
122
123 ninst_loop: do iinst = 1, ninst
124
125 write(0,*) '------------------'
126 write(0,*) trim(instid(iinst))
127
128 nerr = 0
129 total_npixel = 0
130 ips = 0
131 ipe = 0
132
133 allocate ( npixel(0:nproc-1) )
134 allocate ( iunit(0:nproc-1) )
135 allocate ( inpname(0:nproc-1) )
136
137 nproc_loop_1: do iproc = 0, nproc - 1 ! loop for first getting number of pixels from each proc
138
139 if ( nproc > 1 ) then
140 write(unit=inpname(iproc), fmt='(a,i2.2)') &
141 trim(datestr1(itime))//'/'//file_prefix//'_'//trim(instid(iinst))//'.', iproc
142 else
143 write(unit=inpname(iproc), fmt='(a,i2.2)') &
144 trim(datestr1(itime))//'/'//file_prefix//'_'//trim(instid(iinst))
145 end if
146 iunit(iproc) = 101 + iproc
147 inquire(file=trim(inpname(iproc)), exist=isfile)
148 if ( .not. isfile ) Then
149 write(0,*) 'Error opening innovation radiance file ', trim(inpname(iproc))
150 nerr = nerr + 1
151 if ( nerr == nproc ) then
152 write(0,*) 'found no vaild files for ', trim(instid(iinst))
153 deallocate ( npixel )
154 deallocate ( iunit )
155 deallocate ( inpname )
156 cycle ninst_loop
157 end if
158 cycle nproc_loop_1
159 end if
160
161 open(unit=iunit(iproc),file=trim(inpname(iproc)),form='formatted',iostat=ios)
162 read(unit=iunit(iproc),fmt='(a)') buf ! first read in one line
163
164 inst = buf(1:(index(buf,'number-of-pixels')-2)) ! retrieve inst name
165 !
166 ! retrieve number of pixels
167 !
168 numbuf = buf((index(buf,'channel-number')-8):(index(buf,'channel-number')-2))
169 read(numbuf,'(i7)') npixel(iproc)
170
171 total_npixel = total_npixel + npixel(iproc)
172
173 end do nproc_loop_1
174
175 write(0,*) 'total_npixel = ', total_npixel
176
177 ifirst = 1
178 nproc_loop_2: do iproc = 0, nproc - 1
179
180 inquire(file=trim(inpname(iproc)), exist=isfile)
181 if ( .not. isfile ) cycle nproc_loop_2
182 rewind(iunit(iproc))
183 read(unit=iunit(iproc),fmt='(a)') buf
184 !
185 ! retrieve number of channels
186 !
187 numbuf = buf((index(buf,'index-of-channels')-6):(index(buf,'index-of-channels')-2))
188 read(numbuf,'(i5)') nchan
189
190 if ( .not. allocated(ichan) ) allocate ( ichan(1:nchan) )
191
192 read(unit=iunit(iproc),fmt='(10i5)',iostat=ios) ichan
193 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! pixel-info line
194 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! xb-surf-info line
195
196 if ( ifirst == 1 ) then
197
198 allocate ( datestr2(1:total_npixel) )
199 allocate ( scanpos(1:total_npixel) )
200 allocate ( landsea_mask(1:total_npixel) )
201 allocate ( elv(1:total_npixel) )
202 allocate ( lat(1:total_npixel) )
203 allocate ( lon(1:total_npixel) )
204 allocate ( satzen(1:total_npixel) )
205 allocate ( satazi(1:total_npixel) )
206 allocate ( t2m(1:total_npixel) )
207 allocate ( mr2m(1:total_npixel) )
208 allocate ( u10(1:total_npixel) )
209 allocate ( v10(1:total_npixel) )
210 allocate ( ps(1:total_npixel) )
211 allocate ( ts(1:total_npixel) )
212 allocate ( smois(1:total_npixel) )
213 allocate ( tslb(1:total_npixel) )
214 allocate ( snowh(1:total_npixel) )
215 allocate ( isflg(1:total_npixel) )
216 allocate ( soiltyp(1:total_npixel) )
217 allocate ( vegtyp(1:total_npixel) )
218 allocate ( vegfra(1:total_npixel) )
219 allocate ( elev(1:total_npixel) )
220 allocate ( clwp(1:total_npixel) )
221 allocate ( tb_obs(1:nchan,1:total_npixel) )
222 allocate ( tb_bak(1:nchan,1:total_npixel) )
223 allocate ( tb_inv(1:nchan,1:total_npixel) )
224 allocate ( tb_oma(1:nchan,1:total_npixel) )
225 allocate ( tb_err(1:nchan,1:total_npixel) )
226 allocate ( tb_qc(1:nchan,1:total_npixel) )
227 allocate ( ems(1:nchan,1:total_npixel) )
228 allocate ( ems_jac(1:nchan,1:total_npixel) )
229 allocate ( prf_pfull(1:maxlvl,1:total_npixel) )
230 allocate ( prf_phalf(1:maxlvl,1:total_npixel) )
231 allocate ( prf_t(1:maxlvl,1:total_npixel) )
232 allocate ( prf_q(1:maxlvl,1:total_npixel) )
233 allocate ( prf_water(1:maxlvl,1:total_npixel) )
234 allocate ( prf_ice(1:maxlvl,1:total_npixel) )
235 allocate ( prf_rain(1:maxlvl,1:total_npixel) )
236 allocate ( prf_snow(1:maxlvl,1:total_npixel) )
237 allocate ( prf_grau(1:maxlvl,1:total_npixel) )
238 allocate ( prf_hail(1:maxlvl,1:total_npixel) )
239 allocate ( prf_water_reff(1:maxlvl,1:total_npixel) )
240 allocate ( prf_ice_reff(1:maxlvl,1:total_npixel) )
241 allocate ( prf_rain_reff(1:maxlvl,1:total_npixel) )
242 allocate ( prf_snow_reff(1:maxlvl,1:total_npixel) )
243 allocate ( prf_grau_reff(1:maxlvl,1:total_npixel) )
244 allocate ( prf_hail_reff(1:maxlvl,1:total_npixel) )
245 allocate ( prf_t_jac(1:maxlvl,1:nchan,1:total_npixel) )
246 allocate ( prf_q_jac(1:maxlvl,1:nchan,1:total_npixel) )
247 allocate ( prf_water_jac(1:maxlvl,1:nchan,1:total_npixel) )
248 allocate ( prf_ice_jac(1:maxlvl,1:nchan,1:total_npixel) )
249 allocate ( prf_rain_jac(1:maxlvl,1:nchan,1:total_npixel) )
250 allocate ( prf_snow_jac(1:maxlvl,1:nchan,1:total_npixel) )
251 allocate ( prf_grau_jac(1:maxlvl,1:nchan,1:total_npixel) )
252 allocate ( prf_hail_jac(1:maxlvl,1:nchan,1:total_npixel) )
253 allocate ( prf_water_reff_jac(1:maxlvl,1:nchan,1:total_npixel) )
254 allocate ( prf_ice_reff_jac(1:maxlvl,1:nchan,1:total_npixel) )
255 allocate ( prf_rain_reff_jac(1:maxlvl,1:nchan,1:total_npixel) )
256 allocate ( prf_snow_reff_jac(1:maxlvl,1:nchan,1:total_npixel) )
257 allocate ( prf_grau_reff_jac(1:maxlvl,1:nchan,1:total_npixel) )
258 allocate ( prf_hail_reff_jac(1:maxlvl,1:nchan,1:total_npixel) )
259 ! initialize
260 tb_obs = missing_r
261 tb_bak = missing_r
262 tb_inv = missing_r
263 tb_oma = missing_r
264 tb_err = missing_r
265 ncname = 'diags_'//trim(instid(iinst))//"_"//datestr1(itime)//'.nc'
266 ios = NF_CREATE(ncname, NF_CLOBBER, ncid) ! NF_CLOBBER specifies the default behavior of
267 ! overwritting any existing dataset with the
268 ! same file name
269 if ( ios /= 0 ) then
270 write(0,*) 'Error creating netcdf file: ', ncname
271 stop
272 end if
273
274 ifirst = 0
275
276 end if ! end of ifirst if
277 !
278 ! decide the index of pixels to store data in
279 !
280 ips = ipe + 1
281 ipe = ipe + npixel(iproc)
282 write(0,*) 'Processing pixels ', ips, ' to ', ipe
283
284 npixel_loop: do ipixel = ips, ipe
285
286 read(unit=iunit(iproc),fmt='(7x,i7,2x,a19,i3,i3,f6.0,4f8.2)',iostat=ios) &
287 n, datestr2(ipixel), scanpos(ipixel), landsea_mask(ipixel), elv(ipixel), &
288 lat(ipixel), lon(ipixel), satzen(ipixel), satazi(ipixel)
289 if ( scanpos(ipixel) > 360 .or. scanpos(ipixel) == 0 ) then
290 backspace(iunit(iproc))
291 read(unit=iunit(iproc),fmt='(7x,i7,2x,a19,i6,i3,f6.0,4f8.2)',iostat=ios) &
292 n, datestr2(ipixel), scanpos(ipixel), landsea_mask(ipixel), elv(ipixel), &
293 lat(ipixel), lon(ipixel), satzen(ipixel), satazi(ipixel)
294 end if
295 read(unit=iunit(iproc),fmt='(14x,9f10.2,3i3,3f10.2)',iostat=ios) &
296 t2m(ipixel), mr2m(ipixel), u10(ipixel), v10(ipixel), ps(ipixel), ts(ipixel), &
297 smois(ipixel), tslb(ipixel), snowh(ipixel), isflg(ipixel), soiltyp(ipixel), &
298 vegtyp(ipixel), vegfra(ipixel), elev(ipixel), clwp(ipixel)
299 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! OBS
300 read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_obs(:,ipixel)
301 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! BAK
302 read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_bak(:,ipixel)
303 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! IVBC
304 read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_inv(:,ipixel)
305 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! OMA or EMS
306 if ( buf(1:3) == "OMA" ) then
307 read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_oma(:,ipixel)
308 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! OMA or EMS
309 end if
310 read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) ems(:,ipixel)
311 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! EMS_JACOBIAN or ERR
312 if ( buf(1:12) == "EMS_JACOBIAN" ) then
313 jac_found = .true.
314 read(unit=iunit(iproc),fmt='(10f10.3)',iostat=ios) ems_jac(:,ipixel)
315 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! OMA or EMS
316 end if
317 read(unit=iunit(iproc),fmt='(10f11.2)',iostat=ios) tb_err(:,ipixel)
318 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf ! QC
319 read(unit=iunit(iproc),fmt='(10i11)',iostat=ios ) tb_qc(:,ipixel)
320 read(unit=iunit(iproc),fmt='(a)',iostat=ios) buf
321 if ( buf(1:4) == "INFO" ) then
322 backspace(iunit(iproc))
323 cycle npixel_loop
324 else
325 if ( buf(1:6) == " level" ) then ! profiles are available
326 prf_found = .true.
327 ilev = 0
328 do while ( ios == 0 )
329 ilev = ilev + 1
330 read(unit=iunit(iproc),fmt='(3x,2f10.2,f8.2,13f8.3)',iostat=ios) &
331 prf_pfull(ilev,ipixel), prf_phalf(ilev,ipixel), prf_t(ilev,ipixel), &
332 prf_q(ilev,ipixel), prf_water(ilev,ipixel), prf_ice(ilev,ipixel), &
333 prf_rain(ilev,ipixel), prf_snow(ilev,ipixel), prf_grau(ilev,ipixel), &
334 prf_hail(ilev,ipixel),prf_water_reff(ilev,ipixel), &
335 prf_ice_reff(ilev,ipixel), prf_rain_reff(ilev,ipixel), &
336 prf_snow_reff(ilev,ipixel), prf_grau_reff(ilev,ipixel), &
337 prf_hail_reff(ilev,ipixel)
338 end do
339 nlev = ilev - 1
340 ! backspace(iunit(iproc))
341 end if
342 if ( jac_found ) then
343 ios = 0
344 do while ( ios == 0 )
345 ! ilev = ilev + 1
346 read(unit=iunit(iproc),fmt='(i5,i3,10x,14f8.3)',iostat=ios) &
347 ich, ilev, prf_t_jac(ilev,ich,ipixel), prf_q_jac(ilev,ich,ipixel), &
348 prf_water_jac(ilev,ich,ipixel), prf_ice_jac(ilev,ich,ipixel), &
349 prf_rain_jac(ilev,ich,ipixel), prf_snow_jac(ilev,ich,ipixel), &
350 prf_grau_jac(ilev,ich,ipixel), prf_hail_jac(ilev,ich,ipixel), &
351 prf_water_reff_jac(ilev,ich,ipixel), prf_ice_reff_jac(ilev,ich,ipixel), &
352 prf_rain_reff_jac(ilev,ich,ipixel), prf_snow_reff_jac(ilev,ich,ipixel), &
353 prf_grau_reff_jac(ilev,ich,ipixel), prf_hail_reff_jac(ilev,ich,ipixel)
354 nlev = max(nlev,ilev)
355 end do
356 backspace(iunit(iproc))
357 else
358 backspace(iunit(iproc))
359 end if
360 end if
361
362 end do npixel_loop
363
364 end do nproc_loop_2
365
366 write(0,*) 'Writing out data in netCDF format...'
367 !
368 ! define dimensions
369 !
370 ios = NF_DEF_DIM(ncid, 'nchan', nchan, dimid)
371 ios = NF_DEF_DIM(ncid, 'npixel', total_npixel, dimid)
372 ios = NF_DEF_DIM(ncid, 'DateStrLen', datelen2, dimid)
373 if ( prf_found .or. jac_found ) then
374 ios = NF_DEF_DIM(ncid, 'nlev', nlev, dimid)
375 end if
376 !
377 ! define variables
378 !
379 ! define 2-D array for date string
380 !
381 ios = NF_INQ_DIMID(ncid, 'DateStrLen', dimid)
382 ishape(1) = dimid
383 ios = NF_INQ_DIMID(ncid, 'npixel', dimid)
384 ishape(2) = dimid
385 ios = NF_DEF_VAR(ncid, 'date', NF_CHAR, 2, ishape(1:2), varid)
386 !
387 ! define 2-D array with dimensions nchan * total_npixel
388 !
389 ios = NF_INQ_DIMID(ncid, 'nchan', dimid)
390 ishape(1) = dimid
391 ios = NF_INQ_DIMID(ncid, 'npixel', dimid)
392 ishape(2) = dimid
393 ios = NF_DEF_VAR(ncid, 'tb_obs', NF_FLOAT, 2, ishape(1:2), varid)
394 ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
395 ios = NF_DEF_VAR(ncid, 'tb_bak', NF_FLOAT, 2, ishape(1:2), varid)
396 ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
397 ios = NF_DEF_VAR(ncid, 'tb_inv', NF_FLOAT, 2, ishape(1:2), varid)
398 ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
399 ios = NF_DEF_VAR(ncid, 'tb_oma', NF_FLOAT, 2, ishape(1:2), varid)
400 ios = NF_PUT_ATT_REAL(ncid, varid, 'missing_value', NF_FLOAT, 1, missing_r)
401 ios = NF_DEF_VAR(ncid, 'ems', NF_FLOAT, 2, ishape(1:2), varid)
402 if ( jac_found ) then
403 ios = NF_DEF_VAR(ncid, 'ems_jac',NF_FLOAT, 2, ishape(1:2), varid)
404 end if
405 ios = NF_DEF_VAR(ncid, 'tb_err', NF_FLOAT, 2, ishape(1:2), varid)
406 ios = NF_DEF_VAR(ncid, 'tb_qc', NF_INT, 2, ishape(1:2), varid)
407 !
408 ! define 2-D array with dimensions nlev * total_npixel
409 !
410 if ( prf_found ) then
411 ios = NF_INQ_DIMID(ncid, 'nlev', dimid)
412 ishape(1) = dimid
413 ios = NF_INQ_DIMID(ncid, 'npixel', dimid)
414 ishape(2) = dimid
415 ios = NF_DEF_VAR(ncid, 'prf_pfull', NF_FLOAT, 2, ishape(1:2), varid)
416 ios = NF_DEF_VAR(ncid, 'prf_phalf', NF_FLOAT, 2, ishape(1:2), varid)
417 ios = NF_DEF_VAR(ncid, 'prf_t', NF_FLOAT, 2, ishape(1:2), varid)
418 ios = NF_DEF_VAR(ncid, 'prf_q', NF_FLOAT, 2, ishape(1:2), varid)
419 ios = NF_DEF_VAR(ncid, 'prf_water', NF_FLOAT, 2, ishape(1:2), varid)
420 ios = NF_DEF_VAR(ncid, 'prf_ice', NF_FLOAT, 2, ishape(1:2), varid)
421 ios = NF_DEF_VAR(ncid, 'prf_rain', NF_FLOAT, 2, ishape(1:2), varid)
422 ios = NF_DEF_VAR(ncid, 'prf_snow', NF_FLOAT, 2, ishape(1:2), varid)
423 ios = NF_DEF_VAR(ncid, 'prf_grau', NF_FLOAT, 2, ishape(1:2), varid)
424 ios = NF_DEF_VAR(ncid, 'prf_hail', NF_FLOAT, 2, ishape(1:2), varid)
425 ios = NF_DEF_VAR(ncid, 'prf_water_reff', NF_FLOAT, 2, ishape(1:2), varid)
426 ios = NF_DEF_VAR(ncid, 'prf_ice_reff', NF_FLOAT, 2, ishape(1:2), varid)
427 ios = NF_DEF_VAR(ncid, 'prf_rain_reff', NF_FLOAT, 2, ishape(1:2), varid)
428 ios = NF_DEF_VAR(ncid, 'prf_snow_reff', NF_FLOAT, 2, ishape(1:2), varid)
429 ios = NF_DEF_VAR(ncid, 'prf_grau_reff', NF_FLOAT, 2, ishape(1:2), varid)
430 ios = NF_DEF_VAR(ncid, 'prf_hail_reff', NF_FLOAT, 2, ishape(1:2), varid)
431 end if
432 !
433 ! define 3-D array with dimensions nlev * nchan * total_npixel
434 !
435 if ( jac_found ) then
436 ios = NF_INQ_DIMID(ncid, 'nlev', dimid)
437 ishape(1) = dimid
438 ios = NF_INQ_DIMID(ncid, 'nchan', dimid)
439 ishape(2) = dimid
440 ios = NF_INQ_DIMID(ncid, 'npixel', dimid)
441 ishape(3) = dimid
442 ios = NF_DEF_VAR(ncid, 'prf_t_jac', NF_FLOAT, 3, ishape, varid)
443 ios = NF_DEF_VAR(ncid, 'prf_q_jac', NF_FLOAT, 3, ishape, varid)
444 ios = NF_DEF_VAR(ncid, 'prf_water_jac', NF_FLOAT, 3, ishape, varid)
445 ios = NF_DEF_VAR(ncid, 'prf_ice_jac', NF_FLOAT, 3, ishape, varid)
446 ios = NF_DEF_VAR(ncid, 'prf_rain_jac', NF_FLOAT, 3, ishape, varid)
447 ios = NF_DEF_VAR(ncid, 'prf_snow_jac', NF_FLOAT, 3, ishape, varid)
448 ios = NF_DEF_VAR(ncid, 'prf_grau_jac', NF_FLOAT, 3, ishape, varid)
449 ios = NF_DEF_VAR(ncid, 'prf_hail_jac', NF_FLOAT, 3, ishape, varid)
450 ios = NF_DEF_VAR(ncid, 'prf_water_reff_jac', NF_FLOAT, 3, ishape, varid)
451 ios = NF_DEF_VAR(ncid, 'prf_ice_reff_jac', NF_FLOAT, 3, ishape, varid)
452 ios = NF_DEF_VAR(ncid, 'prf_rain_reff_jac', NF_FLOAT, 3, ishape, varid)
453 ios = NF_DEF_VAR(ncid, 'prf_snow_reff_jac', NF_FLOAT, 3, ishape, varid)
454 ios = NF_DEF_VAR(ncid, 'prf_grau_reff_jac', NF_FLOAT, 3, ishape, varid)
455 ios = NF_DEF_VAR(ncid, 'prf_hail_reff_jac', NF_FLOAT, 3, ishape, varid)
456 end if
457 !
458 ! define 1-D array with dimension nchan
459 !
460 ios = NF_INQ_DIMID(ncid, 'nchan', dimid)
461 ishape(1) = dimid
462 ios = NF_DEF_VAR(ncid, 'ichan', NF_INT, 1, ishape(1), varid) ! channel index
463 !
464 ! define 1-D array with dimension total_npixel
465 !
466 ios = NF_INQ_DIMID(ncid, 'npixel', dimid)
467 ishape(1) = dimid
468 ios = NF_DEF_VAR(ncid, 'scanpos', NF_INT, 1, ishape(1), varid)
469 ios = NF_DEF_VAR(ncid, 'landsea_mask', NF_INT, 1, ishape(1), varid)
470 ios = NF_DEF_VAR(ncid, 'elv', NF_FLOAT, 1, ishape(1), varid)
471 ios = NF_DEF_VAR(ncid, 'lat', NF_FLOAT, 1, ishape(1), varid)
472 ios = NF_DEF_VAR(ncid, 'lon', NF_FLOAT, 1, ishape(1), varid)
473 ios = NF_DEF_VAR(ncid, 'satzen', NF_FLOAT, 1, ishape(1), varid)
474 ios = NF_DEF_VAR(ncid, 'satazi', NF_FLOAT, 1, ishape(1), varid)
475 ios = NF_DEF_VAR(ncid, 't2m', NF_FLOAT, 1, ishape(1), varid)
476 ios = NF_DEF_VAR(ncid, 'mr2m', NF_FLOAT, 1, ishape(1), varid)
477 ios = NF_DEF_VAR(ncid, 'u10', NF_FLOAT, 1, ishape(1), varid)
478 ios = NF_DEF_VAR(ncid, 'v10', NF_FLOAT, 1, ishape(1), varid)
479 ios = NF_DEF_VAR(ncid, 'ps', NF_FLOAT, 1, ishape(1), varid)
480 ios = NF_DEF_VAR(ncid, 'ts', NF_FLOAT, 1, ishape(1), varid)
481 ios = NF_DEF_VAR(ncid, 'smois', NF_FLOAT, 1, ishape(1), varid)
482 ios = NF_DEF_VAR(ncid, 'snowh', NF_FLOAT, 1, ishape(1), varid)
483 ios = NF_DEF_VAR(ncid, 'isflg', NF_INT, 1, ishape(1), varid)
484 ios = NF_DEF_VAR(ncid, 'soiltyp', NF_INT, 1, ishape(1), varid)
485 ios = NF_DEF_VAR(ncid, 'vegtyp', NF_INT, 1, ishape(1), varid)
486 ios = NF_DEF_VAR(ncid, 'vegfra', NF_FLOAT, 1, ishape(1), varid)
487 ios = NF_DEF_VAR(ncid, 'elev', NF_FLOAT, 1, ishape(1), varid)
488 ios = NF_DEF_VAR(ncid, 'clwp', NF_FLOAT, 1, ishape(1), varid)
489
490 ios = NF_ENDDEF(ncid)
491 !
492 ! output date string
493 istart(1) = 1
494 istart(2) = 1
495 icount(1) = datelen2
496 icount(2) = total_npixel
497 ios = NF_INQ_VARID (ncid, 'date', varid)
498 ios = NF_PUT_VARA_TEXT(ncid, varid, istart, icount, datestr2)
499 !
500 ! output 2-D array with dimensions nchan * total_npixel
501 !
502 istart(1) = 1
503 istart(2) = 1
504 icount(1) = nchan
505 icount(2) = total_npixel
506 ios = NF_INQ_VARID (ncid, 'tb_obs', varid)
507 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_obs)
508 ios = NF_INQ_VARID (ncid, 'tb_bak', varid)
509 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_bak)
510 ios = NF_INQ_VARID (ncid, 'tb_inv', varid)
511 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_inv)
512 ios = NF_INQ_VARID (ncid, 'tb_oma', varid)
513 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_oma)
514 ios = NF_INQ_VARID (ncid, 'ems', varid)
515 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), ems)
516 if ( jac_found ) then
517 ios = NF_INQ_VARID (ncid, 'ems_jac', varid)
518 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), ems_jac)
519 end if
520 ios = NF_INQ_VARID (ncid, 'tb_err', varid)
521 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), tb_err)
522 ios = NF_INQ_VARID (ncid, 'tb_qc', varid)
523 ios = NF_PUT_VARA_INT(ncid, varid, istart(1:2), icount(1:2), tb_qc)
524 !
525 ! output 2-D array with dimensions nlev * total_npixel
526 !
527 if ( prf_found ) then
528 istart(1) = 1
529 istart(2) = 1
530 icount(1) = nlev
531 icount(2) = total_npixel
532 ios = NF_INQ_VARID (ncid, 'prf_pfull', varid)
533 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_pfull(1:nlev,:))
534 ios = NF_INQ_VARID (ncid, 'prf_phalf', varid)
535 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_phalf(1:nlev,:))
536 ios = NF_INQ_VARID (ncid, 'prf_t', varid)
537 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_t(1:nlev,:))
538 ios = NF_INQ_VARID (ncid, 'prf_q', varid)
539 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_q(1:nlev,:))
540 ios = NF_INQ_VARID (ncid, 'prf_water', varid)
541 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_water(1:nlev,:))
542 ios = NF_INQ_VARID (ncid, 'prf_ice', varid)
543 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_ice(1:nlev,:))
544 ios = NF_INQ_VARID (ncid, 'prf_rain', varid)
545 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_rain(1:nlev,:))
546 ios = NF_INQ_VARID (ncid, 'prf_snow', varid)
547 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_snow(1:nlev,:))
548 ios = NF_INQ_VARID (ncid, 'prf_grau', varid)
549 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_grau(1:nlev,:))
550 ios = NF_INQ_VARID (ncid, 'prf_hail', varid)
551 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_hail(1:nlev,:))
552 ios = NF_INQ_VARID (ncid, 'prf_water_reff', varid)
553 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_water_reff(1:nlev,:))
554 ios = NF_INQ_VARID (ncid, 'prf_ice_reff', varid)
555 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_ice_reff(1:nlev,:))
556 ios = NF_INQ_VARID (ncid, 'prf_rain_reff', varid)
557 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_rain_reff(1:nlev,:))
558 ios = NF_INQ_VARID (ncid, 'prf_snow_reff', varid)
559 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_snow_reff(1:nlev,:))
560 ios = NF_INQ_VARID (ncid, 'prf_grau_reff', varid)
561 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_grau_reff(1:nlev,:))
562 ios = NF_INQ_VARID (ncid, 'prf_hail_reff', varid)
563 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:2), icount(1:2), prf_hail_reff(1:nlev,:))
564 end if
565 !
566 ! output 3-D array with dimensions nlev * nchan * total_npixel
567 !
568 if ( jac_found ) then
569 istart(1) = 1
570 istart(2) = 1
571 istart(3) = 1
572 icount(1) = nlev
573 icount(2) = nchan
574 icount(3) = total_npixel
575 ios = NF_INQ_VARID (ncid, 'prf_t_jac', varid)
576 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_t_jac(1:nlev,:,:))
577 ios = NF_INQ_VARID (ncid, 'prf_q_jac', varid)
578 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_q_jac(1:nlev,:,:))
579 ios = NF_INQ_VARID (ncid, 'prf_water_jac', varid)
580 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_water_jac(1:nlev,:,:))
581 ios = NF_INQ_VARID (ncid, 'prf_ice_jac', varid)
582 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_ice_jac(1:nlev,:,:))
583 ios = NF_INQ_VARID (ncid, 'prf_rain_jac', varid)
584 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_rain_jac(1:nlev,:,:))
585 ios = NF_INQ_VARID (ncid, 'prf_snow_jac', varid)
586 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_snow_jac(1:nlev,:,:))
587 ios = NF_INQ_VARID (ncid, 'prf_grau_jac', varid)
588 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_grau_jac(1:nlev,:,:))
589 ios = NF_INQ_VARID (ncid, 'prf_hail_jac', varid)
590 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_hail_jac(1:nlev,:,:))
591 ios = NF_INQ_VARID (ncid, 'prf_water_reff_jac', varid)
592 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_water_reff_jac(1:nlev,:,:))
593 ios = NF_INQ_VARID (ncid, 'prf_ice_reff_jac', varid)
594 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_ice_reff_jac(1:nlev,:,:))
595 ios = NF_INQ_VARID (ncid, 'prf_rain_reff_jac', varid)
596 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_rain_reff_jac(1:nlev,:,:))
597 ios = NF_INQ_VARID (ncid, 'prf_snow_reff_jac', varid)
598 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_snow_reff_jac(1:nlev,:,:))
599 ios = NF_INQ_VARID (ncid, 'prf_grau_reff_jac', varid)
600 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_grau_reff_jac(1:nlev,:,:))
601 ios = NF_INQ_VARID (ncid, 'prf_hail_reff_jac', varid)
602 ios = NF_PUT_VARA_REAL(ncid, varid, istart(1:3), icount(1:3), prf_hail_reff_jac(1:nlev,:,:))
603 end if
604 !
605 ! output 1-D array with dimension nchan
606 !
607 istart(1) = 1
608 icount(1) = nchan
609 ios = NF_INQ_VARID (ncid, 'ichan', varid)
610 ios = NF_PUT_VARA_INT(ncid, varid, istart(1), icount(1), ichan)
611 !
612 ! output 1-D array with dimension total_npixel
613 !
614 istart(2) = 1
615 icount(2) = total_npixel
616 ios = NF_INQ_VARID (ncid, 'scanpos', varid)
617 ios = NF_PUT_VARA_INT(ncid, varid, istart(2), icount(2), scanpos)
618 ios = NF_INQ_VARID (ncid, 'landsea_mask', varid)
619 ios = NF_PUT_VARA_INT(ncid, varid, istart(2), icount(2), landsea_mask)
620 ios = NF_INQ_VARID (ncid, 'elv', varid)
621 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), elv)
622 ios = NF_INQ_VARID (ncid, 'lat', varid)
623 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), lat)
624 ios = NF_INQ_VARID (ncid, 'lon', varid)
625 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), lon)
626 ios = NF_INQ_VARID (ncid, 'satzen', varid)
627 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), satzen)
628 ios = NF_INQ_VARID (ncid, 'satazi', varid)
629 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), satazi)
630 ios = NF_INQ_VARID (ncid, 't2m', varid)
631 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), t2m)
632 ios = NF_INQ_VARID (ncid, 'mr2m', varid)
633 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), mr2m)
634 ios = NF_INQ_VARID (ncid, 'u10', varid)
635 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), u10)
636 ios = NF_INQ_VARID (ncid, 'v10', varid)
637 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), v10)
638 ios = NF_INQ_VARID (ncid, 'ps', varid)
639 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), ps)
640 ios = NF_INQ_VARID (ncid, 'ts', varid)
641 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), ts)
642 ios = NF_INQ_VARID (ncid, 'smois', varid)
643 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), smois)
644 ios = NF_INQ_VARID (ncid, 'tslb', varid)
645 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), tslb)
646 ios = NF_INQ_VARID (ncid, 'snowh', varid)
647 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), snowh)
648 ios = NF_INQ_VARID (ncid, 'isflg', varid)
649 ios = NF_PUT_VARA_INT(ncid, varid, istart(2), icount(2), isflg)
650 ios = NF_INQ_VARID (ncid, 'soiltyp', varid)
651 ios = NF_PUT_VARA_INT(ncid, varid, istart(2), icount(2), soiltyp)
652 ios = NF_INQ_VARID (ncid, 'vegtyp', varid)
653 ios = NF_PUT_VARA_INT(ncid, varid, istart(2), icount(2), vegtyp)
654 ios = NF_INQ_VARID (ncid, 'vegfra', varid)
655 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), vegfra)
656 ios = NF_INQ_VARID (ncid, 'elev', varid)
657 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), elev)
658 ios = NF_INQ_VARID (ncid, 'clwp', varid)
659 ios = NF_PUT_VARA_REAL(ncid, varid, istart(2), icount(2), clwp)
660 !
661 ios = NF_CLOSE(ncid)
662
663 deallocate ( npixel )
664 deallocate ( iunit )
665 deallocate ( inpname )
666
667 deallocate ( ichan )
668 deallocate ( datestr2 )
669 deallocate ( scanpos )
670 deallocate ( landsea_mask )
671 deallocate ( elv )
672 deallocate ( lat )
673 deallocate ( lon )
674 deallocate ( satzen )
675 deallocate ( satazi )
676 deallocate ( t2m )
677 deallocate ( mr2m )
678 deallocate ( u10 )
679 deallocate ( v10 )
680 deallocate ( ps )
681 deallocate ( ts )
682 deallocate ( smois )
683 deallocate ( tslb )
684 deallocate ( snowh )
685 deallocate ( isflg )
686 deallocate ( soiltyp )
687 deallocate ( vegtyp )
688 deallocate ( vegfra )
689 deallocate ( elev )
690 deallocate ( clwp )
691 deallocate ( tb_obs )
692 deallocate ( tb_bak )
693 deallocate ( tb_inv )
694 deallocate ( tb_oma )
695 deallocate ( ems )
696 deallocate ( ems_jac )
697 deallocate ( tb_err )
698 deallocate ( tb_qc )
699 deallocate ( prf_pfull )
700 deallocate ( prf_phalf )
701 deallocate ( prf_t )
702 deallocate ( prf_q )
703 deallocate ( prf_water )
704 deallocate ( prf_ice )
705 deallocate ( prf_rain )
706 deallocate ( prf_snow )
707 deallocate ( prf_grau )
708 deallocate ( prf_hail )
709 deallocate ( prf_water_reff )
710 deallocate ( prf_ice_reff )
711 deallocate ( prf_rain_reff )
712 deallocate ( prf_snow_reff )
713 deallocate ( prf_grau_reff )
714 deallocate ( prf_hail_reff )
715 deallocate ( prf_t_jac )
716 deallocate ( prf_q_jac )
717 deallocate ( prf_water_jac )
718 deallocate ( prf_ice_jac )
719 deallocate ( prf_rain_jac )
720 deallocate ( prf_snow_jac )
721 deallocate ( prf_grau_jac )
722 deallocate ( prf_hail_jac )
723 deallocate ( prf_water_reff_jac )
724 deallocate ( prf_ice_reff_jac )
725 deallocate ( prf_rain_reff_jac )
726 deallocate ( prf_snow_reff_jac )
727 deallocate ( prf_grau_reff_jac )
728 deallocate ( prf_hail_reff_jac )
729 do i = 0, nproc-1
730 close(iunit(i))
731 end do
732
733 end do ninst_loop
734
735 end do ntime_loop
736
737 deallocate ( datestr1 )
738
739 end program da_rad_diags
740
741 subroutine advance_cymdh(currentdate,dh,newdate)
742
743 implicit none
744
745 character(len=10), intent(in) :: currentdate
746 integer, intent(in) :: dh
747 character(len=10), intent(out) :: newdate
748
749 integer :: ccyy, mm, dd, hh
750
751 read(currentdate(1:10), fmt='(i4, 3i2)') ccyy, mm, dd, hh
752 hh = hh + dh
753 do while (hh < 0)
754 hh = hh + 24
755 call change_date ( ccyy, mm, dd, -1 )
756 end do
757 do while (hh > 23)
758 hh = hh - 24
759 call change_date ( ccyy, mm, dd, 1 )
760 end do
761 write(newdate,'(i4.4,3(i2.2))') ccyy, mm, dd, hh
762
763 end subroutine advance_cymdh
764
765 subroutine change_date( ccyy, mm, dd, delta )
766
767 implicit none
768
769 integer, intent(inout) :: ccyy, mm, dd
770 integer, intent(in) :: delta
771 integer, dimension(12) :: mmday
772
773 mmday = (/31,28,31,30,31,30,31,31,30,31,30,31/)
774 mmday(2) = 28
775 if (mod(ccyy,4) == 0) then
776 mmday(2) = 29
777 if ( mod(ccyy,100) == 0) then
778 mmday(2) = 28
779 endif
780 if(mod(ccyy,400) == 0) then
781 mmday(2) = 29
782 end if
783 endif
784 dd = dd + delta
785 if(dd == 0) then
786 mm = mm - 1
787 if(mm == 0) then
788 mm = 12
789 ccyy = ccyy - 1
790 endif
791 dd = mmday(mm)
792 elseif ( dd .gt. mmday(mm) ) then
793 dd = 1
794 mm = mm + 1
795 if(mm > 12 ) then
796 mm = 1
797 ccyy = ccyy + 1
798 end if
799 end if
800
801 end subroutine change_date