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