da_bias_verif.f90

References to this file elsewhere.
1   PROGRAM da_bias_verif
2 
3   USE RAD_BIAS, only : bias, long, jpchan, jpscan,jband, jpnx, &
4      da_read_biasprep, read_biascoef, get_scorr, qc_amsua, qc_amsub, mask, &
5      write_biascoef, jpny
6 
7 ! PURPOSE.
8 ! --------
9 ! Apply bias correction to radiance data
10 ! to statistically verify algorithm
11 
12   IMPLICIT NONE
13 
14   TYPE(BIAS) :: tovs
15 
16   REAL(KIND=LONG) :: vmean_dep(JPCHAN), vmean_abs(JPCHAN)
17   REAL(KIND=LONG) :: vec_dep(JPCHAN), vec_abs(JPCHAN)
18   REAL(KIND=LONG) :: vstd_dep(JPCHAN), vstd_abs(JPCHAN)
19   REAL(KIND=LONG) :: vmean_dep1(JPCHAN), vmean_abs1(JPCHAN)
20   REAL(KIND=LONG) :: vstd_dep1(JPCHAN), vstd_abs1(JPCHAN)
21 
22   REAL(KIND=LONG) :: airbias(JPCHAN)
23 
24   REAL(KIND=LONG) :: xcoef0(JPCHAN), xcoef(JPCHAN,JPNX)
25 
26   REAL(KIND=LONG) :: vmnrl(JPCHAN,JPSCAN) = 0.0
27   REAL(KIND=LONG) :: vmnrlb(JPCHAN,JPSCAN,JBAND) = 0.0
28 
29   REAL(KIND=LONG) :: SCORR(JPCHAN)
30 
31   INTEGER :: nsel(10)
32   INTEGER :: nobs(JPCHAN),nobs1(JPCHAN)
33 
34   LOGICAL :: LMASK
35 
36   LOGICAL :: lscan = .FALSE.
37 
38   INTEGER :: I, iband, ierr
39   INTEGER :: J, JSCAN, jv
40 
41   INTEGER :: kscanx=90
42   ! INTEGER :: jbandx
43   LOGICAL :: check_limb=.false., check_mask=.false., global
44   REAL    :: FAC = 3.0      ! Number of SD' for QC
45 
46   INTEGER :: nchan,nscan,nband,npred
47 
48   NAMELIST /INPUTS/ global,lscan, kscanx, check_limb, check_mask, FAC
49 
50 !------------------------------------------------------------------------------
51 !        1.   SETUP.
52 !             -----
53   READ(5,INPUTS,END=100)
54   100 CONTINUE
55   WRITE(6,INPUTS)
56 
57 ! 1.0 read bias correction coefficients
58 !   -------------------------------------------------------------
59      open(UNIT=12,file='scor.asc',form='formatted')
60      read (12,'(4i6)') nchan,nscan,nband,npred
61      close(12)
62 
63      call read_biascoef(nchan,nscan,nband,npred,global, &
64                       VMNRL(1:nchan,1:nscan),VMNRLB(1:nchan,1:nscan,1:nband), &
65                       xcoef(1:nchan,1:npred),xcoef0(1:nchan), &
66                       nobs(1:nchan),vmean_abs(1:nchan), &
67                       vstd_abs(1:nchan),vmean_dep(1:nchan), &
68                       vstd_dep(1:nchan) )
69 
70 !----------------------------------------------------------------------
71 !        2.   READ IN DATA, Q.C., CALC MEANS AND VARIANCES.
72 !             ---- -- ----- ----- ---- ----- --- ---------
73 
74   nobs1(:)  = 0
75   nsel(:)  = 0
76   vmean_dep1(:) = 0.0
77   vmean_abs1(:) = 0.0
78   vstd_dep1(:)  = 0.0
79   vstd_abs1(:)  = 0.0
80 
81   READ(UNIT=10,END=300)  tovs%nchan, tovs%npred    ! Read in data
82   REWIND(UNIT=10)
83 
84   allocate(tovs%tb(tovs%nchan))
85   allocate(tovs%omb(tovs%nchan))
86   allocate(tovs%bias(tovs%nchan))
87   allocate(tovs%qc_flag(tovs%nchan))
88   allocate(tovs%cloud_flag(tovs%nchan))
89   allocate(tovs%pred(tovs%npred))
90 
91   300 CONTINUE
92 
93 loop2:&
94   DO
95     call da_read_biasprep(tovs,10,ierr)
96     if (ierr == 0) then      ! not end
97          continue
98     elseif (ierr == 1) then  ! end
99          exit
100     else                     ! error
101          stop 'read error in da_veri_bias'
102     endif
103 
104 ! latitude band
105     iband = INT(tovs%lat/30.000001) + 3
106     IF (lscan) THEN
107         JSCAN = tovs%scanpos
108       if (global) then
109         CALL GET_SCORR(JPCHAN,SCORR(1:JPCHAN),tovs%lat,vmnrlb,JSCAN)
110       else
111         SCORR(1:nchan) = vmnrl(1:nchan,JSCAN)
112       end if
113     ELSE
114         JSCAN = KSCANX/2
115         SCORR(1:JPCHAN) = 0.0
116     ENDIF
117 
118     vec_dep(1:nchan) = tovs%omb(1:nchan) - SCORR(1:nchan)
119     vec_abs(1:nchan) =  tovs%tb(1:nchan) - SCORR(1:nchan)
120 
121 ! 3.2 QC:
122       if (tovs%sensor_id == 3) then
123            call qc_amsua(tovs)
124       elseif(tovs%sensor_id == 4) then
125            call qc_amsub(tovs)
126       elseif(tovs%sensor_id == 15) then
127            call qc_amsub(tovs)
128       endif
129 
130 ! 3.3 limb scan check
131 !---------------------
132    if (check_limb) then
133     IF ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCANX-1))) THEN
134       nsel(2) = nsel(2) + 1
135       CYCLE loop2
136     ENDIF
137    end if
138 
139 ! 3.4 Reject data outside radiosonde mask
140 !------------------------------------------
141     IF (check_mask) THEN
142       CALL MASK(tovs%lat,tovs%lon,LMASK)
143       IF (.NOT. LMASK) THEN
144         nsel(8) = nsel(8) + 1
145         CYCLE loop2
146       ENDIF
147     ENDIF
148 
149 ! 3.5 Reject outliers : facx*sigma, 
150 !--------------------------------------------------------------------------
151     DO j=1, tovs%nchan
152       IF ( (ABS(vec_dep(j)-vmean_dep(j)) > (vstd_dep(j)*FAC)) ) THEN
153         tovs%qc_flag(j) = -1
154       ENDIF
155     ENDDO
156 
157     do i=1,nchan 
158       airbias(i) = xcoef0(i)
159       do j=1,npred
160         airbias(i) = airbias(i) + tovs%pred(j)*xcoef(i,j)
161       end do
162       vec_dep(i) = vec_dep(i) - airbias(i)
163       vec_abs(i) = vec_abs(i) - airbias(i)
164       if ( tovs%omb(i) == -888888.00 ) vec_dep(i)=tovs%omb(i)
165     end do
166 
167 ! mean/std statistics for scan/airmass-bias corrected values                      
168     DO j=1, nchan
169       IF ( tovs%qc_flag(j) == 1 ) THEN
170         jv = j
171         nobs1(j) = nobs1(j) + 1
172         vmean_dep1(j) = vmean_dep1(j) + vec_dep(j)
173          vstd_dep1(j) = vstd_dep1(j) + vec_dep(j)*vec_dep(j)
174         vmean_abs1(j) = vmean_abs1(j) + vec_abs(j)
175          vstd_abs1(j) = vstd_abs1(j) + vec_abs(j)*vec_abs(j)
176       ENDIF
177     ENDDO
178 
179     if (nchan == 15) then
180       write(11,'(15i15)') tovs%qc_flag(1:nchan)
181       write(11,'(15f15.3)') tovs%omb(1:nchan)  ! omb no-biascorrection
182       write(11,'(15f15.3)') vec_dep(1:nchan)   ! omb with bias correction
183     else if(nchan==5) then
184       write(11,'(5i15)') tovs%qc_flag(1:nchan)
185       write(11,'(5f15.3)') tovs%omb(1:nchan)  ! omb no-biascorrection
186       write(11,'(5f15.3)') vec_dep(1:nchan)   ! omb with bias correction
187     else if(nchan==19) then
188       write(11,'(19i15)') tovs%qc_flag(1:nchan)
189       write(11,'(19f15.3)') tovs%omb(1:nchan)  ! omb no-biascorrection
190       write(11,'(19f15.3)') vec_dep(1:nchan)   ! omb with bias correction
191     else
192       write(*,*) 'Unknow sensor'
193     end if
194 
195   ENDDO loop2
196 
197 ! Calculate means, standard deviations and covariances
198 
199   WHERE (nobs1(:) /= 0)
200     vmean_dep1(:) = vmean_dep1(:)/nobs1(:)
201     vstd_dep1(:)  = vstd_dep1(:)/nobs1(:) - vmean_dep1(:)**2
202     vstd_dep1(:)  = SQRT(MAX(0.0,vstd_dep1(:)))
203     vmean_abs1(:) = vmean_abs1(:)/nobs1(:)
204     vstd_abs1(:)  = vstd_abs1(:)/nobs1(:) - vmean_abs1(:)**2
205     vstd_abs1(:)  = SQRT(MAX(0.0,vstd_abs1(:)))
206   ENDWHERE
207 
208 !----------------------------------------------------------------------------
209 
210    deallocate(tovs%tb)
211    deallocate(tovs%omb)
212    deallocate(tovs%qc_flag)
213    deallocate(tovs%cloud_flag)
214    deallocate(tovs%pred)
215 
216 ! out coefs to ASCII file bcor.asc
217   call write_biascoef(nchan,nscan,nband,npred,global, &
218                       VMNRL(1:nchan,1:nscan),VMNRLB(1:nchan,1:nscan,1:nband), &
219                       xcoef(1:nchan,1:npred),xcoef0(1:nchan), &
220                       nobs1(1:nchan),vmean_abs1(1:nchan), &
221                       vstd_abs1(1:nchan),vmean_dep1(1:nchan), &
222                       vstd_dep1(1:nchan) )
223 
224   END PROGRAM da_bias_verif