da_bias_scan.f90
References to this file elsewhere.
1 program da_bias_scan
2
3 use rad_bias, only : long, bias, jpchan, jpscan, da_read_biasprep, print_bias, &
4 qc_amsua, qc_amsub, jband, boff, bdiv
5
6 ! PURPOSE.
7 ! -------
8 ! TO CALCULATE SCAN BIAS CORRECTIONS.
9
10 ! REFERENCE.
11 ! ----------
12 ! The ECMWF Tovs bias Correction, Harris and Kelly..
13
14 implicit none
15
16 ! Local Variables
17
18 type (bias) :: tovs
19
20 integer :: nobs(JPCHAN)
21 integer :: nscan(JPSCAN), nscanch(JPCHAN,JPSCAN)
22 integer :: nscanchb(JPCHAN,JPSCAN,JBAND)
23
24 real (kind=long) :: vmean_dep(JPCHAN), vmean_abs(JPCHAN)
25 real (kind=long) :: vstd_dep(JPCHAN), vstd_abs(JPCHAN)
26 real (kind=long) :: vec_dep(JPCHAN), vec_abs(JPCHAN)
27
28 real (kind=long) :: vmn(JPCHAN), vstd(JPCHAN)
29 real (kind=long) :: vmnsc(JPCHAN,JPSCAN), vstdsc(JPCHAN,JPSCAN), vmnrl(JPCHAN,JPSCAN)
30
31 real (kind=long) :: vmnscb(JPCHAN,JPSCAN,JBAND), vstdscb(JPCHAN,JPSCAN,JBAND)
32 real (kind=long) :: vmnrlb(JPCHAN,JPSCAN,JBAND)
33
34 ! smoothed values
35 integer :: nscanchbt(JPCHAN,JPSCAN,JBAND)
36 real (kind=long) :: vmnscbt(JPCHAN,JPSCAN,JBAND), vstdscbt(JPCHAN,JPSCAN,JBAND)
37 real (kind=long) :: vmnrlbt(JPCHAN,JPSCAN,JBAND)
38
39 ! integer :: ibin
40 integer :: ierr, iscanm
41 integer :: j, jv, jscan, js, k , l
42 ! integer :: path
43 ! integer :: cloud_flag
44
45 integer :: iband
46 ! integer :: ia
47 ! integer :: iab
48 ! integer :: aoff
49 ! integer :: avb
50
51 ! CHARACTER(LEN=2) :: cchan
52
53 real (kind=long) :: vmn0, vmn0b(JBAND)
54 ! real (kind=long) :: VDHL_WINdoW, VDLL_WINdoW
55
56 logical :: global, smoothing
57
58 ! definition of defaut namelist values
59 !-----------------------------------------
60 integer :: sband =1
61 integer :: eband =18
62 integer :: kscan =90
63 integer :: ICOUNT = 0 ! connter
64 real (kind=long) :: fac = 3.0 ! Sdev factor
65 namelist /INPUTS/ kscan,fac,global,sband,eband,smoothing
66
67 !------------------------------------------------------------------------------
68 ! 1. SETUP.
69 ! -----
70
71 read (5,INPUTS,end=100)
72
73 100 continue
74
75 write (6,INPUTS)
76
77 OPEN(unit=10,FORM='UNformatTED') ! Open input file : Input data from select
78 OPEN(unit=11,FORM='UNformatTED') ! Open output file : Scan mean core arrays
79 OPEN(unit=12,FORM='UNformatTED') ! Open output file : Scan bias coef
80
81
82 !------------------------------------------------------------------------------
83 ! 2. READ IN DATA, Q.C., CALC MEANS AND VARIANCES.
84 ! ---- -- ----- ----- ---- ----- --- ---------
85
86 vmean_dep = 0.0 ! Clear matrices
87 vstd_dep = 0.0
88 vmean_abs = 0.0 ! Clear matrices
89 vstd_abs = 0.0
90 nobs = 0
91
92 read (unit=10,end=265) tovs%nchan, tovs%npred ! Read in data
93 rewind(unit=10)
94
95 allocate(tovs%tb(tovs%nchan))
96 allocate(tovs%omb(tovs%nchan))
97 allocate(tovs%bias(tovs%nchan))
98 allocate(tovs%qc_flag(tovs%nchan))
99 allocate(tovs%cloud_flag(tovs%nchan))
100 allocate(tovs%pred(tovs%npred))
101
102 loop1:&
103 do
104 icount=icount+1
105
106 !
107 ! 2.1 read in data
108 !
109 call da_read_biasprep(tovs,10,ierr)
110 if (ierr == 0) then ! not end
111 continue
112 elseif (ierr == 1) then ! end
113 exit
114 else ! error
115 stop 'read error in da_scan_bias'
116 end if
117
118 if (icount < 2) call print_bias(tovs)
119
120 vec_dep(1:tovs%nchan) = tovs%omb(1:tovs%nchan) ! obs minus guess
121 vec_abs(1:tovs%nchan) = tovs%tb(1:tovs%nchan) ! obs
122
123 !
124 ! 2.2 QC: extrme values/extrme departure/window channel
125 !
126 if (tovs%sensor_id == 3) then
127 call qc_amsua(tovs)
128 elseif(tovs%sensor_id == 4) then
129 call qc_amsub(tovs)
130 end if
131
132 do j=1,tovs%nchan
133 !-------------------------------
134 ! 2.3 compute the statistics
135 !-------------------------------
136 if ( tovs%qc_flag(j) == 1 ) then
137 nobs(j) = nobs(j) + 1
138 vmean_dep(j) = vmean_dep(j) + vec_dep(j)
139 vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
140 vmean_abs(j) = vmean_abs(j) + vec_abs(j)
141 vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
142 end if
143 end do
144
145 end do loop1
146
147 265 continue ! MEANS AND STANDARD DEVIATIONS
148
149 write (6,270) icount,nobs(1:tovs%nchan)
150 270 format (/1X,'NUMBER OF DATA Total/ACCEPTED'/1X,i10/1X,15I10)
151
152 !---------------------------------
153 ! 2.8 mean and std for channels
154 !---------------------------------
155 where (nobs(:) /= 0)
156 vmean_dep(:) = vmean_dep(:)/nobs(:)
157 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
158 vstd_dep(:) = SQRT(MAX(0.0,vstd_dep(:)))
159 vmean_abs(:) = vmean_abs(:)/nobs(:)
160 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
161 vstd_abs(:) = SQRT(MAX(0.0,vstd_abs(:)))
162 end where
163
164 vmn(:) = vmean_dep(:) ! used by outlier check later
165 vstd(:) = vstd_dep(:)
166
167 write (6,288)
168 288 format (/1X,'FIRST PASS: MEANS AND STANDARD DEVIATIONS')
169
170 do j=1,tovs%nchan
171 jv = j
172 write (6,289) jv, nobs(j), vmean_abs(j),vstd_abs(j), vmean_dep(j),vstd_dep(j)
173 end do
174 289 format (1X,I5,I10,4F15.2)
175
176 !-------------------------------------------------------------------------------
177 ! 3. SECOND PASS THROUGH DATA, EXTRA Q.C.
178 ! ------ ---- ------- ----- ----- ----
179
180 rewind(unit=10)
181
182 vmean_dep = 0.0 ! Clear matrices
183 vstd_dep = 0.0
184 vmean_abs = 0.0 ! Clear matrices
185 vstd_abs = 0.0
186 vmnsc = 0.0
187 vstdsc = 0.0
188 vmnscb = 0.0
189 vstdscb = 0.0
190 nobs = 0
191 nscan = 0
192 nscanch = 0
193 nscanchb = 0
194
195 loop2:&
196 do
197
198 call da_read_biasprep(tovs,10,ierr)
199 if (ierr == 0) then ! not end
200 continue
201 elseif (ierr == 1) then ! end
202 exit
203 else ! error
204 stop 'read error in da_scan_bias'
205 end if
206
207 vec_dep(1:tovs%nchan) = tovs%omb(1:tovs%nchan)
208 vec_abs(1:tovs%nchan) = tovs%tb(1:tovs%nchan)
209
210 if (tovs%sensor_id == 3) then
211 call qc_amsua(tovs)
212 elseif(tovs%sensor_id == 4) then
213 call qc_amsub(tovs)
214 end if
215
216 !---------------------------------------------------------------------------------
217 ! 3.6 Reject outliers : using fac and stdv calculated in the first pass : loop 1
218 !---------------------------------------------------------------------------------
219 do j=1, tovs%nchan
220 if ( (abs(vec_dep(j)-vmn(j)) > (vstd(j)*FAC)) ) then
221 tovs%qc_flag(j) = -1
222 end if
223 end do
224
225 ! Good data : count and use in the ststistics
226 !----------------------------------------------------
227
228 jscan = tovs%scanpos
229 nscan(jscan) = nscan(jscan) + 1
230
231 iband = FLOOR(tovs%lat/BDIV) + BOFF !! latitude band
232
233 do j=1, tovs%nchan
234 if ( tovs%qc_flag(j) == 1 ) then
235
236 ! statistics for channel
237 !--------------------------
238 nobs(j) = nobs(j) + 1
239 vmean_dep(j) = vmean_dep(j) + vec_dep(j)
240 vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
241 vmean_abs(j) = vmean_abs(j) + vec_abs(j)
242 vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
243
244 ! statistics for channel+scan position
245 !----------------------------------------
246 nscanch(j,jscan) = nscanch(j,jscan) + 1
247 vmnsc(j,jscan) = vmnsc(j,jscan) + vec_dep(j)
248 vstdsc(j,jscan) = vstdsc(j,jscan) + vec_dep(j)*vec_dep(j)
249
250 ! statistics for channel + scan position + lat band
251 !----------------------------------------------------------
252 nscanchb(j,jscan,iband) = nscanchb(j,jscan,iband) + 1
253 vmnscb(j,jscan,iband) = vmnscb(j,jscan,iband) + vec_dep(j)
254 vstdscb(j,jscan,iband) = vstdscb(j,jscan,iband) + vec_dep(j)*vec_dep(j)
255
256 end if
257 end do
258
259 end do loop2
260
261 ! Write scan 'core' arrays NOTE: not divided by nobs
262 !------------------------------------------------------
263 write (11) nobs(:) ! statistics for channel
264 write (11) vmean_dep(:)
265 write (11) vstd_dep(:)
266 write (11) vmean_abs(:)
267 write (11) vstd_abs(:)
268
269 write (11) nscanch(:,:) ! statistics for channel + scan position
270 write (11) vmnsc(:,:)
271 write (11) vstdsc(:,:)
272
273 write (11) nscanchb(:,:,:) ! statistics for channel + scan position + lat band
274 write (11) vmnscb(:,:,:)
275 write (11) vstdscb(:,:,:)
276
277 ! Write out copious amounts of output
278 !----------------------------------------
279 write (6,270) icount,nobs(1:tovs%nchan)
280
281 where (nobs(:) /= 0) ! for channels
282 vmean_dep(:) = vmean_dep(:)/nobs(:)
283 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
284 vstd_dep(:) = SQRT(MAX(0.0,vstd_dep(:)))
285 vmean_abs(:) = vmean_abs(:)/nobs(:)
286 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
287 vstd_abs(:) = SQRT(MAX(0.0,vstd_abs(:)))
288 end where
289
290 do j=1, tovs%nchan
291 where (nscanch(j,:) /= 0) ! for channels + scan
292 vmnsc(j,:) = vmnsc(j,:)/nscanch(j,:)
293 vstdsc(j,:) = vstdsc(j,:)/nscanch(j,:) - vmnsc(j,:)**2
294 vstdsc(j,:) = SQRT(MAX(0.0,vstdsc(j,:)))
295 end where
296 where (nscanchb(j,:,:) /= 0) ! for channels + scan + band
297 vmnscb(j,:,:) = vmnscb(j,:,:)/nscanchb(j,:,:)
298 vstdscb(j,:,:) = vstdscb(j,:,:)/nscanchb(j,:,:) - vmnscb(j,:,:)**2
299 vstdscb(j,:,:) = SQRT(MAX(0.0,vstdscb(j,:,:)))
300 end where
301
302 ! 4.1 compute central scan position mean
303 !--------------------------------------------
304 if (MOD(KSCAN,2) == 0 ) then ! even scan number
305 iscanm = KSCAN/2 !! middle scan position
306
307 if( nscanch(j,iscanm) .ne. 0 .and. nscanch(j,iscanm+1) .ne.0 )then
308 vmn0 = 0.5*(vmnsc(j,iscanm)+vmnsc(j,iscanm+1))
309 vmn0b(:) = 0.5*(vmnscb(j,iscanm,:)+vmnscb(j,iscanm+1,:))
310 end if
311
312 if( nscanch(j,iscanm) .eq. 0 .and. nscanch(j,iscanm+1) .ne.0 )then
313 vmn0 = vmnsc(j,iscanm+1)
314 vmn0b(:) = vmnscb(j,iscanm+1,:)
315 end if
316
317 if( nscanch(j,iscanm) .ne. 0 .and. nscanch(j,iscanm+1) .eq.0 )then
318 vmn0 = vmnsc(j,iscanm)
319 vmn0b(:) = vmnscb(j,iscanm,:)
320 end if
321
322 ELSE
323 iscanm = kscan/2 + 1
324 vmn0 = vmnsc(j,iscanm)
325 vmn0b(:) = vmnscb(j,iscanm,:)
326 end if
327
328 ! 4.2 compute relative bias
329 !------------------------------------
330 do k=1,KSCAN
331 if( nscanch(j,k) .ne. 0 )then
332 vmnrl(j,k)=vmnsc(j,k) - vmn0
333 do l=1, JBAND
334 vmnrlb(j,K,l) = vmnscb(j,K,l) - vmn0b(l)
335 end do
336 end if
337 end do
338
339 end do
340
341 ! prinit output
342
343 !------------------
344 write (6,388)
345 388 format (/1X,'SECOND PASS: MEANS AND STANDARD DEVIATIONS')
346 do j=1, tovs%nchan
347 jv = j
348 write (6,289) jv, nobs(j), vmean_abs(j), vstd_abs(j), vmean_dep(j), vstd_dep(j)
349 end do
350
351 !-------------------
352 write (6,370) (js,js=1,KSCAN)
353 370 format (/5X,'NUMBER AT EACH SCAN POSITION'/5X,30I7)
354
355 do j=1, tovs%nchan
356 jv = j
357 write (6,371) jv, nscanch(j,1:KSCAN)
358 end do
359 371 format(1X,I4,30I7)
360
361 write (6,391) (js,js=1,KSCAN)
362 391 format (/1X,'BIASES FOR EACH SCAN ANGLE'/4X,30I7)
363 do j=1, tovs%nchan
364 jv = j
365 write (6,393) jv, (vmnsc(j,js),js=1,KSCAN)
366 end do
367 393 format (1X,I3,30F7.2)
368
369 write (6,394) (js,js=1,KSCAN)
370 394 format (/1X,'STD DEV FOR EACH SCAN ANGLE'/4X,30I7)
371 do j=1, tovs%nchan
372 jv = j
373 write (6,396) jv, (vstdsc(j,js),js=1,KSCAN)
374 end do
375 396 format (1X,I3,30F7.2)
376
377 !----------------------
378 if (global) then
379 write (6,'(1x,a17,2i4)') 'SCAN COEFFICIENTS', jband,kscan
380 do j=1,tovs%nchan
381 jv = j
382 do iband=1, jband
383 write (6,'(2i3,30f7.2)') jv, iband, vmnrlb(j,1:KSCAN,iband)
384 end do
385 end do
386
387 if (smoothing) then
388
389 close(11)
390
391 OPEN(11,FORM='UNformatTED')
392
393 read (11) nobs(:)
394 read (11) vmean_dep(:)
395 read (11) vstd_dep(:)
396 read (11) vmean_abs(:)
397 read (11) vstd_abs(:)
398
399 read (11) nscanch(:,:)
400 read (11) vmnsc(:,:)
401 read (11) vstdsc(:,:)
402
403 read (11) nscanchb(:,:,:)
404 read (11) vmnscb(:,:,:)
405 read (11) vstdscb(:,:,:)
406
407 ! Perform smoothing on banded corrections.
408 ! for absolute bias/stdev
409
410 vmnscbt = 0.0
411 vstdscbt= 0.0
412 nscanchbt = 0
413
414 do j=1, tovs%nchan
415 do iband=1, JBAND
416
417 ! if (iband <= AVBAND) then
418 ! aoff = AVBAND - iband + 1
419 ! avb = AVBAND
420 ! ELSEif ((iband < 2*AVBAND+1) .AND. (iband > AVBAND)) then
421 ! aoff = 0
422 ! avb = 2*AVBAND+1 - iband
423 ! ELSEif ((iband >= 2*AVBAND+1) .AND. (JBAND-iband+1 >= 2*AVBAND+1)) then
424 ! aoff = 0
425 ! avb = 0
426 ! ELSEif ((JBAND-iband+1 < 2*AVBAND+1) .AND. (JBAND-iband+1 > AVBAND)) then
427 ! aoff = 0
428 ! avb = 2*AVBAND+1 - (JBAND-iband+1)
429 ! ELSEif (JBAND-iband+1 <= AVBAND) then
430 ! aoff = -(AVBAND - (JBAND-iband+1) + 1)
431 ! avb = AVBAND
432 ! end if
433
434 ! do ia=-avb, avb
435 ! iab = iband + ia + aoff
436 ! vmnscbt(j,1:KSCAN,iband) = vmnscbt(j,1:KSCAN,iband) + vmnscb(j,1:KSCAN,IAb)
437 ! vstdscbt(j,1:KSCAN,iband) = vstdscbt(j,1:KSCAN,iband) + vstdscb(j,1:KSCAN,IAb)
438 ! nscanchbt(j,1:KSCAN,iband) = nscanchbt(j,1:KSCAN,iband) + nscanchb(j,1:KSCAN,IAb)
439 ! end do
440 vmnscbt(j,1:KSCAN,iband) = vmnscb(j,1:KSCAN,iband)
441 vstdscbt(j,1:KSCAN,iband) = vstdscb(j,1:KSCAN,iband)
442 nscanchbt(j,1:KSCAN,iband) = nscanchb(j,1:KSCAN,iband)
443 end do
444 end do
445
446 do j=1, tovs%nchan
447 where (nscanchbt(j,1:KSCAN,:) /= 0)
448 vmnscbt(j,1:KSCAN,:) = vmnscbt(j,1:KSCAN,:)/nscanchbt(j,1:KSCAN,:)
449 vstdscbt(j,1:KSCAN,:) = vstdscbt(j,1:KSCAN,:)/nscanchbt(j,1:KSCAN,:) - vmnscbt(j,1:KSCAN,:)**2
450 vstdscbt(j,1:KSCAN,:) = SQRT(MAX(0.0,vstdscbt(j,1:KSCAN,:)))
451 end where
452
453 ! get bias at nadir
454 if (MOD(KSCAN,2) == 0 ) then ! even scan number
455 iscanm = KSCAN/2
456
457 if( nscanch(j,iscanm) .ne. 0 .and. nscanch(j,iscanm+1) .ne.0 )then
458 vmn0b(:) = 0.5*(vmnscbt(j,iscanm,:)+vmnscbt(j,iscanm+1,:))
459 end if
460
461 if( nscanch(j,iscanm) .eq. 0 .and. nscanch(j,iscanm+1) .ne.0 )then
462 vmn0b(:) = vmnscbt(j,iscanm+1,:)
463 end if
464
465 if( nscanch(j,iscanm) .ne. 0 .and. nscanch(j,iscanm+1) .eq.0 )then
466 vmn0b(:) = vmnscbt(j,iscanm,:)
467 end if
468
469 ELSE
470 iscanm = kscan/2 + 1
471 vmn0b(:) = vmnscbt(j,iscanm,:)
472 end if
473
474 ! 4.2 compute relative bias
475 !------------------------------------
476 do k=1,KSCAN
477 if( nscanch(j,k) .ne. 0 )then
478 do l=1, JBAND
479 vmnrlb(j,k,l) = vmnscbt(j,k,l) - vmn0b(l)
480 end do
481 end if
482 end do
483
484 end do
485
486 ! Perform smoothing on banded corrections.
487 ! for relative bias
488
489 ! do j=1, tovs%nchan
490 ! vmnrlbt(j,1:KSCAN,1) = 0.50*vmnrlb(j,1:KSCAN,1) &
491 ! + 0.50*vmnrlb(j,1:KSCAN,2)
492 ! do iband=2, JBAND-1
493 ! vmnrlbt(j,1:KSCAN,iband) = 0.25*vmnrlb(j,1:KSCAN,iband-1) &
494 ! + 0.50*vmnrlb(j,1:KSCAN,iband) &
495 ! + 0.25*vmnrlb(j,1:KSCAN,iband+1)
496 ! end do
497 ! vmnrlbt(j,1:KSCAN,JBAND) = 0.50*vmnrlb(j,1:KSCAN,JBAND-1) &
498 ! + 0.50*vmnrlb(j,1:KSCAN,JBAND)
499 ! end do
500
501 if ( (eband-sband+1) >= 3 ) then
502 do j=1, tovs%nchan
503 vmnrlbt(j,1:KSCAN,sband) = 0.50*vmnrlb(j,1:KSCAN,sband) &
504 + 0.50*vmnrlb(j,1:KSCAN,sband+1)
505 do iband=sband, eband-1
506 vmnrlbt(j,1:KSCAN,iband) = 0.25*vmnrlb(j,1:KSCAN,iband-1) &
507 + 0.50*vmnrlb(j,1:KSCAN,iband) &
508 + 0.25*vmnrlb(j,1:KSCAN,iband+1)
509 end do
510 vmnrlbt(j,1:KSCAN,eband) = 0.50*vmnrlb(j,1:KSCAN,eband-1) &
511 + 0.50*vmnrlb(j,1:KSCAN,eband)
512 end do
513 elseif ( (eband-sband+1) == 2 ) then
514 vmnrlbt(j,1:KSCAN,sband) = 0.50*vmnrlb(j,1:KSCAN,sband) &
515 + 0.50*vmnrlb(j,1:KSCAN,sband+1)
516 vmnrlbt(j,1:KSCAN,eband) = 0.50*vmnrlb(j,1:KSCAN,eband-1) &
517 + 0.50*vmnrlb(j,1:KSCAN,eband)
518 elseif ( (eband-sband+1) == 1 ) then
519 vmnrlbt(j,1:KSCAN,sband) = vmnrlb(j,1:KSCAN,sband)
520 vmnrlbt(j,1:KSCAN,eband) = vmnrlb(j,1:KSCAN,eband)
521 end if
522
523 ! output relative bias
524 do j=1, tovs%nchan
525 jv = j
526 do iband=1,JBAND
527 write (12) jv, vmnrlbt(j,1:KSCAN,iband)
528 end do
529 end do
530
531 !----------------------
532 write (6,'(1x,a30,2i4)') 'SMOOTHED SCAN COEFFICIENTS', jband,kscan
533 do j=1,tovs%nchan
534 jv = j
535 do iband=1, jband
536 write (6,'(2i3,30f7.2)') jv, iband, vmnrlbt(j,1:KSCAN,iband)
537 end do
538 end do
539
540 end if ! end if smoothing
541
542 else ! regional
543 !---------------------
544 write (6,397) (js,js=1,KSCAN)
545 397 format (/1X,'RELATIVE BIASES FOR EACH SCAN ANGLE'/4X,30I7)
546 do j=1, tovs%nchan
547 jv = j
548 write (6,399) jv, vmnrl(j,1:KSCAN)
549 end do
550 399 format (1X,I3,30F7.2)
551
552 do j=1, tovs%nchan
553 jv = j
554 write (12) jv, vmnrl(j,1:KSCAN)
555 end do
556
557 !----------------------
558 end if ! global
559
560 deallocate(tovs%tb)
561 deallocate(tovs%omb)
562 deallocate(tovs%bias)
563 deallocate(tovs%qc_flag)
564 deallocate(tovs%cloud_flag)
565 deallocate(tovs%pred)
566
567 close(unit=10)
568 close(unit=11)
569 close(unit=12)
570
571 end program da_bias_scan