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 ENDIF
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 ENDDO
125 ENDDO
126 else
127 DO J=1, tovs%nchan ! get scan biases (no band dependence)
128 READ(12) JJ, vmnrl(J,1:KSCAN)
129 ENDDO
130 end if
131 ENDIF
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 endif
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 ENDIF
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 elseif(tovs%sensor_id == 15) then ! mhs
172 call qc_amsub(tovs)
173 endif
174
175 !-------------------------
176 ! 2.3 limb sounding check
177 !-------------------------
178 if (check_limb) then
179 IF ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCAN-1))) THEN ! Reject edge of scan
180 nsel(2) = nsel(2) + 1
181 CYCLE loop1
182 ENDIF
183 end if
184
185 !-----------------------------------------
186 ! 2.4 Reject data outside radiosonde mask
187 !-----------------------------------------
188 IF (check_mask) THEN
189 CALL MASK(tovs%lat,tovs%lon,LMASK)
190 IF (.NOT. LMASK) THEN
191 nsel(8) = nsel(8) + 1
192 CYCLE loop1
193 ENDIF
194 ENDIF
195
196 ! Good data : count and use in the statistics
197
198 DO j=1, tovs%nchan
199 IF ( tovs%qc_flag(j) == 1 ) THEN
200 ! statistics for channel
201 !--------------------------
202 nobs(j) = nobs(j) + 1
203 vmean_dep(j) = vmean_dep(j) + vec_dep(j)
204 vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
205 vmean_abs(j) = vmean_abs(j) + vec_abs(j)
206 vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
207 ENDIF
208 ENDDO
209
210 ENDDO loop1
211
212 265 CONTINUE
213
214 !---------------------------------
215 ! 2.8 mean and std for channels
216 !---------------------------------
217 WHERE (nobs(:) /= 0)
218 vmean_dep(:) = vmean_dep(:)/nobs(:)
219 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
220 vstd_dep(:) = SQRT(MAX(0.0,vstd_dep(:)))
221 vmean_abs(:) = vmean_abs(:)/nobs(:)
222 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
223 vstd_abs(:) = SQRT(MAX(0.0,vstd_abs(:)))
224 END WHERE
225
226 WRITE (6,270) icount,nobs(1:tovs%nchan)
227 270 FORMAT (/1X,'NUMBER OF DATA Total/ACCEPTED'/1X,i10/1X,15I10)
228
229 WRITE (6,288)
230 288 FORMAT (/1X,'FIRST PASS: MEANS AND STANDARD DEVIATIONS')
231
232 DO j=1, tovs%nchan
233 jv = j
234 WRITE (6,289) jv, nobs(j), vmean_abs(j), vstd_abs(j), vmean_dep(j), vstd_dep(j)
235 289 FORMAT (1X,I5,I10,4F15.2)
236 ENDDO
237
238 !-----------------------------------------------------------------------------
239 ! 3. SECOND PASS THROUGH DATA, EXTRA Q.C. (sigma-elimination)
240 ! ------ ---- ------- ----- ----- ----
241
242 300 CONTINUE
243
244 vmn(:) = vmean_dep(:)
245 vstd(:)= vstd_dep(:)
246
247 REWIND(UNIT=10)
248
249 vmean_dep(:) = 0.0 ! Clear matrices
250 vmean_abs(:) = 0.0 ! Clear matrices
251 nobs(:) = 0
252 nsel(:) = 0
253 nband(:) = 0
254 nbandl(:,:) = 0
255 nobsy(:) = 0
256
257 xbar(:,:) = 0.0
258 ybar(:) = 0.0
259 xcov(:,:,:) = 0.0
260 ycov(:) = 0.0
261 xycov(:,:) = 0.0
262
263 icount = 0
264 loop2:&
265 DO
266 icount = icount + 1
267
268 call da_read_biasprep(tovs,10,ierr)
269 if (ierr == 0) then ! not end
270 continue
271 elseif (ierr == 1) then ! end
272 exit
273 else ! error
274 stop 'read error in da_airmass_bias'
275 endif
276
277 ! latitude band
278 iband = INT(tovs%lat/30.000001) + 3
279 IF (lscan) THEN
280 JSCAN = tovs%scanpos
281 if (global) then
282 CALL GET_SCORR(JPCHAN,SCORR(1:JPCHAN),tovs%lat,vmnrlb,JSCAN)
283 else
284 SCORR(1:tovs%nchan) = vmnrl(1:tovs%nchan,JSCAN)
285 end if
286 ELSE
287 JSCAN = KSCAN/2
288 SCORR(1:JPCHAN) = 0.0
289 ENDIF
290
291 vec_dep(1:JPCHAN) = tovs%omb(1:JPCHAN) - SCORR(1:JPCHAN)
292 vec_abs(1:JPCHAN) = tovs%tb(1:JPCHAN) - SCORR(1:JPCHAN)
293
294 ! 3.2 QC:
295 if (tovs%sensor_id == 3) then
296 call qc_amsua(tovs)
297 elseif(tovs%sensor_id == 4) then
298 call qc_amsub(tovs)
299 elseif(tovs%sensor_id == 15) then
300 call qc_amsub(tovs)
301 endif
302
303 ! 3.3 limb scan check
304 !---------------------
305 if (check_limb) then
306 IF ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCAN-1))) THEN
307 nsel(2) = nsel(2) + 1
308 CYCLE loop2
309 ENDIF
310 end if
311
312 ! 3.4 Reject data outside radiosonde mask
313 !------------------------------------------
314 IF (check_mask) THEN
315 CALL MASK(tovs%lat,tovs%lon,LMASK)
316 IF (.NOT. LMASK) THEN
317 nsel(8) = nsel(8) + 1
318 CYCLE loop2
319 ENDIF
320 ENDIF
321
322 ! 3.5 Reject outliers : facx*sigma, sigma calculated in first pass : loop1
323 !--------------------------------------------------------------------------
324 DO j=1, tovs%nchan
325 IF ( (ABS(vec_dep(j)-vmn(j)) > (vstd(j)*FAC)) ) THEN
326 tovs%qc_flag(j) = -1
327 ENDIF
328 ENDDO
329
330 ! mean/std statistics for relative scan-bias corrected values
331 DO j=1, tovs%nchan
332 IF ( tovs%qc_flag(j) == 1 ) THEN
333 jv = j
334 nobs(j) = nobs(j) + 1
335 vmean_dep(j) = vmean_dep(j) + vec_dep(j)
336 vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
337 vmean_abs(j) = vmean_abs(j) + vec_abs(j)
338 vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
339 ENDIF
340 ENDDO
341
342 PRED(1:tovs%npred) = tovs%pred(1:tovs%npred)
343
344 ! compute regression variables mean/var/cov: y:departure; x:predictors
345 DO j=1, tovs%nchan
346 IF ( tovs%qc_flag(j) == 1 ) THEN
347 jv = j
348
349 ybar(j) = ybar(j) + vec_dep(j) ! mean of y
350 ycov(j) = ycov(j) + vec_dep(j)*vec_dep(j) ! variance of y
351
352 DO i=1, tovs%npred ! Covariances for regression
353 xbar(j,i) = xbar(j,i) + pred(i) ! mean of x
354 xycov(j,i) = xycov(j,i) + vec_dep(j)*pred(i) ! cov of x and y
355 DO ii=1, tovs%npred
356 xcov(j,i,ii) = xcov(j,i,ii) + pred(i)*pred(ii) ! cov of x and x
357 ENDDO
358 ENDDO
359
360 ENDIF
361 ENDDO
362
363 ENDDO loop2
364
365 365 CONTINUE
366
367 ! Calculate means, standard deviations and covariances
368
369 WHERE (nobs(:) /= 0)
370 vmean_dep(:) = vmean_dep(:)/nobs(:)
371 vstd_dep(:) = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
372 vstd_dep(:) = SQRT(MAX(0.0,vstd_dep(:)))
373 vmean_abs(:) = vmean_abs(:)/nobs(:)
374 vstd_abs(:) = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
375 vstd_abs(:) = SQRT(MAX(0.0,vstd_abs(:)))
376 ENDWHERE
377
378 DO j=1, tovs%nchan
379 IF (nobs(j) /= 0) THEN
380 ybar(j) = ybar(j)/nobs(j)
381 ycov(j) = ycov(j)/nobs(j) - ybar(j)*ybar(j)
382 DO i=1, tovs%npred
383 xbar(j,i) = xbar(j,i)/nobs(j)
384 xycov(j,i) = xycov(j,i)/nobs(j) - xbar(j,i)*ybar(j)
385 ENDDO
386 DO i=1, tovs%npred
387 xcov(j,i,1:tovs%npred) = xcov(j,i,1:tovs%npred)/nobs(j) - xbar(j,i)*xbar(j,1:tovs%npred)
388 ENDDO
389 ENDIF
390 ENDDO
391
392 WRITE (6,270) icount,nobs(1:tovs%nchan)
393
394 WRITE (6,388)
395 388 FORMAT (/1X,'SECOND PASS: MEANS AND STANDARD DEVIATIONS')
396
397 DO j=1, tovs%nchan
398 jv = j
399 WRITE (6,289) jv, nobs(j), vmean_abs(j), vstd_abs(j), vmean_dep(j), vstd_dep(j)
400 ENDDO
401
402 PRINT *, ' '
403 PRINT *, 'PREDICTOR MEANS AND STANDARD DEVIATIONS'
404 DO j=1, tovs%nchan
405 jv = j
406 PRINT *, ' '
407 PRINT *, 'CHANNEL ', jv, ' NOBS = ', nobs(j)
408 DO i=1, tovs%npred
409 WRITE (6,390) i, xbar(j,i), SQRT(xcov(j,i,i))
410 ENDDO
411 ENDDO
412 390 FORMAT (1X,I5,4F15.2)
413
414 !----------------------------------------------------------------------------
415 ! 4. CALCULATE REGRESSION COEFFICIENTS.
416 ! --------- ---------- ------------
417
418 DO j=1, tovs%nchan
419 jv = j
420 if ( nobs(j) >= 10 ) then
421 PRINT *, 'REGRESSION : CHANNEL ', jv
422 CALL REGRESS_ONE(tovs%npred,xbar(j,1:tovs%npred),ybar(j), &
423 xcov(j,1:tovs%npred,1:tovs%npred), &
424 ycov(j),xycov(j,1:tovs%npred), &
425 tovs%npred,coef(j,1:tovs%npred),coef0(j),reserr(j),&
426 rescov(j,j),xvec(j,1:tovs%npred,1:tovs%npred))
427 else
428 PRINT *, 'nobs < 10: ignoring REGRESSION : CHANNEL ', jv
429 end if
430 ENDDO
431
432 PRINT *, 'PREDICTOR EIGENVECTORS'
433 DO j=1, tovs%nchan
434 jv = j
435 PRINT *, ' '
436 PRINT *, 'CHANNEL ', jv
437 DO i=1, tovs%npred
438 PRINT 888, xvec(j,1:tovs%npred,i)
439 ENDDO
440 ENDDO
441 888 FORMAT(1X,6F12.4)
442
443 xcoef0 = 0.0
444 xreser = 0.0
445 xcoef = 0.0
446
447 DO JJ=1,tovs%nchan
448 J = JJ
449 xcoef0(J) = coef0(JJ)
450 xreser(J) = reserr(JJ)
451 DO I=1, tovs%npred
452 xcoef(J,I) = coef(JJ,I)
453 ENDDO
454 ENDDO
455
456 ! 5. OUTPUT RESULTS.
457 ! ------ -------
458
459 coef_year = REAL(CDATE/1000000)
460 coef_month = REAL(MOD(CDATE,1000000)/10000)
461 coef_day = REAL(MOD(CDATE,10000)/100) ! Set date
462 coef_time = REAL(MOD(CDATE,100))
463
464 WRITE (6,*)
465 WRITE (6,502) &
466 coef_year, coef_month, coef_day,coef_time
467 502 FORMAT (/1X, ' YEAR',F5.0,' MONTH',F5.0,' DAY',F5.0,' TIME',F6.0,//&
468 1X,3X,' CH MEAN STDDEV RES.SD COEFFICIENTS')
469
470 DO I=1, tovs%nchan
471 WRITE (6,505) I,vmean_dep(I),vstd_dep(I),xreser(I),(xcoef(I,J),J=1,tovs%npred),xcoef0(I)
472 505 FORMAT (1X,I3,3F8.2,8F12.5)
473 ENDDO
474
475 WRITE (6,508) (I,I=1,tovs%nchan)
476 508 FORMAT (/1X,'RESIDUAL ERROR COVARIANCE'/1X,21I6)
477
478 DO J=1, tovs%nchan
479 WRITE (6,510) (rescov(I,J),I=1,tovs%nchan)
480 ENDDO
481 510 FORMAT (1X,21F6.2)
482
483 WRITE (6,511) (I,I=1,tovs%nchan)
484 511 FORMAT (/1X,'RESIDUAL ERROR CORRELATION'/1X,21I6)
485
486 DO J=1, tovs%nchan
487 DO I=1, tovs%nchan
488 rescov(I,J) = rescov(I,J)/(reserr(I)*reserr(J))
489 ENDDO
490 WRITE (6,510) (rescov(I,J),I=1,tovs%nchan)
491 ENDDO
492
493 !----------------------------------------------------------------------------
494
495 deallocate(tovs%tb)
496 deallocate(tovs%omb)
497 deallocate(tovs%qc_flag)
498 deallocate(tovs%cloud_flag)
499 deallocate(tovs%pred)
500
501 CLOSE(UNIT=10)
502 IF (lscan) CLOSE(UNIT=12)
503
504 ! out coefs to ASCII file bcor.asc
505 call write_biascoef(tovs%nchan,kscan,jband,tovs%npred,global, &
506 VMNRL(1:tovs%nchan,1:kscan),VMNRLB(1:tovs%nchan,1:kscan,1:jband), &
507 xcoef(1:tovs%nchan,1:tovs%npred),xcoef0(1:tovs%nchan), &
508 nobs(1:tovs%nchan),vmean_abs(1:tovs%nchan), &
509 vstd_abs(1:tovs%nchan),vmean_dep(1:tovs%nchan), &
510 vstd_dep(1:tovs%nchan) )
511
512 END PROGRAM da_bias_airmass