da_bias_airmass.f90
References to this file elsewhere.
1 program da_bias_airmass
2
3 use rad_bias
4
5 ! PURPOSE.
6 ! --------
7 ! To calculate spatially varying bias correction COEFFICIENTS
8 ! for TOVS using model bias predictors.
9
10 ! EXTERNALS.
11 ! ----------
12 ! REGRESS_ONE
13
14 ! REFERENCE.
15 ! ----------
16 ! Harris and Kelly 1998.
17
18 implicit none
19
20 type (bias) :: tovs
21
22 real (kind=long) :: vmean_dep(JPCHAN), vmean_abs(JPCHAN)
23 real (kind=long) :: vec_dep(JPCHAN), vec_abs(JPCHAN)
24 real (kind=long) :: vstd_dep(JPCHAN), vstd_abs(JPCHAN)
25 real (kind=long) :: vmn(JPCHAN), vstd(JPCHAN)
26
27 real :: pred(JPNX)
28
29 real (kind=long) :: xbar(JPNY,JPNX), ybar(JPNY)
30 real (kind=long) :: xcov(JPNY,JPNX,JPNX), ycov(JPNY), xycov(JPNY,JPNX)
31
32 real (kind=long) :: coef(JPNY,JPNX), coef0(JPNY)
33 real (kind=long) :: reserr(JPNY), rescov(JPNY,JPNY)
34 real (kind=long) :: xvec(JPNY,JPNX,JPNX)
35 real (kind=long) :: xreser(JPCHAN), xcoef0(JPCHAN), xcoef(JPCHAN,JPNX)
36
37 real (kind=long) :: vmnbd(JPNY,6), vstdbd(JPNY,6)
38 real (kind=long) :: vmnbdl(JPNY,6,0:1), vstdbdl(JPNY,6,0:1)
39 real (kind=long) :: vmnb0(JPNY,6), vstdb0(JPNY,6)
40 real (kind=long) :: vmnbl0(JPNY,6,0:1), vstdbl0(JPNY,6,0:1)
41 real (kind=long) :: dvmnbd(JPNY,6), dvstdb(JPNY,6)
42
43 real (kind=long) :: vmnrl(JPCHAN,JPSCAN) = 0.0
44 real (kind=long) :: vmnrlb(JPCHAN,JPSCAN,JBAND) = 0.0
45
46 real (kind=long) :: SCORR(JPCHAN)
47
48 integer :: nsel(10)
49 integer :: nobs(JPCHAN)
50 integer :: nobsy(JPNY),lchan
51
52 integer :: nband(6), nbandl(6,0:1), n_band_ch(6,JPNY,0:1) = 0
53 integer :: nband0(JPNY,6), nbandd(JPNY,6)
54 integer :: nbandl0(JPNY,6,0:1), nbanddl(JPNY,6,0:1)
55
56 logical :: LMASK
57
58 logical :: lscan = .FALSE., global = .FALSE.
59
60 integer :: IR, ibin, I, iband, II, IV, ib, ierr
61 integer :: JS, J, JJ, JCOUNT, JBOX, JMINI, JSCAN, jv, IIV, JJV, sband
62
63 real (kind=long) :: xcorr(JPCHAN)
64 real (kind=long) :: coef_sat, coef_year, coef_month, coef_day, coef_time
65 integer :: CDATE = 1998010105 ! Current Date
66
67 integer :: icount = 0
68 integer :: KSCAN = JPSCAN
69 logical :: check_limb=.false., check_mask=.false.
70 real :: FAC = 3.0 ! Number of SD' for QC
71
72
73 namelist /INPUTS/ global, lscan, kscan, check_limb, check_mask, &
74 FAC,CDATE
75
76 !------------------------------------------------------------------------------
77 ! 1. SETUP.
78 ! -----
79
80 read (5,INPUTS,end=100)
81 100 continue
82 write (6,INPUTS)
83
84 if (lscan) OPEN(12,FORM='UNformatTED') ! Scan Biases with lat bands
85 OPEN(unit=10,FORM='UNformatTED') ! Input data from SELECT
86
87 ! Read scan biases
88
89 if (lscan) then
90 write (6,112) (JS,JS=1,KSCAN)
91 112 format(/1X,'SCAN DEPendENT BIASES APPLIED'/1X,3X,18I7/(4X,18I7))
92 ELSE
93 write (6,*) 'NO SCAN CORRECTION APPLIED'
94 end if
95
96 !----------------------------------------------------------------------
97 ! 2. READ IN DATA, Q.C., CALC MEANS AND VARIANCES.
98 ! ---- -- ----- ----- ---- ----- --- ---------
99
100 nobs(:) = 0
101 nsel(:) = 0
102 vmean_dep(:) = 0.0
103 vmean_abs(:) = 0.0
104 vstd_dep(:) = 0.0
105 vstd_abs(:) = 0.0
106
107 200 continue
108
109 read (unit=10,end=265) tovs%nchan, tovs%npred ! Read in data
110 rewind(unit=10)
111
112 allocate(tovs%tb(tovs%nchan))
113 allocate(tovs%omb(tovs%nchan))
114 allocate(tovs%bias(tovs%nchan))
115 allocate(tovs%qc_flag(tovs%nchan))
116 allocate(tovs%cloud_flag(tovs%nchan))
117 allocate(tovs%pred(tovs%npred))
118
119 if (lscan) then
120 if (global) then
121 do J=1, tovs%nchan ! get scan biases (latitude bands)
122 do sband=1,JBAND
123 read (12) JJ, vmnrlb(J,1:KSCAN,sband)
124 end do
125 end do
126 else
127 do J=1, tovs%nchan ! get scan biases (no band dependence)
128 read (12) JJ, vmnrl(J,1:KSCAN)
129 end do
130 end if
131 end if
132
133 loop1:&
134 do
135 icount=icount+1
136
137 call da_read_biasprep(tovs,10,ierr)
138 if (ierr == 0) then ! not end
139 continue
140 elseif (ierr == 1) then ! end
141 exit
142 else ! error
143 stop 'read error in da_airmass_bias'
144 end if
145
146 ! select lat band and get scan bias
147 iband = INT(tovs%lat/30.000001) + 3
148 if (lscan) then
149 JSCAN = tovs%scanpos
150 if (global) then
151 call GET_SCORR(JPCHAN,SCORR(1:JPCHAN),tovs%lat,vmnrlb,JSCAN)
152 else
153 SCORR(1:tovs%nchan) = vmnrl(1:tovs%nchan,JSCAN)
154 end if
155 ELSE
156 JSCAN = KSCAN/2
157 SCORR(1:JPCHAN) = 0.0
158 end if
159
160 ! apply scan bias to the departure (from FG or analysis) and the TB
161 vec_dep(1:tovs%nchan) = tovs%omb(1:tovs%nchan) - SCORR(1:tovs%nchan)
162 vec_abs(1:tovs%nchan) = tovs%tb(1:tovs%nchan) - SCORR(1:tovs%nchan)
163
164 !
165 ! 2.2 QC: extrme values/extrme departure/window channel
166 !
167 if (tovs%sensor_id == 3) then
168 call qc_amsua(tovs)
169 elseif(tovs%sensor_id == 4) then
170 call qc_amsub(tovs)
171 end if
172
173 !-------------------------
174 ! 2.3 limb sounding check
175 !-------------------------
176 if (check_limb) then
177 if ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCAN-1))) then ! Reject edge of scan
178 nsel(2) = nsel(2) + 1
179 CYCLE loop1
180 end if
181 end if
182
183 !-----------------------------------------
184 ! 2.4 Reject data outside radiosonde mask
185 !-----------------------------------------
186 if (check_mask) then
187 call MASK(tovs%lat,tovs%lon,LMASK)
188 if (.NOT. LMASK) then
189 nsel(8) = nsel(8) + 1
190 CYCLE loop1
191 end if
192 end if
193
194 ! Good data : count and use in the statistics
195
196 do j=1, tovs%nchan
197 if ( tovs%qc_flag(j) == 1 ) then
198 ! statistics for channel
199 !--------------------------
200 nobs(j) = nobs(j) + 1
201 vmean_dep(j) = vmean_dep(j) + vec_dep(j)
202 vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
203 vmean_abs(j) = vmean_abs(j) + vec_abs(j)
204 vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
205 end if
206 end do
207
208 end do loop1
209
210 265 continue
211
212 !---------------------------------
213 ! 2.8 mean and std for channels
214 !---------------------------------
215 where (nobs(:) /= 0)
216 vmean_dep(:) = vmean_dep(:)/nobs(:)
217 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
218 vstd_dep(:) = SQRT(MAX(0.0,vstd_dep(:)))
219 vmean_abs(:) = vmean_abs(:)/nobs(:)
220 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
221 vstd_abs(:) = SQRT(MAX(0.0,vstd_abs(:)))
222 end where
223
224 write (6,270) icount,nobs(1:tovs%nchan)
225 270 format (/1X,'NUMBER OF DATA Total/ACCEPTED'/1X,i10/1X,15I10)
226
227 write (6,288)
228 288 format (/1X,'FIRST PASS: MEANS AND STANDARD DEVIATIONS')
229
230 do j=1, tovs%nchan
231 jv = j
232 write (6,289) jv, nobs(j), vmean_abs(j), vstd_abs(j), vmean_dep(j), vstd_dep(j)
233 289 format (1X,I5,I10,4F15.2)
234 end do
235
236 !-----------------------------------------------------------------------------
237 ! 3. SECOND PASS THROUGH DATA, EXTRA Q.C. (sigma-elimination)
238 ! ------ ---- ------- ----- ----- ----
239
240 300 continue
241
242 vmn(:) = vmean_dep(:)
243 vstd(:)= vstd_dep(:)
244
245 rewind(unit=10)
246
247 vmean_dep(:) = 0.0 ! Clear matrices
248 vmean_abs(:) = 0.0 ! Clear matrices
249 nobs(:) = 0
250 nsel(:) = 0
251 nband(:) = 0
252 nbandl(:,:) = 0
253 nobsy(:) = 0
254
255 xbar(:,:) = 0.0
256 ybar(:) = 0.0
257 xcov(:,:,:) = 0.0
258 ycov(:) = 0.0
259 xycov(:,:) = 0.0
260
261 icount = 0
262 loop2:&
263 do
264 icount = icount + 1
265
266 call da_read_biasprep(tovs,10,ierr)
267 if (ierr == 0) then ! not end
268 continue
269 elseif (ierr == 1) then ! end
270 exit
271 else ! error
272 stop 'read error in da_airmass_bias'
273 end if
274
275 ! latitude band
276 iband = INT(tovs%lat/30.000001) + 3
277 if (lscan) then
278 JSCAN = tovs%scanpos
279 if (global) then
280 call GET_SCORR(JPCHAN,SCORR(1:JPCHAN),tovs%lat,vmnrlb,JSCAN)
281 else
282 SCORR(1:tovs%nchan) = vmnrl(1:tovs%nchan,JSCAN)
283 end if
284 ELSE
285 JSCAN = KSCAN/2
286 SCORR(1:JPCHAN) = 0.0
287 end if
288
289 vec_dep(1:JPCHAN) = tovs%omb(1:JPCHAN) - SCORR(1:JPCHAN)
290 vec_abs(1:JPCHAN) = tovs%tb(1:JPCHAN) - SCORR(1:JPCHAN)
291
292 ! 3.2 QC:
293 if (tovs%sensor_id == 3) then
294 call qc_amsua(tovs)
295 elseif(tovs%sensor_id == 4) then
296 call qc_amsub(tovs)
297 end if
298
299 ! 3.3 limb scan check
300 !---------------------
301 if (check_limb) then
302 if ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCAN-1))) then
303 nsel(2) = nsel(2) + 1
304 CYCLE loop2
305 end if
306 end if
307
308 ! 3.4 Reject data outside radiosonde mask
309 !------------------------------------------
310 if (check_mask) then
311 call MASK(tovs%lat,tovs%lon,LMASK)
312 if (.NOT. LMASK) then
313 nsel(8) = nsel(8) + 1
314 CYCLE loop2
315 end if
316 end if
317
318 ! 3.5 Reject outliers : facx*sigma, sigma calculated in first pass : loop1
319 !--------------------------------------------------------------------------
320 do j=1, tovs%nchan
321 if ( (abs(vec_dep(j)-vmn(j)) > (vstd(j)*FAC)) ) then
322 tovs%qc_flag(j) = -1
323 end if
324 end do
325
326 ! mean/std statistics for relative scan-bias corrected values
327 do j=1, tovs%nchan
328 if ( tovs%qc_flag(j) == 1 ) then
329 jv = j
330 nobs(j) = nobs(j) + 1
331 vmean_dep(j) = vmean_dep(j) + vec_dep(j)
332 vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
333 vmean_abs(j) = vmean_abs(j) + vec_abs(j)
334 vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
335 end if
336 end do
337
338 PRED(1:tovs%npred) = tovs%pred(1:tovs%npred)
339
340 ! compute regression variables mean/var/cov: y:departure; x:predictors
341 do j=1, tovs%nchan
342 if ( tovs%qc_flag(j) == 1 ) then
343 jv = j
344
345 ybar(j) = ybar(j) + vec_dep(j) ! mean of y
346 ycov(j) = ycov(j) + vec_dep(j)*vec_dep(j) ! variance of y
347
348 do i=1, tovs%npred ! Covariances for regression
349 xbar(j,i) = xbar(j,i) + pred(i) ! mean of x
350 xycov(j,i) = xycov(j,i) + vec_dep(j)*pred(i) ! cov of x and y
351 do ii=1, tovs%npred
352 xcov(j,i,ii) = xcov(j,i,ii) + pred(i)*pred(ii) ! cov of x and x
353 end do
354 end do
355
356 end if
357 end do
358
359 end do loop2
360
361 365 continue
362
363 ! Calculate means, standard deviations and covariances
364
365 where (nobs(:) /= 0)
366 vmean_dep(:) = vmean_dep(:)/nobs(:)
367 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
368 vstd_dep(:) = SQRT(MAX(0.0,vstd_dep(:)))
369 vmean_abs(:) = vmean_abs(:)/nobs(:)
370 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
371 vstd_abs(:) = SQRT(MAX(0.0,vstd_abs(:)))
372 end where
373
374 do j=1, tovs%nchan
375 if (nobs(j) /= 0) then
376 ybar(j) = ybar(j)/nobs(j)
377 ycov(j) = ycov(j)/nobs(j) - ybar(j)*ybar(j)
378 do i=1, tovs%npred
379 xbar(j,i) = xbar(j,i)/nobs(j)
380 xycov(j,i) = xycov(j,i)/nobs(j) - xbar(j,i)*ybar(j)
381 end do
382 do i=1, tovs%npred
383 xcov(j,i,1:tovs%npred) = xcov(j,i,1:tovs%npred)/nobs(j) - xbar(j,i)*xbar(j,1:tovs%npred)
384 end do
385 end if
386 end do
387
388 write (6,270) icount,nobs(1:tovs%nchan)
389
390 write (6,388)
391 388 format (/1X,'SECOND PASS: MEANS AND STANDARD DEVIATIONS')
392
393 do j=1, tovs%nchan
394 jv = j
395 write (6,289) jv, nobs(j), vmean_abs(j), vstd_abs(j), vmean_dep(j), vstd_dep(j)
396 end do
397
398 print *, ' '
399 print *, 'PREDICTOR MEANS AND STANDARD DEVIATIONS'
400 do j=1, tovs%nchan
401 jv = j
402 print *, ' '
403 print *, 'CHANNEL ', jv, ' NOBS = ', nobs(j)
404 do i=1, tovs%npred
405 write (6,390) i, xbar(j,i), SQRT(xcov(j,i,i))
406 end do
407 end do
408 390 format (1X,I5,4F15.2)
409
410 !----------------------------------------------------------------------------
411 ! 4. CALCULATE REGRESSION COEFFICIENTS.
412 ! --------- ---------- ------------
413
414 do j=1, tovs%nchan
415 jv = j
416 if ( nobs(j) >= 10 ) then
417 print *, 'REGRESSION : CHANNEL ', jv
418 call REGRESS_ONE(tovs%npred,xbar(j,1:tovs%npred),ybar(j), &
419 xcov(j,1:tovs%npred,1:tovs%npred), &
420 ycov(j),xycov(j,1:tovs%npred), &
421 tovs%npred,coef(j,1:tovs%npred),coef0(j),reserr(j),&
422 rescov(j,j),xvec(j,1:tovs%npred,1:tovs%npred))
423 else
424 print *, 'nobs < 10: ignoring REGRESSION : CHANNEL ', jv
425 end if
426 end do
427
428 print *, 'PREDICTOR EIGENVECTORS'
429 do j=1, tovs%nchan
430 jv = j
431 print *, ' '
432 print *, 'CHANNEL ', jv
433 do i=1, tovs%npred
434 print 888, xvec(j,1:tovs%npred,i)
435 end do
436 end do
437 888 format(1X,6F12.4)
438
439 xcoef0 = 0.0
440 xreser = 0.0
441 xcoef = 0.0
442
443 do JJ=1,tovs%nchan
444 J = JJ
445 xcoef0(J) = coef0(JJ)
446 xreser(J) = reserr(JJ)
447 do I=1, tovs%npred
448 xcoef(J,I) = coef(JJ,I)
449 end do
450 end do
451
452 ! 5. OUTPUT RESULTS.
453 ! ------ -------
454
455 coef_year = real(CDATE/1000000)
456 coef_month = real(MOD(CDATE,1000000)/10000)
457 coef_day = real(MOD(CDATE,10000)/100) ! Set date
458 coef_time = real(MOD(CDATE,100))
459
460 write (6,*)
461 write (6,502) &
462 coef_year, coef_month, coef_day,coef_time
463 502 format (/1X, ' YEAR',F5.0,' MONTH',F5.0,' DAY',F5.0,' TIME',F6.0,//&
464 1X,3X,' CH MEAN STDDEV RES.SD COEFFICIENTS')
465
466 do I=1, tovs%nchan
467 write (6,505) I,vmean_dep(I),vstd_dep(I),xreser(I),(xcoef(I,J),J=1,tovs%npred),xcoef0(I)
468 505 format (1X,I3,3F8.2,8F12.5)
469 end do
470
471 write (6,508) (I,I=1,tovs%nchan)
472 508 format (/1X,'RESIDUAL ERROR COVARIANCE'/1X,21I6)
473
474 do J=1, tovs%nchan
475 write (6,510) (rescov(I,J),I=1,tovs%nchan)
476 end do
477 510 format (1X,21F6.2)
478
479 write (6,511) (I,I=1,tovs%nchan)
480 511 format (/1X,'RESIDUAL ERROR CORRELATION'/1X,21I6)
481
482 do J=1, tovs%nchan
483 do I=1, tovs%nchan
484 rescov(I,J) = rescov(I,J)/(reserr(I)*reserr(J))
485 end do
486 write (6,510) (rescov(I,J),I=1,tovs%nchan)
487 end do
488
489 !----------------------------------------------------------------------------
490
491 deallocate(tovs%tb)
492 deallocate(tovs%omb)
493 deallocate(tovs%qc_flag)
494 deallocate(tovs%cloud_flag)
495 deallocate(tovs%pred)
496
497 close(unit=10)
498 if (lscan) close(unit=12)
499
500 ! out coefs to ASCII file bcor.asc
501 call write_biascoef(tovs%nchan,kscan,jband,tovs%npred,global, &
502 VMNRL(1:tovs%nchan,1:kscan),VMNRLB(1:tovs%nchan,1:kscan,1:jband), &
503 xcoef(1:tovs%nchan,1:tovs%npred),xcoef0(1:tovs%nchan), &
504 nobs(1:tovs%nchan),vmean_abs(1:tovs%nchan), &
505 vstd_abs(1:tovs%nchan),vmean_dep(1:tovs%nchan), &
506 vstd_dep(1:tovs%nchan) )
507
508 end program da_bias_airmass