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 end if
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 end if
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 end if
127
128 ! 3.3 limb scan check
129 !---------------------
130 if (check_limb) then
131 if ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCANX-1))) then
132 nsel(2) = nsel(2) + 1
133 CYCLE loop2
134 end if
135 end if
136
137 ! 3.4 Reject data outside radiosonde mask
138 !------------------------------------------
139 if (check_mask) then
140 call MASK(tovs%lat,tovs%lon,LMASK)
141 if (.NOT. LMASK) then
142 nsel(8) = nsel(8) + 1
143 CYCLE loop2
144 end if
145 end if
146
147 ! 3.5 Reject outliers : facx*sigma,
148 !--------------------------------------------------------------------------
149 do j=1, tovs%nchan
150 if ( (abs(vec_dep(j)-vmean_dep(j)) > (vstd_dep(j)*FAC)) ) then
151 tovs%qc_flag(j) = -1
152 end if
153 end do
154
155 do i=1,nchan
156 airbias(i) = xcoef0(i)
157 do j=1,npred
158 airbias(i) = airbias(i) + tovs%pred(j)*xcoef(i,j)
159 end do
160 vec_dep(i) = vec_dep(i) - airbias(i)
161 vec_abs(i) = vec_abs(i) - airbias(i)
162 if ( tovs%omb(i) == -888888.00 ) vec_dep(i)=tovs%omb(i)
163 end do
164
165 ! mean/std statistics for scan/airmass-bias corrected values
166 do j=1, nchan
167 if ( tovs%qc_flag(j) == 1 ) then
168 jv = j
169 nobs1(j) = nobs1(j) + 1
170 vmean_dep1(j) = vmean_dep1(j) + vec_dep(j)
171 vstd_dep1(j) = vstd_dep1(j) + vec_dep(j)*vec_dep(j)
172 vmean_abs1(j) = vmean_abs1(j) + vec_abs(j)
173 vstd_abs1(j) = vstd_abs1(j) + vec_abs(j)*vec_abs(j)
174 end if
175 end do
176
177 if (nchan == 15) then
178 write (11,'(15i15)') tovs%qc_flag(1:nchan)
179 write (11,'(15f15.3)') tovs%omb(1:nchan) ! omb no-biascorrection
180 write (11,'(15f15.3)') vec_dep(1:nchan) ! omb with bias correction
181 else if(nchan==5) then
182 write (11,'(5i15)') tovs%qc_flag(1:nchan)
183 write (11,'(5f15.3)') tovs%omb(1:nchan) ! omb no-biascorrection
184 write (11,'(5f15.3)') vec_dep(1:nchan) ! omb with bias correction
185 else if(nchan==19) then
186 write (11,'(19i15)') tovs%qc_flag(1:nchan)
187 write (11,'(19f15.3)') tovs%omb(1:nchan) ! omb no-biascorrection
188 write (11,'(19f15.3)') vec_dep(1:nchan) ! omb with bias correction
189 else
190 write (*,*) 'Unknow sensor'
191 end if
192
193 end do loop2
194
195 ! Calculate means, standard deviations and covariances
196
197 where (nobs1(:) /= 0)
198 vmean_dep1(:) = vmean_dep1(:)/nobs1(:)
199 vstd_dep1(:) = vstd_dep1(:)/nobs1(:) - vmean_dep1(:)**2
200 vstd_dep1(:) = SQRT(MAX(0.0,vstd_dep1(:)))
201 vmean_abs1(:) = vmean_abs1(:)/nobs1(:)
202 vstd_abs1(:) = vstd_abs1(:)/nobs1(:) - vmean_abs1(:)**2
203 vstd_abs1(:) = SQRT(MAX(0.0,vstd_abs1(:)))
204 end where
205
206 !----------------------------------------------------------------------------
207
208 deallocate(tovs%tb)
209 deallocate(tovs%omb)
210 deallocate(tovs%qc_flag)
211 deallocate(tovs%cloud_flag)
212 deallocate(tovs%pred)
213
214 ! out coefs to ASCII file bcor.asc
215 call write_biascoef(nchan,nscan,nband,npred,global, &
216 VMNRL(1:nchan,1:nscan),VMNRLB(1:nchan,1:nscan,1:nband), &
217 xcoef(1:nchan,1:npred),xcoef0(1:nchan), &
218 nobs1(1:nchan),vmean_abs1(1:nchan), &
219 vstd_abs1(1:nchan),vmean_dep1(1:nchan), &
220 vstd_dep1(1:nchan) )
221
222 end program da_bias_verif