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 endif
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 ENDIF
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 ENDIF
192 ENDIF
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 ENDIF
206 ENDDO
207
208 ENDDO 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 ENDDO
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 endif
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 ENDIF
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 endif
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 ENDIF
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 ENDIF
316 ENDIF
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 ENDIF
324 ENDDO
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 ENDIF
336 ENDDO
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 ENDDO
354 ENDDO
355
356 ENDIF
357 ENDDO
358
359 ENDDO 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 ENDWHERE
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 ENDDO
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 ENDDO
385 ENDIF
386 ENDDO
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 ENDDO
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 ENDDO
407 ENDDO
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 ENDDO
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 ENDDO
436 ENDDO
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 ENDDO
450 ENDDO
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 ENDDO
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 ENDDO
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 ENDDO
486 WRITE (6,510) (rescov(I,J),I=1,tovs%nchan)
487 ENDDO
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