3 use rad_bias
, only
: long
, bias
, jpchan
, jpscan
, da_read_biasprep
, print_bias
, &
4 qc_amsua
, qc_amsub
, qc_ssmis
, jband
, boff
, bdiv
8 ! TO CALCULATE SCAN BIAS CORRECTIONS.
12 ! The ECMWF Tovs bias Correction, Harris and Kelly..
20 integer :: nobs(JPCHAN
)
21 integer :: nscan(JPSCAN
), nscanch(JPCHAN
,JPSCAN
)
22 integer :: nscanchb(JPCHAN
,JPSCAN
,JBAND
)
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
)
28 real (kind
=long
) :: vmn(JPCHAN
), vstd(JPCHAN
)
29 real (kind
=long
) :: vmnsc(JPCHAN
,JPSCAN
), vstdsc(JPCHAN
,JPSCAN
), vmnrl(JPCHAN
,JPSCAN
)
31 real (kind
=long
) :: vmnscb(JPCHAN
,JPSCAN
,JBAND
), vstdscb(JPCHAN
,JPSCAN
,JBAND
)
32 real (kind
=long
) :: vmnrlb(JPCHAN
,JPSCAN
,JBAND
)
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
)
40 integer :: ierr
, iscanm
41 integer :: j
, jv
, jscan
, js
, k
, l
43 ! integer :: cloud_flag
51 ! CHARACTER(LEN=2) :: cchan
53 real (kind
=long
) :: vmn0
, vmn0b(JBAND
)
54 ! real (kind=long) :: VDHL_WINdoW, VDLL_WINdoW
56 logical :: global
, smoothing
58 ! definition of defaut namelist values
59 !-----------------------------------------
63 integer :: ICOUNT
= 0 ! connter
64 real (kind
=long
) :: fac
= 3.0 ! Sdev factor
65 namelist /INPUTS
/ kscan
,fac
,global
,sband
,eband
,smoothing
67 !------------------------------------------------------------------------------
71 read (5,INPUTS
,end=100)
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
82 !------------------------------------------------------------------------------
83 ! 2. READ IN DATA, Q.C., CALC MEANS AND VARIANCES.
84 ! ---- -- ----- ----- ---- ----- --- ---------
86 vmean_dep
= 0.0 ! Clear matrices
88 vmean_abs
= 0.0 ! Clear matrices
92 read (unit
=10,end=265) tovs
%nchan
, tovs
%npred
! Read in data
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
))
109 call da_read_biasprep(tovs
,10,ierr
)
110 if (ierr
== 0) then ! not end
112 elseif (ierr
== 1) then ! end
115 stop 'read error in da_scan_bias'
118 if (icount
< 2) call print_bias(tovs
)
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
124 ! 2.2 QC: extrme values/extrme departure/window channel
126 if (tovs
%sensor_id
== 3) then
128 elseif(tovs
%sensor_id
== 4) then
130 elseif(tovs
%sensor_id
== 15) then ! mhs
132 elseif(tovs
%sensor_id
== 10) then ! ssmis
137 !-------------------------------
138 ! 2.3 compute the statistics
139 !-------------------------------
140 if ( tovs
%qc_flag(j
) == 1 ) then
141 nobs(j
) = nobs(j
) + 1
142 vmean_dep(j
) = vmean_dep(j
) + vec_dep(j
)
143 vstd_dep(j
) = vstd_dep(j
) + vec_dep(j
)*vec_dep(j
)
144 vmean_abs(j
) = vmean_abs(j
) + vec_abs(j
)
145 vstd_abs(j
) = vstd_abs(j
) + vec_abs(j
)*vec_abs(j
)
151 265 continue ! MEANS AND STANDARD DEVIATIONS
153 write (6,270) icount
,nobs(1:tovs
%nchan
)
154 270 format (/1X
,'NUMBER OF DATA Total/ACCEPTED'/1X
,i10
/1X
,15I10
)
156 !---------------------------------
157 ! 2.8 mean and std for channels
158 !---------------------------------
160 vmean_dep(:) = vmean_dep(:)/nobs(:)
161 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
162 vstd_dep(:) = SQRT(MAX(0.0_8
,vstd_dep(:)))
163 vmean_abs(:) = vmean_abs(:)/nobs(:)
164 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
165 vstd_abs(:) = SQRT(MAX(0.0_8
,vstd_abs(:)))
168 vmn(:) = vmean_dep(:) ! used by outlier check later
169 vstd(:) = vstd_dep(:)
172 288 format (/1X
,'FIRST PASS: MEANS AND STANDARD DEVIATIONS')
176 write (6,289) jv
, nobs(j
), vmean_abs(j
),vstd_abs(j
), vmean_dep(j
),vstd_dep(j
)
178 289 format (1X
,I5
,I10
,4F15
.2
)
180 !-------------------------------------------------------------------------------
181 ! 3. SECOND PASS THROUGH DATA, EXTRA Q.C.
182 ! ------ ---- ------- ----- ----- ----
186 vmean_dep
= 0.0 ! Clear matrices
188 vmean_abs
= 0.0 ! Clear matrices
202 call da_read_biasprep(tovs
,10,ierr
)
203 if (ierr
== 0) then ! not end
205 elseif (ierr
== 1) then ! end
208 stop 'read error in da_scan_bias'
211 vec_dep(1:tovs
%nchan
) = tovs
%omb(1:tovs
%nchan
)
212 vec_abs(1:tovs
%nchan
) = tovs
%tb(1:tovs
%nchan
)
214 if (tovs
%sensor_id
== 3) then
216 elseif(tovs
%sensor_id
== 4) then
218 elseif(tovs
%sensor_id
== 15) then ! mhs
220 elseif(tovs
%sensor_id
== 10) then ! ssmis
224 !---------------------------------------------------------------------------------
225 ! 3.6 Reject outliers : using fac and stdv calculated in the first pass : loop 1
226 !---------------------------------------------------------------------------------
228 if ( (abs(vec_dep(j
)-vmn(j
)) > (vstd(j
)*FAC
)) ) then
233 ! Good data : count and use in the ststistics
234 !----------------------------------------------------
237 nscan(jscan
) = nscan(jscan
) + 1
239 iband
= FLOOR(tovs
%lat
/BDIV
) + BOFF
!! latitude band
242 if ( tovs
%qc_flag(j
) == 1 ) then
244 ! statistics for channel
245 !--------------------------
246 nobs(j
) = nobs(j
) + 1
247 vmean_dep(j
) = vmean_dep(j
) + vec_dep(j
)
248 vstd_dep(j
) = vstd_dep(j
) + vec_dep(j
)*vec_dep(j
)
249 vmean_abs(j
) = vmean_abs(j
) + vec_abs(j
)
250 vstd_abs(j
) = vstd_abs(j
) + vec_abs(j
)*vec_abs(j
)
252 ! statistics for channel+scan position
253 !----------------------------------------
254 nscanch(j
,jscan
) = nscanch(j
,jscan
) + 1
255 vmnsc(j
,jscan
) = vmnsc(j
,jscan
) + vec_dep(j
)
256 vstdsc(j
,jscan
) = vstdsc(j
,jscan
) + vec_dep(j
)*vec_dep(j
)
258 ! statistics for channel + scan position + lat band
259 !----------------------------------------------------------
260 nscanchb(j
,jscan
,iband
) = nscanchb(j
,jscan
,iband
) + 1
261 vmnscb(j
,jscan
,iband
) = vmnscb(j
,jscan
,iband
) + vec_dep(j
)
262 vstdscb(j
,jscan
,iband
) = vstdscb(j
,jscan
,iband
) + vec_dep(j
)*vec_dep(j
)
269 ! Write scan 'core' arrays NOTE: not divided by nobs
270 !------------------------------------------------------
271 write (11) nobs(:) ! statistics for channel
272 write (11) vmean_dep(:)
273 write (11) vstd_dep(:)
274 write (11) vmean_abs(:)
275 write (11) vstd_abs(:)
277 write (11) nscanch(:,:) ! statistics for channel + scan position
278 write (11) vmnsc(:,:)
279 write (11) vstdsc(:,:)
281 write (11) nscanchb(:,:,:) ! statistics for channel + scan position + lat band
282 write (11) vmnscb(:,:,:)
283 write (11) vstdscb(:,:,:)
285 ! Write out copious amounts of output
286 !----------------------------------------
287 write (6,270) icount
,nobs(1:tovs
%nchan
)
289 where (nobs(:) /= 0) ! for channels
290 vmean_dep(:) = vmean_dep(:)/nobs(:)
291 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
292 vstd_dep(:) = SQRT(MAX(0.0_8
,vstd_dep(:)))
293 vmean_abs(:) = vmean_abs(:)/nobs(:)
294 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
295 vstd_abs(:) = SQRT(MAX(0.0_8
,vstd_abs(:)))
299 where (nscanch(j
,:) /= 0) ! for channels + scan
300 vmnsc(j
,:) = vmnsc(j
,:)/nscanch(j
,:)
301 vstdsc(j
,:) = vstdsc(j
,:)/nscanch(j
,:) - vmnsc(j
,:)**2
302 vstdsc(j
,:) = SQRT(MAX(0.0_8
,vstdsc(j
,:)))
304 where (nscanchb(j
,:,:) /= 0) ! for channels + scan + band
305 vmnscb(j
,:,:) = vmnscb(j
,:,:)/nscanchb(j
,:,:)
306 vstdscb(j
,:,:) = vstdscb(j
,:,:)/nscanchb(j
,:,:) - vmnscb(j
,:,:)**2
307 vstdscb(j
,:,:) = SQRT(MAX(0.0_8
,vstdscb(j
,:,:)))
310 ! 4.1 compute central scan position mean
311 !--------------------------------------------
312 if (MOD(KSCAN
,2) == 0 ) then ! even scan number
313 iscanm
= KSCAN
/2 !! middle scan position
315 if( nscanch(j
,iscanm
) .ne
. 0 .and
. nscanch(j
,iscanm
+1) .ne
.0 )then
316 vmn0
= 0.5*(vmnsc(j
,iscanm
)+vmnsc(j
,iscanm
+1))
317 vmn0b(:) = 0.5*(vmnscb(j
,iscanm
,:)+vmnscb(j
,iscanm
+1,:))
320 if( nscanch(j
,iscanm
) .eq
. 0 .and
. nscanch(j
,iscanm
+1) .ne
.0 )then
321 vmn0
= vmnsc(j
,iscanm
+1)
322 vmn0b(:) = vmnscb(j
,iscanm
+1,:)
325 if( nscanch(j
,iscanm
) .ne
. 0 .and
. nscanch(j
,iscanm
+1) .eq
.0 )then
326 vmn0
= vmnsc(j
,iscanm
)
327 vmn0b(:) = vmnscb(j
,iscanm
,:)
332 vmn0
= vmnsc(j
,iscanm
)
333 vmn0b(:) = vmnscb(j
,iscanm
,:)
336 ! 4.2 compute relative bias
337 !------------------------------------
339 if( nscanch(j
,k
) .ne
. 0 )then
340 vmnrl(j
,k
)=vmnsc(j
,k
) - vmn0
342 vmnrlb(j
,K
,l
) = vmnscb(j
,K
,l
) - vmn0b(l
)
353 388 format (/1X
,'SECOND PASS: MEANS AND STANDARD DEVIATIONS')
356 write (6,289) jv
, nobs(j
), vmean_abs(j
), vstd_abs(j
), vmean_dep(j
), vstd_dep(j
)
360 write (6,370) (js
,js
=1,KSCAN
)
361 370 format (/5X
,'NUMBER AT EACH SCAN POSITION'/5X
,90I7
)
365 write (6,371) jv
, nscanch(j
,1:KSCAN
)
367 371 format(1X
,I4
,90I7
)
369 write (6,391) (js
,js
=1,KSCAN
)
370 391 format (/1X
,'BIASES FOR EACH SCAN ANGLE'/4X
,90I7
)
373 write (6,393) jv
, (vmnsc(j
,js
),js
=1,KSCAN
)
375 393 format (1X
,I3
,90F7
.2
)
377 write (6,394) (js
,js
=1,KSCAN
)
378 394 format (/1X
,'STD DEV FOR EACH SCAN ANGLE'/4X
,90I7
)
381 write (6,396) jv
, (vstdsc(j
,js
),js
=1,KSCAN
)
383 396 format (1X
,I3
,90F7
.2
)
385 !----------------------
387 write (6,'(1x,a17,2i4)') 'SCAN COEFFICIENTS', jband
,kscan
391 write (6,'(2i3,90f7.2)') jv
, iband
, vmnrlb(j
,1:KSCAN
,iband
)
399 OPEN(11,FORM
='UNformatTED')
402 read (11) vmean_dep(:)
403 read (11) vstd_dep(:)
404 read (11) vmean_abs(:)
405 read (11) vstd_abs(:)
407 read (11) nscanch(:,:)
409 read (11) vstdsc(:,:)
411 read (11) nscanchb(:,:,:)
412 read (11) vmnscb(:,:,:)
413 read (11) vstdscb(:,:,:)
415 ! Perform smoothing on banded corrections.
416 ! for absolute bias/stdev
425 ! if (iband <= AVBAND) then
426 ! aoff = AVBAND - iband + 1
428 ! ELSEif ((iband < 2*AVBAND+1) .AND. (iband > AVBAND)) then
430 ! avb = 2*AVBAND+1 - iband
431 ! ELSEif ((iband >= 2*AVBAND+1) .AND. (JBAND-iband+1 >= 2*AVBAND+1)) then
434 ! ELSEif ((JBAND-iband+1 < 2*AVBAND+1) .AND. (JBAND-iband+1 > AVBAND)) then
436 ! avb = 2*AVBAND+1 - (JBAND-iband+1)
437 ! ELSEif (JBAND-iband+1 <= AVBAND) then
438 ! aoff = -(AVBAND - (JBAND-iband+1) + 1)
443 ! iab = iband + ia + aoff
444 ! vmnscbt(j,1:KSCAN,iband) = vmnscbt(j,1:KSCAN,iband) + vmnscb(j,1:KSCAN,IAb)
445 ! vstdscbt(j,1:KSCAN,iband) = vstdscbt(j,1:KSCAN,iband) + vstdscb(j,1:KSCAN,IAb)
446 ! nscanchbt(j,1:KSCAN,iband) = nscanchbt(j,1:KSCAN,iband) + nscanchb(j,1:KSCAN,IAb)
448 vmnscbt(j
,1:KSCAN
,iband
) = vmnscb(j
,1:KSCAN
,iband
)
449 vstdscbt(j
,1:KSCAN
,iband
) = vstdscb(j
,1:KSCAN
,iband
)
450 nscanchbt(j
,1:KSCAN
,iband
) = nscanchb(j
,1:KSCAN
,iband
)
455 where (nscanchbt(j
,1:KSCAN
,:) /= 0)
456 vmnscbt(j
,1:KSCAN
,:) = vmnscbt(j
,1:KSCAN
,:)/nscanchbt(j
,1:KSCAN
,:)
457 vstdscbt(j
,1:KSCAN
,:) = vstdscbt(j
,1:KSCAN
,:)/nscanchbt(j
,1:KSCAN
,:) - vmnscbt(j
,1:KSCAN
,:)**2
458 vstdscbt(j
,1:KSCAN
,:) = SQRT(MAX(0.0_8
,vstdscbt(j
,1:KSCAN
,:)))
462 if (MOD(KSCAN
,2) == 0 ) then ! even scan number
465 if( nscanch(j
,iscanm
) .ne
. 0 .and
. nscanch(j
,iscanm
+1) .ne
.0 )then
466 vmn0b(:) = 0.5*(vmnscbt(j
,iscanm
,:)+vmnscbt(j
,iscanm
+1,:))
469 if( nscanch(j
,iscanm
) .eq
. 0 .and
. nscanch(j
,iscanm
+1) .ne
.0 )then
470 vmn0b(:) = vmnscbt(j
,iscanm
+1,:)
473 if( nscanch(j
,iscanm
) .ne
. 0 .and
. nscanch(j
,iscanm
+1) .eq
.0 )then
474 vmn0b(:) = vmnscbt(j
,iscanm
,:)
479 vmn0b(:) = vmnscbt(j
,iscanm
,:)
482 ! 4.2 compute relative bias
483 !------------------------------------
485 if( nscanch(j
,k
) .ne
. 0 )then
487 vmnrlb(j
,k
,l
) = vmnscbt(j
,k
,l
) - vmn0b(l
)
494 ! Perform smoothing on banded corrections.
498 ! vmnrlbt(j,1:KSCAN,1) = 0.50*vmnrlb(j,1:KSCAN,1) &
499 ! + 0.50*vmnrlb(j,1:KSCAN,2)
500 ! do iband=2, JBAND-1
501 ! vmnrlbt(j,1:KSCAN,iband) = 0.25*vmnrlb(j,1:KSCAN,iband-1) &
502 ! + 0.50*vmnrlb(j,1:KSCAN,iband) &
503 ! + 0.25*vmnrlb(j,1:KSCAN,iband+1)
505 ! vmnrlbt(j,1:KSCAN,JBAND) = 0.50*vmnrlb(j,1:KSCAN,JBAND-1) &
506 ! + 0.50*vmnrlb(j,1:KSCAN,JBAND)
509 if ( (eband
-sband
+1) >= 3 ) then
511 vmnrlbt(j
,1:KSCAN
,sband
) = 0.50*vmnrlb(j
,1:KSCAN
,sband
) &
512 + 0.50*vmnrlb(j
,1:KSCAN
,sband
+1)
513 do iband
=sband
, eband
-1
514 vmnrlbt(j
,1:KSCAN
,iband
) = 0.25*vmnrlb(j
,1:KSCAN
,iband
-1) &
515 + 0.50*vmnrlb(j
,1:KSCAN
,iband
) &
516 + 0.25*vmnrlb(j
,1:KSCAN
,iband
+1)
518 vmnrlbt(j
,1:KSCAN
,eband
) = 0.50*vmnrlb(j
,1:KSCAN
,eband
-1) &
519 + 0.50*vmnrlb(j
,1:KSCAN
,eband
)
521 elseif ( (eband
-sband
+1) == 2 ) then
522 vmnrlbt(j
,1:KSCAN
,sband
) = 0.50*vmnrlb(j
,1:KSCAN
,sband
) &
523 + 0.50*vmnrlb(j
,1:KSCAN
,sband
+1)
524 vmnrlbt(j
,1:KSCAN
,eband
) = 0.50*vmnrlb(j
,1:KSCAN
,eband
-1) &
525 + 0.50*vmnrlb(j
,1:KSCAN
,eband
)
526 elseif ( (eband
-sband
+1) == 1 ) then
527 vmnrlbt(j
,1:KSCAN
,sband
) = vmnrlb(j
,1:KSCAN
,sband
)
528 vmnrlbt(j
,1:KSCAN
,eband
) = vmnrlb(j
,1:KSCAN
,eband
)
531 ! output relative bias
535 write (12) jv
, vmnrlbt(j
,1:KSCAN
,iband
)
539 !----------------------
540 write (6,'(1x,a30,2i4)') 'SMOOTHED SCAN COEFFICIENTS', jband
,kscan
544 write (6,'(2i3,90f7.2)') jv
, iband
, vmnrlbt(j
,1:KSCAN
,iband
)
548 end if ! end if smoothing
551 !---------------------
552 write (6,397) (js
,js
=1,KSCAN
)
553 397 format (/1X
,'RELATIVE BIASES FOR EACH SCAN ANGLE'/4X
,90I7
)
556 write (6,399) jv
, vmnrl(j
,1:KSCAN
)
558 399 format (1X
,I3
,90F7
.2
)
562 write (12) jv
, vmnrl(j
,1:KSCAN
)
565 !----------------------
570 deallocate(tovs
%bias
)
571 deallocate(tovs
%qc_flag
)
572 deallocate(tovs
%cloud_flag
)
573 deallocate(tovs
%pred
)
579 end program da_bias_scan