rad_bias.f90
References to this file elsewhere.
1 MODULE rad_bias
2
3 ! PRINT_BIAS, MASK
4
5 IMPLICIT NONE
6
7 ! INTEGER, PARAMETER :: LONG = SELECTED_REAL_KIND(9,99)
8 INTEGER, PARAMETER :: LONG = SELECTED_REAL_KIND(15,37)
9
10 INTEGER, PARAMETER :: JPNX=8
11 INTEGER, PARAMETER :: JPXEIG=10
12
13 INTEGER, PARAMETER :: JPNY=20
14 INTEGER, PARAMETER :: JPCHAN=20
15 INTEGER, PARAMETER :: JPSCAN=90
16
17 INTEGER, PARAMETER :: JBAND=18
18 INTEGER, PARAMETER :: BOFF = JBAND/2 + 1
19 REAL, PARAMETER :: BDIV = 180.0/JBAND + 0.0001
20 INTEGER, PARAMETER :: AVBAND=2
21
22 REAL(KIND=LONG), PARAMETER :: PI = 3.141592654
23 REAL(KIND=LONG), PARAMETER :: VMAX = 350.0, VMIN = 50.0, VDMAX = 20.0
24 CHARACTER(LEN=10), PARAMETER :: labsel(10) = &
25 (/' WRONG SAT',' PATH 2+3',' LAND', &
26 ' ROGUE TB',' ROGUE DTB','BAD Win CH', &
27 ' OUTLIERS',' MASKED',' ', &
28 ' GOOD'/)
29
30 TYPE BIAS
31 INTEGER :: nchan ! number of channels
32 INTEGER :: npred ! number of predictors
33 INTEGER :: platform_id,satellite_id,sensor_id
34 INTEGER :: year, month, day, hour, min, sec
35 INTEGER :: scanline,scanpos
36 INTEGER :: landmask
37 INTEGER :: surf_flag
38 INTEGER, pointer :: qc_flag(:) ! 1/0:good/bad
39 INTEGER, pointer :: cloud_flag(:) ! 1/0:no-cloud/cloud
40 REAL :: elevation,lat,lon,ps, t2m, q2m, tsk, clwp
41 REAL, pointer :: tb(:), omb(:), bias(:)
42 REAL, pointer :: pred(:)
43 END TYPE BIAS
44
45 CONTAINS
46
47 SUBROUTINE PRINT_BIAS(rads)
48
49 IMPLICIT NONE
50
51 TYPE(BIAS), INTENT(IN) :: rads
52
53 INTEGER :: lev
54
55 WRITE(6,*) ' '
56
57 WRITE(6,'(a20,3i4)') 'Instrument Triplet :', rads%platform_id, &
58 rads%satellite_id,rads%sensor_id
59 WRITE(6,'(a20,6i5)') 'DATE :', rads%year, &
60 rads%month,rads%day,rads%hour,rads%min,rads%sec
61
62 WRITE(6,'(a20,2f10.3)') 'Latitude Longitude :', rads%lat,rads%lon
63
64 WRITE(6,'(a20,2i5)') 'Scan Line/Position :', rads%scanline,rads%scanpos
65 WRITE(6,'(a20,i3,f10.2)') 'Land/Sea Elevation :', rads%landmask, rads%elevation
66
67 WRITE(6,'(a20,4f10.2)') 'Ps Ts T2m Q2m :', rads%ps, rads%tsk, rads%t2m, rads%q2m
68
69 WRITE(6,'(a20,i5)') 'Number of Channels :', rads%nchan
70 DO lev=1, rads%nchan
71 WRITE(6,'(10x,i3,2f10.2)') lev, rads%tb(lev), rads%omb(lev)
72 ENDDO
73
74 WRITE(6,'(a60)') ' predictors : 1000-300/200-50mb thickness T_skin TPW'
75 WRITE(6,'(4f10.2)') rads%pred(1:rads%npred)
76
77 RETURN
78 END SUBROUTINE PRINT_BIAS
79
80 SUBROUTINE MASK(PLAT,PLONG,LMASK)
81
82 ! ROUTINE TO ACCESS RADIOSONDE MASK:
83 ! RETURNS LMASK=.TRUE. IF PLAT,PLONG CLOSE TO RADIOSONDE
84 ! LOCATION. DISTANCE DEFINED IN CREATION OF MASK.
85
86 IMPLICIT NONE
87
88 REAL, INTENT(INOUT) :: PLAT, PLONG
89 LOGICAL, INTENT(OUT) :: LMASK
90
91 INTEGER :: IMASK(360,181) = 0
92 INTEGER :: IERR, ICL, ILEN
93 ! INTEGER :: CRAYOPEN, CRAYREAD, CRAYCLOSE
94 INTEGER :: IUMASK = 39, LOOP = 0
95
96 INTEGER, SAVE :: INIT = 0
97 INTEGER :: IX, IY
98
99 CHARACTER(LEN=1) :: C(600000)
100
101 IF (INIT == 0) THEN
102
103 INIT = 1
104 open(169,FILE='mask_asc')
105 print * ,' open radiosonde mask file mask_asc '
106 read(169,99)imask
107 99 FORMAT (1X,90I1)
108 close(169)
109
110 ENDIF
111
112 LOOP = LOOP + 1
113 IF (PLONG < 0.0) PLONG = PLONG + 360.0
114 IX = MOD(NINT(PLONG),360) + 1
115 IY = NINT(90.0-PLAT) + 1
116 IF (IMASK(IX,IY) == 1) THEN
117 LMASK = .TRUE.
118 ELSE
119 LMASK = .FALSE.
120 ENDIF
121
122 ! IF (LOOP.LE.100) WRITE(6,*) LOOP,PLAT,PLONG,IX,IY,IMASK(IX,IY),LMASK
123
124 RETURN
125 END SUBROUTINE MASK
126
127 LOGICAL FUNCTION USE_CHAN(j,land,path,ICLR,ICLD,LNDCH)
128
129 IMPLICIT NONE
130
131 INTEGER, INTENT(IN) :: j, land, path
132 INTEGER, INTENT(IN) :: ICLR(JPNY), ICLD(JPNY), LNDCH(JPNY)
133
134 INTEGER :: Lmask, Cmask
135
136 SELECT CASE(land)
137 CASE(0) ! over land
138 Lmask = LNDCH(J) ! 0:not use; 1: use
139 CASE(1) ! over sea
140 Lmask = 1
141 CASE DEFAULT
142 Lmask = 0
143 END SELECT
144
145 SELECT CASE(path)
146 CASE(0) ! clear sky
147 Cmask = ICLR(J) ! 0:not use; 1: use
148 CASE(1) ! cloudy sky
149 Cmask = ICLD(J) ! 0:not use; 1: use
150 CASE DEFAULT
151 Cmask = 0
152 END SELECT
153
154 IF (Lmask*Cmask == 0) THEN
155 USE_CHAN = .FALSE.
156 ELSE
157 USE_CHAN = .TRUE.
158 ENDIF
159
160 RETURN
161 END FUNCTION USE_CHAN
162
163 SUBROUTINE GET_SCORR(JPCHAN,SCORR,LAT,vmnrlb,JSCAN)
164
165 IMPLICIT NONE
166
167 ! Assume JBAND, JPCHAN, JPSCAN, BDIV and BOFF are inherited from main program.
168 ! Will need to be explicitly defined for 1DVAR, presumably #include in future.
169
170 INTEGER, INTENT(IN) :: JPCHAN
171 REAL(KIND=LONG), INTENT(OUT) :: SCORR(JPCHAN)
172 REAL, INTENT(INOUT) :: LAT
173 REAL(KIND=LONG), INTENT(INOUT) :: vmnrlb(JPCHAN,JPSCAN,JBAND)
174 INTEGER, INTENT(IN) :: JSCAN
175
176 INTEGER :: sband,i
177 REAL :: BLAT
178
179 sband = FLOOR(LAT/BDIV) + BOFF
180 BLAT = FLOOR(LAT/BDIV)*BDIV
181
182 IF (LAT >= BLAT+BDIV/2) THEN
183
184 IF (sband < JBAND) THEN
185 SCORR(1:JPCHAN) = (LAT - (BLAT+BDIV/2)) * vmnrlb(1:JPCHAN,JSCAN,sband+1) / BDIV &
186 + ((BLAT+3*BDIV/2) - LAT) * vmnrlb(1:JPCHAN,JSCAN,sband ) / BDIV
187 ELSE
188 SCORR(1:JPCHAN) = vmnrlb(1:JPCHAN,JSCAN,sband)
189 ENDIF
190
191 ELSEIF (LAT < BLAT+BDIV/2) THEN
192
193 IF (sband > 1) THEN
194 SCORR(1:JPCHAN) = (LAT - (BLAT-BDIV/2)) * vmnrlb(1:JPCHAN,JSCAN,sband ) / BDIV &
195 + ((BLAT+BDIV/2) - LAT) * vmnrlb(1:JPCHAN,JSCAN,sband-1) / BDIV
196 ELSE
197 SCORR(1:JPCHAN) = vmnrlb(1:JPCHAN,JSCAN,sband)
198 ENDIF
199
200 ENDIF
201
202 END SUBROUTINE GET_SCORR
203
204 SUBROUTINE QC_AMSUA(tovs)
205
206 TYPE(bias), intent(inout) :: tovs
207
208 integer :: j, jv
209
210 tovs%qc_flag(:) = 1
211 !--------------------------------------------------------
212 ! 2.1 extrem values
213 !--------------------------------------------------------
214 DO j=1, tovs%nchan ! Reject silly values
215 jv = j
216 IF ((tovs%tb(jv) > VMAX) .OR. (tovs%tb(jv) < VMIN) ) then
217 tovs%qc_flag(jv) = -1
218 ENDIF
219 ENDDO
220
221 !------------------------------------
222 ! 2.2 departure extrem values test
223 !------------------------------------
224 DO j=1, tovs%nchan
225 jv = j
226 IF (ABS(tovs%omb(jv)) > VDMAX) THEN
227 tovs%qc_flag(jv) = -1
228 ENDIF
229 ENDDO
230
231 !------------------------------------
232 ! 2.3 window channel
233 !------------------------------------
234 IF (tovs%surf_flag > 0) THEN ! not over sea
235 tovs%qc_flag(1:3) = -1
236 tovs%qc_flag(15) = -1
237 ENDIF
238
239 END SUBROUTINE QC_AMSUA
240
241 SUBROUTINE QC_AMSUB(tovs)
242
243 TYPE(bias), intent(inout) :: tovs
244 integer :: j, jv
245
246 tovs%qc_flag(:) = 1
247 !--------------------------------------------------------
248 ! 2.1 extrem values
249 !--------------------------------------------------------
250 DO j=1, tovs%nchan ! Reject silly values
251 jv = j
252 IF ((tovs%tb(jv) > VMAX) .OR. (tovs%tb(jv) < VMIN) ) then
253 tovs%qc_flag(jv) = -1
254 ENDIF
255 ENDDO
256
257 !------------------------------------
258 ! 2.2 departure extrem values test
259 !------------------------------------
260 DO j=1, tovs%nchan
261 jv = j
262 IF (ABS(tovs%omb(jv)) > VDMAX) THEN
263 tovs%qc_flag(jv) = -1
264 ENDIF
265 ENDDO
266
267 !------------------------------------
268 ! 2.3 window channel
269 !------------------------------------
270 IF (tovs%surf_flag > 0) THEN ! not over sea
271 tovs%qc_flag(1:2) = -1
272 ENDIF
273
274 END SUBROUTINE QC_AMSUB
275
276 subroutine da_read_biasprep(radbias,biasprep_unit,ierr)
277 TYPE(bias), INTENT(INOUT) :: radbias
278 integer, intent(in) :: biasprep_unit
279 integer, intent(out) :: ierr
280 read(UNIT=biasprep_unit,END=990,ERR=995) radbias%nchan,radbias%npred
281 read(UNIT=biasprep_unit,END=990,ERR=995) radbias%platform_id , &
282 radbias%satellite_id, &
283 radbias%sensor_id, &
284 radbias%year,radbias%month,&
285 radbias%day, radbias%hour, &
286 radbias%min, radbias%sec, &
287 radbias%scanline, &
288 radbias%scanpos, &
289 radbias%landmask, &
290 radbias%elevation,&
291 radbias%lat,radbias%lon, &
292 radbias%ps,radbias%t2m, &
293 radbias%q2m,radbias%tsk, &
294 radbias%tb(1:radbias%nchan), &
295 radbias%omb(1:radbias%nchan), &
296 radbias%bias(1:radbias%nchan), &
297 radbias%pred(1:radbias%npred), &
298 radbias%qc_flag(1:radbias%nchan), &
299 radbias%cloud_flag(1:radbias%nchan), &
300 radbias%surf_flag, radbias%clwp
301
302 ierr = 0
303 RETURN
304
305 990 CONTINUE
306 ierr = 1
307 RETURN
308
309 995 CONTINUE
310 ierr = 2
311 RETURN
312
313 end subroutine da_read_biasprep
314
315 subroutine da_write_biasprep(radbias,biasprep_unit)
316 TYPE(bias), INTENT(IN) :: radbias
317 integer, INTENT(IN) :: biasprep_unit
318
319 write(UNIT=biasprep_unit) radbias%nchan,radbias%npred
320 write(UNIT=biasprep_unit) radbias%platform_id , &
321 radbias%satellite_id, &
322 radbias%sensor_id, &
323 radbias%year,radbias%month,&
324 radbias%day, radbias%hour, &
325 radbias%min, radbias%sec, &
326 radbias%scanline, &
327 radbias%scanpos, &
328 radbias%landmask, &
329 radbias%elevation,&
330 radbias%lat,radbias%lon, &
331 radbias%ps,radbias%t2m, &
332 radbias%q2m,radbias%tsk, &
333 radbias%tb(1:radbias%nchan), &
334 radbias%omb(1:radbias%nchan), &
335 radbias%bias(1:radbias%nchan), &
336 radbias%pred(1:radbias%npred), &
337 radbias%qc_flag(1:radbias%nchan), &
338 radbias%cloud_flag(1:radbias%nchan), &
339 radbias%surf_flag, radbias%clwp
340 end subroutine da_write_biasprep
341
342 subroutine write_biascoef(nchan,nscan,nband,npred,global, &
343 scanbias,scanbias_b,coef,coef0, &
344 nobs,vmean_abs,vstd_abs,vmean_dep,vstd_dep)
345
346 integer, intent(in) :: nchan,nscan,nband,npred,nobs(nchan)
347 logical, intent(in) :: global
348 real(KIND=LONG), intent(in) :: scanbias(nchan,nscan), &
349 scanbias_b(nchan,nscan,nband), &
350 coef(nchan,npred),coef0(nchan), &
351 vmean_abs(nchan), vstd_abs(nchan), &
352 vmean_dep(nchan),vstd_dep(nchan)
353 integer :: iunit,j,i,stdout
354
355 iunit=99
356 stdout=6
357 open(UNIT=iunit,file='bcor.asc',form='formatted')
358
359 write (iunit,'(4i6)') nchan,nscan,nband,npred
360 do i=1, nchan
361 write (iunit,'(i5,i7,4F8.2)') i,nobs(i),vmean_abs(i),vstd_abs(i),vmean_dep(i),vstd_dep(i)
362 end do
363
364 do i=1, nchan
365 write (iunit,'(i5,5F12.5)') i,(coef(i,j),j=1,npred),coef0(i)
366 end do
367
368 if (global) then
369 if (nscan == 30) then ! amsua
370 write (iunit,'(a,/8X,30I7)') 'RELATIVE SCAN BIASES ',(j,j=1,nscan)
371 else if (nscan == 90) then ! amsub
372 write (iunit,'(a,/8X,90I7)') 'RELATIVE SCAN BIASES ',(j,j=1,nscan)
373 else if (nscan == 56) then ! hirs
374 write (iunit,'(a,/8X,56I7)') 'RELATIVE SCAN BIASES ',(j,j=1,nscan)
375 else
376 write(UNIT=stdout,FMT=*) 'Unknown Sensor Type'
377 end if
378 do j=1, nchan
379 do i=1, nband
380 if (nscan == 30) then ! amsua
381 write(iunit,'(i5,i3,30F7.2)') j,i, scanbias_b(j,1:nscan,i)
382 else if (nscan == 90) then ! amsub
383 write(iunit,'(i5,i3,90F7.2)') j,i, scanbias_b(j,1:nscan,i)
384 else if (nscan == 56) then ! hirs
385 write(iunit,'(i5,i3,56F7.2)') j,i, scanbias_b(j,1:nscan,i)
386 else
387 write(UNIT=stdout,FMT=*) 'Unknown Sensor Type'
388 end if
389 end do
390 end do
391 else
392 if (nscan == 30) then ! amsua
393 write (iunit,'(a,/5X,30I7)') 'RELATIVE SCAN BIASES ',(j,j=1,nscan)
394 else if (nscan == 90) then ! amsub
395 write (iunit,'(a,/5X,90I7)') 'RELATIVE SCAN BIASES ',(j,j=1,nscan)
396 else if (nscan == 56) then ! hirs
397 write (iunit,'(a,/5X,56I7)') 'RELATIVE SCAN BIASES ',(j,j=1,nscan)
398 else
399 write(UNIT=stdout,FMT=*) 'Unknown Sensor Type'
400 end if
401 do j=1, nchan
402 if (nscan == 30) then ! amsua
403 write(iunit,'(i5,30F7.2)') j, scanbias(j,1:nscan)
404 else if (nscan == 90) then ! amsub
405 write(iunit,'(i5,90F7.2)') j, scanbias(j,1:nscan)
406 else if (nscan == 56) then ! hirs
407 write(iunit,'(i5,56F7.2)') j, scanbias(j,1:nscan)
408 else
409 write(UNIT=stdout,FMT=*) 'Unknown Sensor Type'
410 end if
411 end do
412 end if
413
414 close(iunit)
415 end subroutine write_biascoef
416
417 subroutine read_biascoef(nchan,nscan,nband,npred,global, &
418 scanbias,scanbias_b,coef,coef0, &
419 nobs,vmean_abs,vstd_abs,vmean_dep,vstd_dep)
420
421 integer, intent(in) :: nchan,nscan,nband,npred
422 logical, intent(in) :: global
423 integer, intent(out) :: nobs(nchan)
424 real(KIND=LONG), intent(out) :: scanbias(nchan,nscan), &
425 scanbias_b(nchan,nscan,nband), &
426 coef(nchan,npred),coef0(nchan), &
427 vmean_abs(nchan), vstd_abs(nchan), &
428 vmean_dep(nchan),vstd_dep(nchan)
429
430 integer :: iunit,j,i,stdout, ii,jj
431
432 iunit=99
433 stdout=6
434 open(UNIT=iunit,file='bcor.asc',form='formatted')
435
436 !read (iunit,'(4i6)') nchan,nscan,nband,npred
437 read (iunit,'(4i6)')
438 do i=1, nchan
439 read (iunit,'(i5,i7,4F8.2)') ii,nobs(i),vmean_abs(i),vstd_abs(i),vmean_dep(i),vstd_dep(i)
440 end do
441
442 do i=1, nchan
443 read (iunit,'(i5,5F12.5)') ii,(coef(i,j),j=1,npred),coef0(i)
444 end do
445
446 read (iunit,*)
447 read (iunit,*)
448 if (global) then
449 do j=1, nchan
450 do i=1, nband
451 if (nscan == 30) then ! amsua
452 read(iunit,'(i5,i3,30F7.2)') jj,ii, scanbias_b(j,1:nscan,i)
453 else if (nscan == 90) then ! amsub
454 read(iunit,'(i5,i3,90F7.2)') jj,ii, scanbias_b(j,1:nscan,i)
455 else if (nscan == 56) then ! hirs
456 read(iunit,'(i5,i3,56F7.2)') jj,ii, scanbias_b(j,1:nscan,i)
457 else
458 write(UNIT=stdout,FMT=*) 'Unknown Sensor Type'
459 end if
460 end do
461 end do
462 else
463 do j=1, nchan
464 if (nscan == 30) then ! amsua
465 read(iunit,'(i5,30F7.2)') jj, scanbias(j,1:nscan)
466 else if (nscan == 90) then ! amsub
467 read(iunit,'(i5,90F7.2)') jj, scanbias(j,1:nscan)
468 else if (nscan == 56) then ! hirs
469 read(iunit,'(i5,56F7.2)') jj, scanbias(j,1:nscan)
470 else
471 write(UNIT=stdout,FMT=*) 'Unknown Sensor Type'
472 end if
473 end do
474 end if
475
476 close(iunit)
477 end subroutine read_biascoef
478
479 END MODULE rad_bias