1module bgckSSMIS_mod
2 ! MODULE bgckSSMIS_mod (prefix='ssbg' category='1. High-level functionality')
3 !
4 !:Purpose: To perform background check and quality control for SSMIS radiance observations.
5 !
6 use midasMpi_mod
7 use MathPhysConstants_mod
8 use utilities_mod
9 use obsSpaceData_mod
10 use tovsNL_mod
11 use obsErrors_mod
12
13 implicit none
14 save
15 private
16
17 ! Public functions/subroutines
18 public :: ssbg_computeSsmisSurfaceType
19 public :: ssbg_bgCheckSSMIS
20 real :: ssbg_clwQcThreshold
21 logical :: ssbg_debug
22
23 real, parameter :: ssbg_realMissing=-99.
24 integer, parameter :: ssbg_intMissing=-1
25 ! Other variables:
26 real, parameter :: ssbg_rmisg=-999.0
27 real, parameter :: ssbg_clwThresh=0.02
28 integer, parameter :: ssbg_mxval=30
29 integer, parameter :: ssbg_maxObsNum=2500
30 real, parameter :: clw_amsu_rej=0.3
31 real, parameter :: clw_amsu_rej_ch3=0.1
32 ! Highest peaking AMSU-A like SSMIS channel for ocean-only and CLW filtering
33 ! 3 = mid-troposphere (AMSU/operations -- AMSU chan. 5)
34 ! 4 = upper-troposphere (scat. index used in AMSU/operations -- AMSU chan. 6)
35 ! (AMSU-A scat. index cannot be computed here; need AMSU-A channels 1,2)
36 integer, parameter :: ipc=4
37 ! Module variable
38
39 character(len=128), parameter :: fileMgLg='fstglmg' ! glace de mer file
40 character(len=128), parameter :: fileGlace='bicefil' ! binaire 0.1degre ice file
41 character(len=128), parameter :: fileWentz='wentz_surf.std' ! surface wentz file
42 character(len=128), parameter :: algOption = 'fwentz'
43 ! Other NRL thresholds
44
45 integer, parameter :: ssbg_maxNumSat = 4
46 integer, parameter :: ssbg_maxNumChan = 24
47 integer, parameter :: ssbg_maxNumTest = 16
48
49 ! namelist variables
50 logical :: RESETQC ! reset Qc flags option
51 logical :: debug ! debug mode
52
53
54 namelist /nambgck/debug, RESETQC
55
56contains
57
58 subroutine ssbg_init()
59 !:Purpose: This subroutine reads the namelist section NAMBGCK
60 ! for the module.
61
62 implicit none
63
64 ! Locals:
65 integer :: ierr, nulnam
66 ! External functions
67 integer, external :: fclos, fnom
68
69 ! Default values for namelist variables
70 debug = .false.
71 RESETQC = .false.
72
73 nulnam = 0
74 ierr = fnom(nulnam, './flnml','FTN+SEQ+R/O', 0)
75 read(nulnam, nml=nambgck, iostat=ierr)
76 if (ierr /= 0) call utl_abort('ssbg_init: Error reading namelist')
77 if (mmpi_myid == 0) write(*, nml=nambgck)
78 ierr = fclos(nulnam)
79
80 ssbg_debug = debug
81
82 end subroutine ssbg_init
83
84 !--------------------------------------------------------------------
85 ! ssmis_tb2ta
86 !--------------------------------------------------------------------
87 subroutine ssmis_tb2ta(numObsToProcess, grossRej, ztb, zta)
88 !:Purpose: Convert Tbs received from UKMO to Tas, by reversing Ta to Tb
89 ! spillover correction applied in B. Bell's pre-processing.
90
91 implicit none
92
93 ! Arguments:
94 integer, intent(in) :: numObsToProcess ! Number of obs points to process
95 logical, intent(in) :: grossRej(:) ! Gross rejection indicator
96 real, intent(in) :: ztb(:) ! Tbs from input BURP file
97 real, intent(out) :: zta(:) ! Tas after conversion
98
99 ! Locals:
100 integer :: hiIndex
101 integer :: loIndex
102 integer :: obsIndex
103 real :: spillCoeffs(ssbg_maxNumChan) ! Spillover correction coefficients
104
105 ! Define spillover correction coefficients
106 !! Row1 ch1 ch2 ch3 ch4 ch5 ch6
107 !! Row2 ch7 ch8 ch9 ch10 ch11 ch12
108 !! Row3 ch13 ch14 ch15 ch16 ch17 ch18
109 !! Row4 ch19 ch20 ch21 ch22 ch23 ch24
110
111 ! Spillover coeff for all channels (from Steve Swadley/NRL 9 June 2010)
112
113 ! spillCoeffs = (/ 0.9850, 0.9850, 0.9850, 0.9850, 0.9850, 0.9815, &
114 ! 0.9815, 0.9949, 0.9934, 0.9934, 0.9934, 0.9680, &
115 ! 0.9720, 0.9820, 0.9810, 0.9850, 0.9820, 0.9780, &
116 ! 0.9815, 0.9815, 0.9815, 0.9815, 0.9815, 0.9815 /)
117 !
118 ! Spillover coeff for 7 SSM/I-like channels (ch. 12-18) only
119 !
120 spillCoeffs = (/ 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, &
121 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 0.9680, &
122 0.9720, 0.9820, 0.9810, 0.9850, 0.9820, 0.9780, &
123 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000 /)
124
125 ! Apply Tb -> Ta conversion
126
127 loIndex = 1
128 do obsIndex = 1, numObsToProcess
129
130 hiIndex = obsIndex*ssbg_maxNumChan
131 if ( .not. grossRej(obsIndex) ) then
132 zta(loIndex:hiIndex) = spillCoeffs(:) * ztb(loIndex:hiIndex)
133 else
134 zta(loIndex:hiIndex) = ztb(loIndex:hiIndex)
135 end if
136
137 loIndex = hiIndex + 1
138
139 end do
140
141 end subroutine ssmis_tb2ta
142
143 !--------------------------------------------------------------------------
144 ! f16tdr_remapping
145 !--------------------------------------------------------------------------
146 subroutine f16tdr_remapping(satId, SSMIS_Ta, Remapped_SSMI_Ta)
147 !:Purpose: Remap SSMIS imaging channel antenna temperature to SSMI Ta
148
149 ! SSMIS C_Freq SSMI
150 ! ------- ----------------- -------
151 ! Chan12 19.35h Chan2
152 ! Chan13 19.35v Chan1
153 ! Chan14 22.235v Chan3
154 ! Chan15 37.0h Chan5
155 ! Chan16 37.0v Chan4
156 ! Chan17 91.65v -> 85.5v Chan6
157 ! Chan18 91.65h -> 85.5h Chan7
158
159 implicit none
160
161 ! Arguments:
162 integer, intent(in) :: satId ! Satellite ID
163 real, intent(in) :: SSMIS_Ta(ssbg_maxNumChan) ! SSMIS antenna temperature
164 real, intent(out) :: Remapped_SSMI_Ta(ssbg_maxNumChan) ! Remapped SSMI antenna temperature
165
166 ! Locals:
167 integer, parameter :: f16_id = 1
168 integer, parameter :: f17_id = 2
169 integer, parameter :: f18_id = 3
170 integer(2) :: channelIndex
171 real(8) :: CP(ssbg_maxNumChan)
172 real :: tbx(ssbg_maxNumChan)
173
174 ! NESDIS Intercept/Slope for F16 SSMIS 7 IMG channels to SSMI linear remapping
175 ! ch1 ch2 ch3 ch4 ch5 ch6 ch7 ch8 ch9 ch10 ch11 (LAS/ENV)
176 real(8), parameter :: AP(ssbg_maxNumChan)=(/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
177 7.44254,7.80472,6.76383,8.55426,7.34409,6.57813,6.45397, &
178 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
179
180 real(8), parameter :: BP(ssbg_maxNumChan)=(/1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0,&
181 0.969424,0.967519,0.959808,0.954316,0.958955,0.980339,0.978795, &
182 1.0, 1.0, 1.0, 1.0, 1.0, 1.0/)
183
184 ! NESDIS F17 and F18 Tb biases with respect to F16
185 real(8), parameter :: CP_F17(ssbg_maxNumChan)=(/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
186 -0.779, -1.446, -1.013, -0.522, -0.240, 0.735, 0.521, &
187 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
188
189 real(8), parameter :: CP_F18(ssbg_maxNumChan)=(/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
190 -0.773, -0.688, -1.031, -0.632, -0.411, 0.171, 0.928, &
191 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
192
193 ! Initialization
194 tbx(1:ssbg_maxNumChan) = SSMIS_Ta(1:ssbg_maxNumChan)
195
196 if ( satId == f17_id ) then
197 CP = CP_F17
198 else if ( satId == f18_id ) then
199 CP = CP_F18
200 else
201 CP = 0.0
202 end if
203
204 do channelIndex=1, ssbg_maxNumChan
205 Remapped_SSMI_Ta(channelIndex) = AP(channelIndex) + BP(channelIndex)*(tbx(channelIndex)+CP(channelIndex))
206 end do
207
208 end subroutine f16tdr_remapping
209
210
211 !--------------------------------------------------------------------------
212 ! ssmi_ta2tb_fweng
213 !--------------------------------------------------------------------------
214 subroutine ssmi_ta2tb_fweng(Ta, Tb)
215 !:Purpose: To convert antenna temperature(Ta) to brightness temperature(Tb).
216
217 ! (1) All channel antenna gain spill-over correction
218 ! (2) Imaging channel Cross-polarization correction
219 ! (3) Doppler correction (to be developed)
220
221 implicit none
222
223 ! Arugments
224 real, intent(in) :: Ta(24) ! Antenna temperature
225 real, intent(out) :: Tb(24) ! Brightness temperature
226
227 ! Locals:
228 real(8), parameter :: AP(24)=(/0.9850,0.9850,0.9850,0.9850,0.9850,0.9790,0.9815,&
229 0.9949,0.9934,0.9934,0.9934, &
230 0.9690,0.9690,0.9740,0.9860,0.9860,0.9880,0.9880,&
231 0.9815,0.9815,0.9815,0.9815,0.9815,0.9815/)
232 real(8), parameter :: BP(24)=(/0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
233 0.0, 0.0, 0.0, 0.0, &
234 0.00415,0.00473,0.0107,0.02612,0.0217,0.01383,0.01947,&
235 0.0, 0.0, 0.0, 0.0, 0.0, 0.0/)
236
237 integer(4) :: channelIndex
238 real(4) :: CP(24), DP(24)
239
240 ! All channel antenna gain correction. Note the cross-polarization effects on antenna gain correction.
241 do channelIndex=1, 24
242 CP(channelIndex) = 1.0/( AP(channelIndex)*(1.0 - BP(channelIndex)) )
243 DP(channelIndex) = CP(channelIndex) * BP(channelIndex)
244 Tb(channelIndex) = CP(channelIndex)*Ta(channelIndex)
245 end do
246
247 ! SSMI IMG channel cross polarization correction
248 Tb(12) = Tb(12) - DP(12)*Ta(13) ! 19H
249 Tb(13) = Tb(13) - DP(13)*Ta(12) ! 19V
250 Tb(14) = Tb(14) - DP(13)*(0.65*Ta(13)+96.6) ! 22V
251 Tb(15) = Tb(15) - DP(15)*Ta(16) ! 37H
252 Tb(16) = Tb(16) - DP(16)*Ta(15) ! 37V
253 Tb(17) = Tb(17) - DP(17)*Ta(18) ! 85V
254 Tb(18) = Tb(18) - DP(18)*Ta(17) ! 85H
255
256 end subroutine ssmi_ta2tb_fweng
257
258 !--------------------------------------------------------------------------
259 ! ssmi_ta2tb_fwengtz
260 !--------------------------------------------------------------------------
261 subroutine ssmi_ta2tb_fwentz(Ta, Tb)
262 !:Purpose: Convert antenna temperatures to brightness temperatures.
263
264 implicit none
265
266 ! Arguments:
267 real, intent(in) :: Ta(:) ! Antenna temperature
268 real, intent(out) :: Tb(:) ! Brightness temperature
269
270 ! Locals:
271 real :: AV19V,AH19V,A019V,AH19H,AV19H,A019H,AV22V,A022V
272 real :: AV37V,AH37V,A037V,AH37H,AV37H,A037H
273 real :: AV85V,AH85V,A085V,AH85H,AV85H,A085H
274 real :: TA19V,TA19H,TA22V,TA37V,TA37H,TA85V,TA85H
275 real :: TB19V,TB19H,TB22V,TB37V,TB37H,TB85V,TB85H
276
277 ! Define APC coefficients
278
279 AV19V = 1.0369830
280 AH19V = -0.0039359
281 A019V = -0.0892273
282
283 AH19H = 1.0384912
284 AV19H = -0.0054442
285 A019H = -0.0892271
286
287 AV22V = 1.01993
288 A022V = 1.994
289
290 AV37V = 1.0368094
291 AH37V = -0.0222607
292 A037V = -0.0392815
293
294 AH37H = 1.0421693
295 AV37H = -0.0276206
296 A037H = -0.0392816
297
298 AV85V = 1.0263188
299 AH85V = -0.0143165
300 A085V = -0.0324062
301
302 AH85H = 1.0321901
303 AV85H = -0.0201877
304 A085H = -0.0324065
305
306 ! Extract list of Tas from parent array.
307
308 TA19H = Ta(12)
309 TA19V = Ta(13)
310 TA22V = Ta(14)
311 TA37H = Ta(15)
312 TA37V = Ta(16)
313 TA85V = Ta(17)
314 TA85H = Ta(18)
315
316 ! Apply APC corrections: convert TDR's to SDR's.
317
318 TB19V= (AV19V * TA19V) + (AH19V * TA19H) + (A019V)
319 TB19H= (AH19H * TA19H) + (AV19H * TA19V) + (A019H)
320 TB22V= (AV22V * TA22V) + A022V
321 TB37V= (AV37V * TA37V) + (AH37V * TA37H) + (A037V)
322 TB37H= (AH37H * TA37H) + (AV37H * TA37V) + (A037H)
323 TB85V= (AV85V * TA85V) + (AH85V * TA85H) + (A085V)
324 TB85H= (AH85H * TA85H) + (AV85H * TA85V) + (A085H)
325
326 ! Insert list of new Tbs into parent array.
327
328 Tb(12) = TB19H
329 Tb(13) = TB19V
330 Tb(14) = TB22V
331 Tb(15) = TB37H
332 Tb(16) = TB37V
333 Tb(17) = TB85V
334 Tb(18) = TB85H
335
336 end subroutine ssmi_ta2tb_fwentz
337
338 !--------------------------------------------------------------------------
339 ! compute_iwv_101
340 !--------------------------------------------------------------------------
341 subroutine compute_iwv_101(Tb, iwv)
342 !:Purpose: Compute integrated water vapor from SSMI brightness temperatures.
343
344 implicit none
345
346 ! Arguments:
347 real, intent(in) :: Tb(24) ! Brightness temperature
348 real, intent(out) :: iwv ! Integrated water vapor (kg/m**2)
349
350 ! Locals:
351 integer :: precipScreen ! = 1: possible presence of precipitation does
352 ! not allow retrieval of IWV and CLW.
353 ! = 0: retrieval is possible
354 real :: IWV_alishouse ! estimated total precipitable water, with cubic polynomial correction (=iwv)
355 real :: IWV_alishouse0 ! estimated total precipitable water, without cubic polynomial correction
356 real :: precipThresh
357 real :: tb19v,tb22v,tb37v,tb37h
358 real :: ciwv(0:4)
359
360 ! Initializations.
361
362 ciwv = (/232.89393,-0.148596,-1.829125,-0.36954,0.006193/)
363
364 tb19v = Tb(13)
365 tb22v = Tb(14)
366 tb37h = Tb(15)
367 tb37v = Tb(16)
368
369 ! Compute IWV_alishouse.
370
371 IWV_alishouse0 = ciwv(0) + ciwv(1)*tb19v + ciwv(2)*tb22v + &
372 & ciwv(3)*tb37v + ciwv(4)*tb22v**2
373
374 ! Apply Cubic Polynomial Correction to the CAL/VAL water vapour
375 ! algorithm (G.W. Petty).
376
377 IWV_alishouse = -3.75 + 1.507*IWV_alishouse0 - 0.01933*IWV_alishouse0**2 + &
378 & (2.191e-04)*IWV_alishouse0**3
379
380 ! Compute Precipitation Screen.
381
382 precipThresh = -11.7939 - 0.02727*tb37v + 0.09920*tb37h
383 if ( precipThresh > 0.0 ) then
384 precipScreen = 1
385 else
386 precipScreen = 0
387 end if
388
389 ! Apply precipitation screen to IWV_alishouse (units are kg/m**2).
390 ! Bounds now applied in cld_filter_fweng.ftn90 routinge.
391
392 if( precipScreen == 1 ) then
393 IWV_alishouse = ssbg_rmisg
394 end if
395
396 ! Store the IWV value.
397
398 iwv = IWV_alishouse
399
400 end subroutine compute_iwv_101
401
402 !--------------------------------------------------------------------------
403 ! determ_tpw
404 !--------------------------------------------------------------------------
405 subroutine determ_tpw(Tb, sType, seaIce, TPW)
406 !:Purpose: To calculate total precipitable water (in mm).
407
408 implicit none
409
410 ! Arguments:
411 real, intent(in) :: Tb(24) ! Brightness temperature
412 integer, intent(in) :: sType ! Surface type
413 real, intent(in) :: seaIce ! Sea ice coverage
414 real, intent(out) :: TPW ! Total precipitable water (mm)
415
416 ! Locals:
417 integer, parameter :: ocean=0
418 real :: SCT
419 real :: Tb19v, Tb22v, Tb37v, Tb85v
420
421 ! Extract 4 IMG channels from Tb
422 Tb19v = Tb(13)
423 Tb22v = Tb(14)
424 Tb37v = Tb(16)
425 Tb85v = Tb(17)
426
427 TPW = ssbg_rmisg
428
429 ! No TWP over land and sea ice
430 if ( sType == ocean ) then
431
432 SCT = -182.7 + 0.75*Tb19v + 2.543*Tb22v - 0.00543*Tb22v*Tb22v - Tb85v
433
434 if ( abs(seaIce) < 70.0 ) then
435 TPW = 232.89393 - 0.148596*Tb19v - 1.829125*Tb22v + 0.006193*Tb22v**2 - 0.36954*Tb37v
436 TPW = -3.753 + 1.507*TPW - 0.01933*TPW**2 + 0.0002191*TPW**3
437
438 if ( TPW < 0.0 ) TPW = 0.0
439 if ( TPW > 80.0 ) TPW = 80.0
440 end if
441
442 end if
443
444 end subroutine determ_tpw
445
446 !--------------------------------------------------------------------------
447 ! determ_sea_ice
448 !--------------------------------------------------------------------------
449 subroutine determ_sea_ice(ocean, Ta, sType, seaIce, latitude)
450 !:Purpose: To calculate sea ice cover (in %).
451
452 implicit none
453
454 ! Arguments:
455 integer, intent(in) :: ocean ! Ocean surface type index
456 real, intent(in) :: Ta(24) ! Antenna temperature
457 integer, intent(in) :: sType ! Surface type
458 real, intent(out) :: seaIce ! Sea ice coverage
459 real, intent(in) :: latitude ! Latitude of observation
460
461 ! Locals:
462 real :: Ta19v, Ta19h, Ta22v, Ta37v, Ta37h, Ta85v, Ta85h
463
464 ! Extract 7 IMG channels from Ta
465 Ta19v = Ta(13)
466 Ta19h = Ta(12)
467 Ta22v = Ta(14)
468 Ta37v = Ta(16)
469 Ta37h = Ta(15)
470 Ta85v = Ta(17)
471 Ta85h = Ta(18)
472
473 ! Calculate Sea Ice Coverage
474 if ( sType == ocean ) then ! Over Ocean
475 if (latitude > 44.4 .or. latitude < -52.0 ) then
476 seaIce = 91.9 - 2.994*Ta22v + 2.846*Ta19v - 0.386*Ta37v + 0.495*Ta85v &
477 + 1.005*Ta19h - 0.904*Ta37h
478 if ( seaIce >= 70.0 ) then
479 seaIce = 100.0
480 else
481 seaIce = 0.0
482 end if
483 else
484 seaIce = 0.0
485 end if
486 else ! Over Land, there is no sea ice!
487 seaIce = ssbg_rmisg
488 end if
489
490 end subroutine determ_sea_ice
491
492 !--------------------------------------------------------------------------
493 ! determ_clw
494 !--------------------------------------------------------------------------
495 subroutine determ_clw(algOption, Ta, Tb, sType, CLW, IWV, latitude)
496 !:Purpose: To calculate cloud liquid water for a single data point (in kg/m**2).
497
498 ! Normally called when sType = 0 (open water).
499 ! Also retrieves sea-ice to see if there is sea-ice at water point.
500
501 implicit none
502
503 ! Arguments:
504 character(len=6), intent(in) :: algOption ! Algorithm option (fweng, fwentz or nsun)
505 real, intent(in) :: Ta(24) ! Antenna temperature
506 real, intent(in) :: Tb(24) ! Brightness temperature
507 integer, intent(inout) :: sType ! Surface type
508 real, intent(out) :: CLW ! Cloud liquid water (in kg/m**2)
509 real, intent(inout) :: IWV ! Integrated water vapor
510 real, intent(in) :: latitude ! Latitude of observation
511
512 ! Locals:
513 integer, parameter :: isSeaIce = 1
514 integer, parameter :: ocean = 0
515 integer(4), parameter :: RT = 285
516 real, parameter :: clwLimit = 6.0
517 real :: ALG1, ALG2, ALG3
518 real :: seaIce
519 real :: Ta19v, Ta19h, Ta22v, Ta37v, Ta37h, Ta85v, Ta85h
520 real :: Tb19v, Tb22v, Tb37v
521 real :: TPW
522
523 CLW = ssbg_rmisg
524
525 ! Extract 7 IMG channels from Ta and Tb
526 Ta19v = Ta(13)
527 Ta19h = Ta(12)
528 Ta22v = Ta(14)
529 Ta37v = Ta(16)
530 Ta37h = Ta(15)
531 Ta85v = Ta(17)
532 Ta85h = Ta(18)
533
534 Tb19v = Tb(13)
535 Tb22v = Tb(14)
536 Tb37v = Tb(16)
537
538 ! Call determ_sea_ice to find the seaIce
539 ! -- seaIce = 100.0 when sea ice >= 70%
540 ! = 0.0 when sea ice < 70%
541 call determ_sea_ice(ocean, Ta, sType, seaIce, latitude)
542
543 ! Calculate CLW Over Ocean
544 if ( (sType == ocean) .and. (seaIce /= 100.0) ) then
545 ALG1 = ssbg_rmisg
546 ALG2 = ssbg_rmisg
547 ALG3 = ssbg_rmisg
548
549 if ( trim(algOption) == 'fweng' ) then
550 ! Compute IWV using F. Weng algorithm.
551 IWV = 232.89 - 0.1486*Tb19v - 0.3695*Tb37v - (1.8291 - 0.006193*Tb22v)*Tb22v
552 if ( IWV < 0.0 ) IWV = 0.0
553 end if
554 if ( trim(algOption) == 'nsun') then
555 call determ_tpw(Tb,sType,seaIce,TPW)
556 IWV = TPW
557 end if
558
559 if ( (Ta19v < RT) .and. (Ta22v < RT) ) then
560 ALG1 = -3.20 * ( ALOG(290.0-Ta19v) - 2.80 - 0.42*ALOG(290.0-Ta22v) ) !TA
561 ! ALG1 = -3.20 * ( ALOG(290.0-Tb19v) - 2.84 - 0.40*ALOG(290.0-Tb22v) ) !TB
562 end if
563
564 if ( (Ta37v < RT) .and. (Ta22v < RT) ) then
565 ALG2 = -1.66 * ( ALOG(290.0-Ta37v) - 2.90 - 0.349*ALOG(290.0-Ta22v) ) !TA
566 ! ALG2 = -1.66 * ( ALOG(290.0-Tb37v) - 2.99 - 0.32*ALOG(290.0-Tb22v) ) !TB
567 end if
568
569 if ( (Ta85h < RT) .and. (Ta22v < RT) ) then
570 ALG3 = -0.44 * ( ALOG(290.0-Ta85h) + 1.60 - 1.354*ALOG(290.0-Ta22v) ) !TA
571 ! ALG3 = -0.44 * ( ALOG(290.0-Tb85h) + 1.11 - 1.26*ALOG(290.0-Tb22v) ) !TB
572 end if
573
574 if ( ALG1 > 0.70 ) then
575 CLW = ALG1
576 else if ( ALG2 > 0.28 ) then
577 CLW = ALG2
578 else if ( IWV < 30.0 ) then
579 CLW = ALG3
580 else
581 CLW = ALG2
582 end if
583
584 ! Verify CLW is within acceptable upper limit.
585 if ( CLW > clwLimit ) CLW = ssbg_rmisg
586
587 ! Force negative CLW values to zero.
588 if ( CLW < 0.0 .and. CLW /= ssbg_rmisg ) CLW = 0.0
589
590 else
591 ! Sea Ice (>70%) detected from s/r determ_sea_ice but sType was 0 = waterobs (on call)
592 CLW = -500.0
593 IWV = ssbg_rmisg
594 sType = isSeaIce
595
596 end if
597
598 end subroutine determ_clw
599
600 !--------------------------------------------------------------------------
601 ! cld_filter_fweng
602 !--------------------------------------------------------------------------
603 subroutine cld_filter_fweng(numObsToProcess, obsTb, algOption, waterObs, grossRej, &
604 & cloudObs, iwvReject, precipObs, rclw, riwv, iSatId, &
605 & obsLatitude, numSeaIceObs)
606 !:Purpose: Compute the cloud liquid water (CLW) from SSMIS channels using the
607 ! regression algorithm of Fuzhong Weng and Ninghai Sun.
608 ! Retrieve CLW path from F16 SSMIS TDR data
609
610 implicit none
611
612 ! Arguments:
613 integer, intent(in) :: numObsToProcess ! Number of obs points to process
614 real, intent(in) :: obsTb(:) ! Brightness temperature of observations
615 character(len=6), intent(in) :: algOption ! Algorithm option (fweng, fwentz or nsun)
616 logical, intent(inout) :: waterObs(:) ! Open water identifier for each obs
617 logical, intent(in) :: grossRej(:) ! Logical array of obs with gross error (obs to reject)
618 logical, intent(inout) :: cloudObs(:) ! Logical array of obs for which CLW > 0.01 kg/m**2 or with precipitations
619 logical, intent(inout) :: iwvReject(:) ! Logical array of obs for which IWV > 80 kg/m**2
620 logical, intent(inout) :: precipObs(:) ! Logical array of obs with precipitations (CLW missing)
621 real, intent(inout) :: rclw(:) ! Real array of CLW
622 real, intent(inout) :: riwv(:) ! Real array of integrated water vapor (IWV)
623 integer, intent(in) :: iSatId ! Satellite identifier
624 real, intent(in) :: obsLatitude(:) ! Observation latitudes
625 integer, intent(inout) :: numSeaIceObs ! Number of observations with sea ice
626
627 ! Locals:
628 real, parameter :: iwvThresh = 80.0 ! Upper bound for IWV in kg/m**2
629 integer, parameter :: ocean = 0
630 integer, parameter :: seaIce = 1
631 integer :: hiIndex
632 integer :: loIndex
633 integer :: obsIndex
634 integer :: sType
635 real :: clw
636 real :: iwv
637 real :: latitude
638 real :: F16TDR(ssbg_maxNumChan)
639 real :: remappedTa(ssbg_maxNumChan)
640 real :: Tb(ssbg_maxNumChan)
641 real :: zta(ssbg_mxval*ssbg_maxObsNum)
642
643 !---------------------------------------------------
644
645 ! Convert Tbs received from UKMO to Tas, by reversing Ta to Tb
646 ! spillover correction applied in B. Bell's pre-processing.
647 ! Missing Tbs are also checked for here.
648 ! Apply to all obs of current record.
649
650 call ssmis_tb2ta(numObsToProcess,grossRej,obsTb,zta)
651
652 rclw(:) = ssbg_rmisg
653 riwv(:) = ssbg_rmisg
654
655 loIndex = 1
656 HEADER: do obsIndex = 1, numObsToProcess
657
658 hiIndex = obsIndex*ssbg_maxNumChan
659 F16TDR(:) = zta(loIndex:hiIndex)
660 latitude = obsLatitude(obsIndex)
661 loIndex = hiIndex + 1
662
663 ! Obtain CLW and IWV for obs points over open ocean, and where
664 ! Tb values have passed gross filter check.
665 ! Initialize IWV and CLW to missing for all other cases.
666
667 clw = ssbg_rmisg
668 iwv = ssbg_rmisg
669
670 if ( waterObs(obsIndex) .and. (.not. grossRej(obsIndex)) ) then
671
672 sType = ocean ! Surface-type=Ocean
673
674 ! Call SSMIS TDR to SSMI TDR remapping subroutine
675
676 call f16tdr_remapping(iSatId, F16TDR, remappedTa)
677
678 ! Call SSM/I Ta to Tb conversion subroutine
679
680 if ( trim(algOption) == 'fweng' ) then
681 call ssmi_ta2tb_fweng(remappedTa, Tb)
682 ! IWV computed in determ_clw subroutine below.
683 else if ( trim(algOption) == 'fwentz' ) then
684 call ssmi_ta2tb_fwentz(remappedTa, Tb)
685 ! Compute IWV using Alishouse and Petty,
686 ! because it won't be computed in determ_clw.
687 ! Missing value for IWV means possible precipitation
688 ! and so CLW will not be computed
689 call compute_iwv_101(Tb,iwv)
690 if ( iwv == ssbg_rmisg ) precipObs(obsIndex) = .true.
691 else if ( trim(algOption) == 'nsun' ) then
692 call ssmi_ta2tb_fweng(remappedTa, Tb)
693 ! IWV computed in determ_clw subroutine below.
694 else
695 write(*,*) ' cld_filter_fweng: Invalid algorithm option !! '
696 call abort()
697 end if
698
699 ! Call CLW retrieval algorithm subroutine.
700 ! -- also computes and returns (output) IWV if algOption=fweng or nsun
701 ! -- sType is also changed to seaIce value if sea-ice is detected
702
703 if ( trim(algOption) /= 'fwentz' ) then
704 call determ_clw(algOption, remappedTa, Tb, sType, clw, iwv, latitude)
705 else
706 if ( .not. precipObs(obsIndex) ) then
707 call determ_clw(algOption, remappedTa, Tb, sType, clw, iwv, latitude)
708 end if
709 end if
710
711 ! Check for newly detected sea-ice (sType changed)
712
713 if ( sType == seaIce ) numSeaIceObs = numSeaIceObs + 1
714
715 ! Reject obs with precipitation or cloud amount more than threshold.
716 ! Also, reject obs if CLW is missing (undetermined)
717 ! Set waterObs flag to false if deter_clw returns "seaIce" value (-500)
718 if ( (clw > ssbg_clwThresh) .or. precipObs(obsIndex) ) cloudObs(obsIndex) = .true.
719 if ( clw == ssbg_rmisg ) cloudObs(obsIndex) = .true.
720 if ( clw == -500.0 ) then
721 waterObs(obsIndex) = .false.
722 clw = ssbg_rmisg
723 end if
724 ! Reject obs with IWV value more than threshold and set IWV to missing.
725 ! First, set IWV values below zero to 0.0.
726 if ( iwv < 0.0 .and. iwv /= ssbg_rmisg ) iwv = 0.0
727 if ( iwv > iwvThresh ) then
728 if ( .not. cloudObs(obsIndex) ) iwvReject(obsIndex) = .true.
729 iwv = ssbg_rmisg
730 end if
731
732 ! Store CLW and IWV
733 if (ssbg_debug) then
734 write(*,*)'cld_filter_fweng: CLOUD BY DETERM_CLW = ', clw
735 write(*,*)'cld_filter_fweng: IWV BY DETERM_CLW = ', iwv
736 end if
737 rclw(obsIndex) = clw
738 riwv(obsIndex) = iwv
739
740 end if
741
742 end do HEADER
743
744 120 format(' obsIndex ',2x,i3,2x,' clw ',f4.2,2x,' CLW threshold = ',f4.2,2x,' cloudObs ',l2)
745 130 format(' obsIndex ',2x,i3,2x,' iwv ',f4.2,2x,' IWV threshold = ',f4.2,2x,' iwvReject ',l2)
746
747 end subroutine cld_filter_fweng
748
749 !--------------------------------------------------------------------------
750 ! copy1Dimto2DimRealArray
751 !--------------------------------------------------------------------------
752 subroutine copy1Dimto2DimRealArray(oneDimArray, firstDim, secondDim, twoDimArray)
753 !
754 !:Purpose: copy 1D real array into 2D real array given firstDim and secondDim
755 !
756 implicit none
757
758 ! Arguments:
759 integer, intent(in) :: firstDim ! First dimension
760 integer, intent(in) :: secondDim ! Second dimension
761 real, intent(in) :: oneDimArray(firstDim*secondDim) ! 1D real array
762 real, intent(inout) :: twoDimArray(firstDim,secondDim) ! 2D real array
763
764 ! Locals:
765 integer :: firstDimIndex
766 integer :: productDimIndex
767 integer :: secondDimIndex
768
769 ! copy the original input 1D array to 2D array. The 2D arrays are used in this s/r.
770 do secondDimIndex=1,secondDim
771 do firstDimIndex=1,firstDim
772 productDimIndex = (secondDimIndex-1)*firstDim + firstDimIndex
773 twoDimArray(firstDimIndex,secondDimIndex) = oneDimArray(productDimIndex)
774 end do
775 end do
776
777 end subroutine copy1Dimto2DimRealArray
778
779 !--------------------------------------------------------------------------
780 ! copy1Dimto2DimIntegerArray
781 !--------------------------------------------------------------------------
782 subroutine copy1Dimto2DimIntegerArray(oneDimArray, firstDim, secondDim, twoDimArray)
783 !
784 !:Purpose: copy 1D integer array into 2D Integer array given firstDim and secondDim
785 !
786 implicit none
787
788 ! Arguments:
789 integer, intent(in) :: firstDim ! First dimension
790 integer, intent(in) :: secondDim ! Second dimension
791 integer, intent(in) :: oneDimArray(firstDim*secondDim) ! 1D integer array
792 integer, intent(inout) :: twoDimArray(firstDim,secondDim) ! 2D integer array
793
794 ! Locals:
795 integer :: firstDimIndex
796 integer :: productDimIndex
797 integer :: secondDimIndex
798
799 ! copy the original input 1D array to 2D array. The 2D arrays are used in this s/r.
800 do secondDimIndex=1,secondDim
801 do firstDimIndex=1,firstDim
802 productDimIndex = (secondDimIndex-1)*firstDim + firstDimIndex
803 twoDimArray(firstDimIndex,secondDimIndex) = oneDimArray(productDimIndex)
804 end do
805 end do
806
807 end subroutine copy1Dimto2DimIntegerArray
808
809 !------------------------------------------------------------------------------------
810 ! bennartz
811 !------------------------------------------------------------------------------------
812 subroutine bennartz (ier, numObsToProcess, tb89, tb150, satZenithAngle, landSeaQualifier, scatL, scatW)
813 !:Purpose: Compute the following parameters using 2 AMSU-B channels:
814 ! - scattering index (over land and ocean).
815 ! The two channels used are: 89Ghz, 150Ghz.
816 !
817 ! REFERENCES: Bennartz, R., A. Thoss, A. Dybbroe and D. B. Michelson,
818 ! 1999: Precipitation Analysis from AMSU, Nowcasting SAF,
819 ! Swedish Meteorologicali and Hydrological Institute,
820 ! Visiting Scientist Report, November 1999.
821 !
822 implicit none
823
824 ! Arguments:
825 integer, intent(in) :: numObsToProcess ! Number of obs points to process
826 integer, intent(out) :: ier(numObsToProcess) ! Error return code
827 real, intent(in) :: tb89(:) ! 89Ghz AMSU-B brightness temperature (K)
828 real, intent(in) :: tb150(:) ! 150Ghz AMSU-B brightness temperature (K)
829 real, intent(in) :: satZenithAngle(:) ! Satellite zenith angle (deg.)
830 integer, intent(in) :: landSeaQualifier(:) ! Land/sea indicator (0=land; 1=ocean)
831 real, intent(out) :: scatL(:) ! Scattering index over land
832 real, intent(out) :: scatW(:) ! Scattering index over water
833
834 ! Locals:
835 ! Notes: In the case where an output parameter cannot be calculated, the
836 ! value of this parameter is the missing value, i.e. -99.
837 integer :: obsIndex
838
839 ! ____1) Initialise output parameters:** -------------------------------*
840
841 do obsIndex = 1, numObsToProcess
842 scatL(obsIndex) = ssbg_realMissing
843 scatW(obsIndex) = ssbg_realMissing
844 end do
845
846 !____2) Validate input parameters:** -----------------------------*
847 do obsIndex = 1, numObsToProcess
848 if ( tb89(obsIndex) < 120. .or. &
849 tb89(obsIndex) > 350. .or. &
850 tb150(obsIndex) < 120. .or. &
851 tb150(obsIndex) > 350. .or. &
852 satZenithAngle(obsIndex) < -90. .or. &
853 satZenithAngle(obsIndex) > 90. .or. &
854 landSeaQualifier(obsIndex) < 0 .or. &
855 landSeaQualifier(obsIndex) > 1 ) then
856 ier(obsIndex) = 1
857 else
858 ier(obsIndex) = 0
859 end if
860 end do
861
862 !____3) Compute parameters:** ----------------------*
863 do obsIndex = 1, numObsToProcess
864 if ( ier(obsIndex) == 0 ) then
865 if (landSeaQualifier(obsIndex) == 1 ) then
866 scatW(obsIndex) = (tb89(obsIndex)-tb150(obsIndex)) - &
867 (-39.2010+0.1104*satZenithAngle(obsIndex))
868 else
869 scatL(obsIndex) = (tb89(obsIndex)-tb150(obsIndex)) - &
870 (0.158+0.0163*satZenithAngle(obsIndex))
871 end if
872 else if ( (ier(obsIndex) /= 0 ) .and. (obsIndex <= 100 ) .and. (ssbg_debug)) then
873 write(*,*), 'bennartz: Input Parameters are not all valid: '
874 write(*,*), ' obsIndex,tb89(obsIndex),tb150(obsIndex),satZenithAngle(obsIndex),landSeaQualifier(obsIndex) = ', &
875 obsIndex,tb89(obsIndex),tb150(obsIndex),satZenithAngle(obsIndex),landSeaQualifier(obsIndex)
876 write(*,*), ' ier(obsIndex),scatL(obsIndex),scatW(obsIndex)=', &
877 ier(obsIndex),scatL(obsIndex),scatW(obsIndex)
878 end if
879 end do
880
881 end subroutine bennartz
882
883 !--------------------------------------------------------------------------
884 ! ssbg_readGeophysicFieldsAndInterpolate
885 !--------------------------------------------------------------------------
886 subroutine ssbg_readGeophysicFieldsAndInterpolate(obsLatitude, obsLongitude, modelInterpTer)
887 !
888 !:Purpose: Reads geophysical model variable (GZ) and saves for the first time.
889 ! GZ is geopotential height (GZ at surface = surface height in dam).
890 ! Then interpolates those variables to observation location.
891 !
892 implicit none
893
894 ! Arguments:
895 real, intent(in) :: obsLatitude(:) ! Observation latitudes
896 real, intent(in) :: obsLongitude(:) ! Observation longitudes
897 real, allocatable, intent(out) :: modelInterpTer(:) ! Filtered and interpolated topography (in m)
898
899 ! Locals:
900 integer, parameter :: mxLat = 5
901 integer, parameter :: mxLon = 5
902 real, parameter :: dLat = 0.4
903 real, parameter :: dLon = 0.6
904 character(len=12) :: etikxx
905 character(len=1) :: grtyp
906 character(len=4) :: nomvxx
907 character(len=2) :: typxx
908 integer :: boxPointIndex
909 integer :: boxPointNum
910 integer :: dataIndex
911 integer :: dataNum
912 integer :: ezQkDef
913 integer :: ezSetOpt
914 integer :: gdllsval
915 integer, save :: gdgz ! topo interpolation param
916 integer :: idum1, idum2, idum3, idum4
917 integer :: idum5, idum6, idum7, idum8
918 integer :: idum9, idum10, idum11, idum12
919 integer :: idum13, idum14, idum15
920 integer :: idum16, idum17, idum18
921 integer :: ier
922 integer :: ig1, ig2, ig3, ig4
923 integer :: irec
924 integer :: iUnGeo
925 integer :: latIndex
926 integer :: lonIndex
927 integer :: ni, nj, nk
928 integer :: nLat
929 integer :: nLon
930 integer :: nObsLat
931 integer :: nObsLon
932 logical, save :: ifFirstCall = .true. ! If .True. we read GL, GZ and MG
933 real, allocatable, save :: GZ(:) ! Modele Topographie (GZ)
934 real, allocatable :: GZIntBox(:,:)
935 real, allocatable :: obsLatBox (:,:)
936 real, allocatable :: obsLonBox (:,:)
937 real :: topoFact ! Facteur x topo pour avoir des unites en metre
938 real :: xLat
939 real :: xLon
940 ! External functions
941 integer, external :: fclos
942 integer, external :: fnom
943 integer, external :: fstfrm
944 integer, external :: fstinf
945 integer, external :: fstlir
946 integer, external :: fstouv
947 integer, external :: fstprm
948 integer, external :: ip1_all
949
950 ! STEP 1: CHECK if obsLatitude AND obsLongitude ARE SAME DIMENSION
951 nObsLat = size(obsLatitude)
952 nObsLon = size(obsLongitude)
953 if (nObsLat /= nObsLon) then
954 call utl_abort ('ssbg_readGeophysicFieldsAndInterpolate: OBSERVATION obsLatitude and obsLongitude should have SAME LENGTH')
955 else
956 dataNum = nObsLat
957 end if
958
959 ! STEP 2: READ GZ from the FST FILE
960 if(ifFirstCall) then
961 iUnGeo = 0
962 ier = fnom(iUnGeo,'trlm_01','STD+RND+R/O',0)
963 ier = fstouv(iUnGeo,'RND')
964
965 ! Using hybrid coordinates
966 irec = fstinf(iUnGeo,ni,nj,nk,-1,' ',ip1_all(1.0,5),-1,-1,' ','GZ')
967 if (irec < 0) then
968 call utl_abort('ssbg_readGeophysicFieldsAndInterpolate: LA TOPOGRAPHIE EST INEXISTANTE')
969 end if
970 topoFact = 10.0 ! dam --> m
971
972 if (allocated(GZ)) deallocate(GZ)
973 allocate ( GZ(ni*nj), stat=ier)
974 if ( ier /= 0 ) then
975 call utl_abort('ssbg_readGeophysicFieldsAndInterpolate: Allocation of array GZ failed')
976 end if
977 ier = fstlir(GZ,iUnGeo,ni,nj,nk,-1,' ',ip1_all(1.0,5),-1,-1,' ','GZ')
978
979 GZ(:) = GZ(:)*topoFact
980
981 ier = fstprm ( irec, idum1, idum2, idum3, idum4, &
982 idum5, idum6, idum7, idum8, idum9, idum10, &
983 idum11, typxx, nomvxx, etikxx, grtyp, ig1, &
984 ig2, ig3, ig4, idum12, idum13, idum14, &
985 idum15, idum16, idum17, idum18 )
986 write (*,*) ' GRILLE GZ : ',grtyp,ni,nj, &
987 ig1,ig2,ig3,ig4
988 ier = ezSetOpt('INTERP_DEGREE','LINEAR')
989 gdgz = ezQkDef(ni,nj,grtyp,ig1,ig2,ig3,ig4,iUnGeo)
990
991 ier = fstfrm(iUnGeo)
992 ier = fclos(iUnGeo)
993 ifFirstCall = .false.
994 end if
995
996 ! STEP 3: Interpolation de la glace et le champ terre/mer du modele aux pts TOVS.
997 ! N.B.: on examine ces champs sur une boite centree sur chaque obs.
998 boxPointNum = mxLat*mxLon
999 if(allocated(obsLatBox)) deallocate(obsLatBox)
1000 allocate (obsLatBox(boxPointNum, dataNum) , stat=ier)
1001 if(allocated(obsLonBox)) deallocate(obsLonBox)
1002 allocate (obsLonBox(boxPointNum, dataNum) , stat=ier)
1003 if(allocated(GZIntBox)) deallocate(GZIntBox)
1004 allocate (GZIntBox(boxPointNum, dataNum) , stat=ier)
1005
1006 nLat = (mxLat-1)/2
1007 nLon = (mxLon-1)/2
1008 do dataIndex = 1, dataNum
1009 boxPointIndex = 0
1010 do latIndex = -nLat, nLat
1011 xLat = obsLatitude(dataIndex) +latIndex*dLat
1012 xLat = max(-90.0,min(90.0,xLat))
1013 do lonIndex = -nLon, nLon
1014 boxPointIndex = boxPointIndex + 1
1015 xLon = obsLongitude(dataIndex) +lonIndex*dLon
1016 if ( xLon < -180. ) xLon = xLon + 360.
1017 if ( xLon > 180. ) xLon = xLon - 360.
1018 if ( xLon < 0. ) xLon = xLon + 360.
1019 obsLatBox(boxPointIndex,dataIndex) = xLat
1020 obsLonBox(boxPointIndex,dataIndex) = xLon
1021 end do
1022 end do
1023 end do
1024
1025 ier = gdllsval(gdgz,GZIntBox,GZ,obsLatBox,obsLonBox,boxPointNum*dataNum)
1026 if (ier < 0) then
1027 call utl_abort ('ssbg_readGeophysicFieldsAndInterpolate: ERROR in the interpolation of GZ')
1028 end if
1029
1030 if(allocated(modelInterpTer)) deallocate(modelInterpTer)
1031 allocate (modelInterpTer(dataNum) , stat=ier)
1032
1033 modelInterpTer(:) = 0.0
1034 do dataIndex = 1, dataNum
1035 if (ssbg_debug) then
1036 write(*,*), 'ssbg_readGeophysicFieldsAndInterpolate: infos'
1037 write(*,*), ' '
1038 write(*,*), ' dataIndex = ', dataIndex
1039 write(*,*), ' '
1040 write(*,*), ' obsLatitude, obsLongitude = ', obsLatitude(dataIndex), obsLongitude(dataIndex)
1041 write(*,*), ' '
1042 write(*,*), ' obsLatBox = '
1043 write(*,*), (obsLatBox(boxPointIndex,dataIndex),boxPointIndex=1,boxPointNum)
1044 write(*,*), ' obsLonBox = '
1045 write(*,*), (obsLonBox(boxPointIndex,dataIndex),boxPointIndex=1,boxPointNum)
1046 write(*,*), ' GZIntBox = '
1047 write(*,*), (GZIntBox(boxPointIndex,dataIndex),boxPointIndex=1,boxPointNum)
1048 end if
1049 do boxPointIndex=1, boxPointNum
1050 modelInterpTer(dataIndex) = max(modelInterpTer(dataIndex),GZIntBox(boxPointIndex,dataIndex))
1051 end do
1052 if (ssbg_debug) then
1053 write(*,*), 'ssbg_readGeophysicFieldsAndInterpolate: modelInterpTer = ', modelInterpTer(dataIndex)
1054 end if
1055 end do
1056
1057 end subroutine ssbg_readGeophysicFieldsAndInterpolate
1058
1059 !--------------------------------------------------------------------
1060 ! land_ice_mask_ssmis
1061 !--------------------------------------------------------------------
1062 subroutine land_ice_mask_ssmis(numObsToProcess, obsLatitude, obsLongitude, landSeaQualifier, &
1063 & terrainType, waterObs)
1064 !:Purpose: Determine for each observation point the ice mask value from
1065 ! the binary file copied to the local work directory.
1066
1067 ! This may be a user-specified file or it is copied from
1068 ! /users/dor/afsi/sio/datafiles/data2/ade.maskice10
1069 ! The ice mask is updated every day at 00 UTC. The binary file has
1070 ! a resolution of 0.1 deg. Observations with an ice mask value of 0
1071 ! (=ice; land or sea=-1) are removed. The ice mask value also
1072 ! determines the terrain-type qualifier (element 13039) for each obs
1073 ! pt, which is required when writing output to boxed format.
1074
1075 ! Finally, this routine performs a land/ice consistency check using
1076 ! the MG and LG fields used by the model which produces the trial field.
1077 ! The purpose of this check is to remove obs that reside close to coasts
1078 ! or ice, and so whose TBs may be contaminated.
1079 ! For the GEM global model, any of the analysis files provides the MG
1080 ! field, while for the meso-global model a user-specified file is required
1081 ! to define MG. In either case, the GEM Global analysis file provides the
1082 ! LG field, and it is interpolated to the resolution of the model.
1083 ! - MG from climate file : /data/gridpt/dbase/anal/glbeta2/DATE_000
1084 ! - LG from analysis file: ex. /data/gridpt/dbase/anal/glbpres2/DATE_000
1085
1086 ! In the application of this check, a 5x5 mesh, with spacing defined by rLatKm and
1087 ! rLonKm, is positioned with its center over an obs pt (2 grid pts on either side
1088 ! of the obs pt; size of mesh is equal to 4*rLatKm x 4*rLonKm). The values of MG
1089 ! and LG are evaluated at the grid points of this mesh. The maximum value of each
1090 ! determines whether the obs pt is too close to ice or land to be retained.
1091 ! **NOTE: the threshold value for MG has a very strong effect on the distance
1092 ! from land that is permitted for an obs to be retained
1093
1094
1095 ! Maximum FOV x---x---x---x---x ^
1096 ! = 75km x 75km | | | | | |
1097 ! for Meso-sphere CHs x---x---x---x---x |
1098 ! = 74km x 47km | | | | | |
1099 ! for 19 GHz x---x---o---x---x | = 4*rLatKm
1100 ! | | | | | | = 4*40 km
1101 ! ^ x---x---x---x---x | = 160 km = 80 km north & south
1102 ! rLatKm | | | | | | |
1103 ! v x---x---x---x---x v
1104 ! <--->
1105 ! rLonKm
1106 !
1107 ! <--------------->
1108 ! = 4*rLonKm
1109 ! = 4*40 km
1110 ! = 160 km = 80 km east & west
1111
1112
1113 ! MG value = 1.0 ==> LAND MG value = 0.0 ==> OCEAN
1114 ! LG value = 1.0 ==> ICE LG value = 0.0 ==> NO ICE
1115
1116 !--------------------------------------------------------------------
1117 ! Variable Definitions
1118 ! --------------------
1119 ! fileMgLg - input - name of file holding model MG and LG fields (external)
1120 ! numObsToProcess - input - number of input obs pts in record
1121 ! obsLatitude - input - array holding lat values for all obs pts in record
1122 ! obsLongitude - input - array holding lon values for all obs pts in record
1123 ! landSeaQualifier - in/out - array holding land/sea qualifier values for all obs
1124 ! pts of record (0 = land, 1 = sea)
1125 ! terrainType - output - array holding terrain-type values for all obs pts
1126 ! of current record
1127 ! waterObs - output - logical array identifying for each obs in current record
1128 ! whether it is over open water, far from coast/ice
1129 ! iMask -internal- value determined by interpolating obs pt to
1130 ! binary ice mask field
1131 ! for land/sea: iMask = -1
1132 ! for ice: iMask = 0
1133 ! mxLat -internal- number of grid pts in lat. direction for mesh
1134 ! mxLon -internal- number of grid pts in lon. direction for mesh
1135 ! rLatKm -internal- spacing desired between mesh grid points in km
1136 ! along lat. direction
1137 ! rLonKm -internal- spacing desired between mesh grid points in km
1138 ! along lon. direction
1139 ! dLat -internal- spacing between mesh grid points along lon. direction
1140 ! in degrees computed from rLatKm
1141 ! dLon -internal- spacing between mesh grid points along lon. direction
1142 ! in degrees computed from rLonKm
1143 ! rKmPerDeg -internal- distance in km per degree
1144 ! = Earth radius * PI/180.0
1145 ! = 6371.01 km * PI/180.0
1146 ! = 111.195 km
1147 ! nLat,nLon -internal- used to define the lat/lon of the grid pts of mesh
1148 ! obsLatBox -internal- lat values at all grid pts of mesh for all obs pts
1149 ! obsLonBox -internal- lon values at all grid pts of mesh for all obs pts
1150 ! latMesh -internal- lat values at all grid pts of mesh for 1 obs pt
1151 ! lonMesh -internal- lon values at all grid pts of mesh for 1 obs pt
1152 ! mgIntOb -internal- interpolated MG values at all grid pts of mesh for 1 obs pt
1153 ! lgIntOb -internal- interpolated LG values at all grid pts of mesh for 1 obs pt
1154 ! mgIntrp -internal- max. interpolated MG value on mesh for all obs pts
1155 ! lgIntrp -internal- max. interpolated LG value on mesh for all obs pts
1156 ! MGthresh -internal- maximum allowable land fraction for obs to be kept
1157 ! LGthresh -internal- maximum allowable ice fraction for obs to be kept
1158 !--------------------------------------------------------------------
1159 !
1160 implicit none
1161
1162 ! Arguments:
1163 integer, intent(in) :: numObsToProcess ! Number of obs points to process
1164 real, intent(in) :: obsLatitude(:) ! Observation latitudes
1165 real, intent(in) :: obsLongitude(:) ! Observation longitudes
1166 integer, intent(inout) :: landSeaQualifier(:) ! Land/sea indicator (0=land; 1=ocean)
1167 integer, allocatable, intent(out) :: terrainType(:) ! Terrain type qualifier
1168 logical, allocatable, intent(out) :: waterObs(:) ! Open water identifier for each obs
1169
1170 ! Locals:
1171 integer, parameter :: mxLat=5
1172 integer, parameter :: mxLon=5
1173 real, parameter :: LGthresh=0.01
1174 real, parameter :: MGthresh=0.01
1175 real, parameter :: pi=3.141592654
1176 real, parameter :: rKmPerDeg=111.195
1177 real, parameter :: rLatKm=40.0
1178 real, parameter :: rLonKm=40.0
1179 character(len=12) :: etikxx
1180 character(len=1) :: grtyp
1181 character(len=1) :: grtyplg
1182 character(len=4) :: nomvxx
1183 character(len=2) :: typxx
1184 integer :: boxPointIndex
1185 integer :: ezQkDef
1186 integer :: ezSetOpt
1187 integer, save :: gdId
1188 integer, save :: gdIdlg
1189 integer :: gdllsval
1190 integer :: idum1, idum2, idum3, idum4
1191 integer :: idum5, idum6, idum7, idum8
1192 integer :: idum9, idum10, idum11, idum12
1193 integer :: idum13, idum14, idum15
1194 integer :: idum16, idum17, idum18
1195 integer :: ier
1196 integer :: ig1, ig2, ig3, ig4
1197 integer :: ig1lg, ig2lg, ig3lg, ig4lg
1198 integer :: iMask
1199 integer :: irec
1200 integer :: iUnGeo
1201 integer :: latIndex
1202 integer :: lonIndex
1203 integer :: ni, nj, nk
1204 integer :: nilg, njlg
1205 integer :: nLat
1206 integer :: nLon
1207 integer :: obsIndex
1208 logical, save :: firstCall=.true.
1209 real, allocatable, save :: lg(:)
1210 real, allocatable, save :: mg(:)
1211 real, allocatable :: latMesh(:)
1212 real, allocatable :: lgIntOb(:)
1213 real, allocatable :: lgIntrp(:)
1214 real, allocatable, save :: lonMesh(:)
1215 real, allocatable :: mgIntOb(:)
1216 real, allocatable :: mgIntrp(:)
1217 real, allocatable :: obsLatBox(:,:)
1218 real, allocatable :: obsLonBox(:,:)
1219 real :: dLat
1220 real :: dLon
1221 real :: rLatIndex
1222 real :: rLonIndex
1223 real :: xLat
1224 real :: xLatRad
1225 real :: xLon
1226 ! External functions
1227 integer, external :: fclos
1228 integer, external :: fnom
1229 integer, external :: fstfrm
1230 integer, external :: fstinf
1231 integer, external :: fstlir
1232 integer, external :: fstouv
1233 integer, external :: fstprm
1234
1235 ! Allocate space for arrays holding values on mesh grid pts.
1236 call utl_reAllocate(latMesh, mxLat*mxLon)
1237 call utl_reAllocate(lonMesh, mxLat*mxLon)
1238 call utl_reAllocate(mgIntOb, mxLat*mxLon)
1239 call utl_reAllocate(lgIntOb, mxLat*mxLon)
1240 call utl_reAllocate(obsLatBox, mxLat*mxLon, numObsToProcess)
1241 call utl_reAllocate(obsLonBox, mxLat*mxLon, numObsToProcess)
1242 call utl_reAllocate(terrainType, numObsToProcess)
1243 call utl_reAllocate(waterObs, numObsToProcess)
1244
1245 if (firstCall) then
1246
1247 firstCall = .false.
1248
1249 ! Open FST file.
1250 iUnGeo = 0
1251 ier = fnom( iUnGeo,fileMgLg,'STD+RND+R/O',0 )
1252 ier = fstouv( iUnGeo,'RND' )
1253
1254 ! Read MG field.
1255 irec = fstinf(iUnGeo,ni,nj,nk,-1,' ',-1,-1,-1,' ' ,'MG')
1256 if ( irec < 0 ) then
1257 call utl_abort('land_ice_mask_ssmis: The MG field is MISSING')
1258 end if
1259
1260 call utl_reAllocate(mg, ni*nj)
1261
1262 ier = fstlir(mg,iUnGeo,ni,nj,nk,-1,' ',-1,-1,-1,' ','MG')
1263
1264 ier = fstprm(irec,idum1,idum2,idum3,idum4,idum5,idum6,idum7,idum8, &
1265 idum9,idum10,idum11,typxx,nomvxx,etikxx,grtyp,ig1,ig2, &
1266 ig3,ig4,idum12,idum13,idum14,idum15,idum16,idum17, &
1267 idum18)
1268
1269 ! Read LG field.
1270
1271 irec = fstinf(iUnGeo,nilg,njlg,nk,-1,' ',-1,-1,-1,' ' ,'LG')
1272 if ( irec < 0 ) then
1273 call utl_abort('land_ice_mask_ssmis: The LG field is MISSING ')
1274 end if
1275 call utl_reAllocate(lg, nilg*njlg)
1276 ier = fstlir(lg,iUnGeo,nilg,njlg,nk,-1,' ',-1,-1,-1,' ','LG')
1277
1278 ier = fstprm(irec,idum1,idum2,idum3,idum4,idum5,idum6,idum7,idum8, &
1279 & idum9,idum10,idum11,typxx,nomvxx,etikxx,grtyplg,ig1lg,ig2lg, &
1280 & ig3lg,ig4lg,idum12,idum13,idum14,idum15,idum16,idum17, &
1281 & idum18)
1282
1283 gdId = ezQkDef(ni,nj,grtyp,ig1,ig2,ig3,ig4,iUnGeo)
1284 gdIdlg = ezQkDef(nilg,njlg,grtyplg,ig1lg,ig2lg,ig3lg,ig4lg,iUnGeo)
1285
1286 ier = fstfrm(iUnGeo)
1287 ier = fclos(iUnGeo)
1288
1289 end if ! firstCall
1290
1291
1292 ! For each obs pt, define a grid of artificial pts surrounding it.
1293
1294 nLat = (mxLat-1)/2
1295 nLon = (mxLon-1)/2
1296
1297 dLat = rLatKm / rKmPerDeg
1298 do obsIndex = 1, numObsToProcess
1299 boxPointIndex = 0
1300
1301 do latIndex = -nLat, nLat
1302 rlatIndex = float(latIndex)
1303 xLat = obsLatitude(obsIndex) + rLatIndex*dLat
1304 xLat = max( -90.0, min(90.0,xLat) )
1305 xLatRad = xLat*pi/180.0
1306
1307 do lonIndex = -nLon, nLon
1308 dLon = rLonKm / ( rKmPerDeg*cos(xLatRad) )
1309 rLonIndex = float(lonIndex)
1310 boxPointIndex = boxPointIndex + 1
1311 xLon = obsLongitude(obsIndex) + rLonIndex*dLon
1312 if ( xLon < -180. ) xLon = xLon + 360.
1313 if ( xLon > 180. ) xLon = xLon - 360.
1314 if ( xLon < 0. ) xLon = xLon + 360.
1315 obsLatBox(boxPointIndex,obsIndex) = xLat
1316 obsLonBox(boxPointIndex,obsIndex) = xLon
1317 end do
1318
1319 end do
1320 end do
1321
1322
1323 ! Interpolate values from MG and LG field to grid pts of mesh centred over each obs pt.
1324 ! Determine for each obs pt, the max interpolated MG and LG value within the box
1325 ! surrounding it.
1326
1327 ier = ezSetOpt('INTERP_DEGREE','LINEAR')
1328
1329 call utl_reAllocate(mgIntrp, numObsToProcess)
1330 call utl_reAllocate(lgIntrp, numObsToProcess)
1331
1332 mgIntrp(:) = 0.0
1333 lgIntrp(:) = 0.0
1334 do obsIndex = 1, numObsToProcess
1335
1336 latMesh = obsLatBox(:,obsIndex)
1337 lonMesh = obsLonBox(:,obsIndex)
1338
1339 ier = gdllsval(gdId,mgIntOb,mg,latMesh,lonMesh,mxLat*mxLon)
1340 ier = gdllsval(gdIdlg,lgIntOb,lg,latMesh,lonMesh,mxLat*mxLon)
1341
1342 mgIntrp(obsIndex) = maxval(mgIntOb(:))
1343 lgIntrp(obsIndex) = maxval(lgIntOb(:))
1344
1345 end do
1346
1347 ! Initialize all obs as being over land and free of ice or snow.
1348 ! Determine which obs are over open water.
1349
1350 waterObs(:) = .false. ! not over open water
1351 terrainType(:) = -1 ! no ice or snow
1352 HEADER: do obsIndex = 1, numObsToProcess
1353
1354 ! Determine for each obs pt, the value of iMask from the 0.1 deg binary file.
1355 ! - open land/sea: iMask = -1 (missing)
1356 ! - ice or snow: iMask = 0 (from LG field --> binary ice mask)
1357 ! Define the terrain-type qualifier terrainType for each point based on the ice mask values.
1358 ! terrainType = -1 not defined (open water or snow free land)
1359 ! 0 sea-ice (landSeaQualifier = 1)
1360 ! 1 snow-covered land (landSeaQualifier = 0)
1361
1362 if (lgintrp(obsIndex) < LGthresh ) then
1363 iMask = -1
1364 else
1365 iMask = 0
1366 end if
1367
1368 if ( iMask == 0 ) then ! if ice or snow
1369 terrainType(obsIndex) = 1 - landSeaQualifier(obsIndex)
1370 end if
1371
1372 ! If iMask is -1 (no ice/snow), and this is consistent with the model ice
1373 ! LG value (ie. < LGthresh), and the max MG value indicates ocean (ie.
1374 ! < MGthresh), then this is a WATEROBS point.
1375
1376 if ( iMask == -1 .and. lgintrp(obsIndex) < LGthresh .and. mgintrp(obsIndex) < MGthresh ) then
1377 !!!! TO RUN WITHOUT BINARY ICE MASK CHECK (I.E. RELY ON LG ONLY TO REMOVE ICE PTS):
1378 !if ( lgintrp(obsIndex) < LGthresh .and. mgintrp(obsIndex) < MGthresh ) then
1379 waterObs(obsIndex) = .true.
1380 end if
1381
1382 ! Modify land/sea quailifier if not consistent with waterObs (applies to remote small
1383 ! islands that should be treated as sea points (for RTTOV)):
1384 ! -- if waterObs=true, land/sea qualifier should be "sea" value (1)
1385
1386 if ( waterObs(obsIndex) .and. (landSeaQualifier(obsIndex) == 0) ) landSeaQualifier(obsIndex) = 1
1387
1388 end do HEADER
1389
1390 end subroutine land_ice_mask_ssmis
1391
1392 !--------------------------------------------------------------------------
1393 ! wentz_sfctype_ssmis
1394 !--------------------------------------------------------------------------
1395 subroutine wentz_sfctype_ssmis(numObsToProcess, obsLatitude, obsLongitude, landSeaQualifier)
1396 !:Purpose: Determine for each observation point the wentz surface value
1397 ! from the FST file wentz_surf.std.
1398
1399 ! This file has a resolution of 0.25 deg, and it discriminates
1400 ! between the following surface types:
1401 ! a) land (wentz value = 0)
1402 ! b) ice (wentz value = 4)
1403 ! c) sea (wentz value = 5)
1404 ! d) coast (wentz value = 6)
1405 ! This information is used to define the land/sea qualifier
1406 ! (element 008012) for each obs pt, which is needed by 3D-Var.
1407 ! Also, the land/sea qualifier is used later to flag obs pts
1408 ! over land and coast.
1409 !
1410 !--------------------------------------------------------------------
1411 ! Variable Definitions
1412 ! --------------------
1413 ! numObsToProcess - input - number of input obs pts in record
1414 ! obsLatitude - input - array of size ssbg_maxObsNum holding lat values for obs pts
1415 ! in record plus undefined pts
1416 ! obsLongitude - input - array of size ssbg_maxObsNum holding lon values for obs pts
1417 ! in record plus undefined pts
1418 ! landSeaQualifier - output - array holding land/sea qualifier values for all obs pts of record
1419 ! xLat -internal- array of size numObsToProcess holding lat values for obs pts in record
1420 ! xLon -internal- array of size numObsToProcess holding lon values for obs pts in record
1421 ! lm -internal- array of size 1440x720 holding gridded wentz surface values
1422 ! wenTyp -internal- array of size numObsToProcess holding wentz surface values interpolated to obs pts
1423 !--------------------------------------------------------------------
1424 !
1425 implicit none
1426
1427 ! Arguments:
1428 integer, intent(in) :: numObsToProcess ! Number of obs points to process
1429 real, intent(in) :: obsLatitude(:) ! Observation latitudes
1430 real, intent(in) :: obsLongitude(:) ! Observation longitudes
1431 integer, allocatable, intent(out) :: landSeaQualifier(:) ! Land/sea indicator (0=land; 1=ocean)
1432
1433 ! Locals:
1434 character(len=12) :: etikxx
1435 character(len=1) :: grtyp
1436 character(len=4) :: nomvxx
1437 character(len=2) :: typxx
1438 integer :: ezQkDef
1439 integer :: ezSetOpt
1440 integer :: gdId
1441 integer :: gdllsval
1442 integer :: idum1, idum2, idum3, idum4
1443 integer :: idum5, idum6, idum7, idum8
1444 integer :: idum9, idum10, idum11, idum12
1445 integer :: idum13, idum14, idum15
1446 integer :: idum16, idum17, idum18
1447 integer :: ier
1448 integer :: ig1, ig2, ig3, ig4
1449 integer :: irec
1450 integer :: iUnIn
1451 integer :: ni, nj, nk
1452 integer :: obsIndex
1453 real, allocatable :: xLat(:)
1454 real, allocatable :: xLon(:)
1455 real, allocatable :: wenTyp(:)
1456 real, allocatable :: lm(:)
1457 ! External functions
1458 integer, external :: fclos
1459 integer, external :: fnom
1460 integer, external :: fstfrm
1461 integer, external :: fstinf
1462 integer, external :: fstlir
1463 integer, external :: fstouv
1464 integer, external :: fstprm
1465
1466 ! Open Wentz surface field if first call
1467 iUnIn = 0
1468 ier = fnom( iUnIn,fileWentz,'STD+RND+R/O',0 )
1469 ier = fstouv( iUnIn,'RND' )
1470
1471 irec = fstinf(iUnIn,ni,nj,nk,-1,' ',0,0,0,' ','LM')
1472 if ( irec < 0 ) then
1473 call utl_abort('wentz_sfctype_ssmis: The LM field is MISSING ')
1474 else
1475 call utl_reAllocate( lm, ni*nj )
1476 ier = fstlir(lm,iUnIn,ni,nj,nk,-1,' ',-1,-1,-1,' ','LM')
1477 end if
1478
1479 ier = fstprm(irec,idum1,idum2,idum3,idum4,idum5,idum6,idum7,idum8, &
1480 & idum9,idum10,idum11,typxx,nomvxx,etikxx,grtyp,ig1,ig2, &
1481 & ig3,ig4,idum12,idum13,idum14,idum15,idum16,idum17, &
1482 & idum18)
1483
1484 ! Re-define the grid type.
1485 ! - L-grid: IG's SHOULD define grid spacing over a REGIONAL
1486 ! domain. Here, IG's are =0 though.
1487
1488 ! - A-grid: is for global grid. See RPN documentation on
1489 ! grid-types for proper IG values.
1490 grtyp ='A'
1491
1492 ! Transfer data from obsLatitude,obsLongitude to xLat,xLon.
1493 ! Ensure proper range for lon values.
1494
1495 call utl_reAllocate( xLat, numObsToProcess)
1496 call utl_reAllocate( xLon, numObsToProcess)
1497 call utl_reAllocate( wenTyp, numObsToProcess)
1498 call utl_reAllocate( landSeaQualifier, numObsToProcess)
1499
1500 landSeaQualifier(:) = 0
1501 xLat(:) = obsLatitude(1:numObsToProcess)
1502 xLon(:) = obsLongitude(1:numObsToProcess)
1503 do obsIndex = 1, numObsToProcess
1504 if ( xLon(obsIndex) < 0. ) xLon(obsIndex) = xLon(obsIndex) + 360.
1505 end do
1506
1507 ! Interpolate values from LM field to all obs pts of record.
1508 ier = ezSetOpt('INTERP_DEGREE','VOISIN')
1509 gdId = ezQkDef(ni,nj,grTyp,ig1,ig2,ig3,ig4,iUnIn)
1510 ier = gdllsval(gdId,wenTyp,lm,xLat,xLon,numObsToProcess)
1511
1512 ! Define the land/sea qualifier for each point based on wentz surface values.
1513 do obsIndex = 1, numObsToProcess
1514 if ( wenTyp(obsIndex) == 0. .or. wenTyp(obsIndex) == 6. ) then
1515 ! wentz = land/coast --> land
1516 landSeaQualifier(obsIndex) = 0
1517 else if ( wenTyp(obsIndex) == 4. .or. wenTyp(obsIndex) == 5. ) then
1518 ! wentz = sea/sea-ice --> sea
1519 landSeaQualifier(obsIndex) = 1
1520 else
1521 call utl_abort('wentz_sfctype_ssmis: Unexpected Wentz value ')
1522 end if
1523 end do
1524 ier = fstfrm(iUnIn)
1525 ier = fclos(iUnIn)
1526
1527 end subroutine wentz_sfctype_ssmis
1528
1529 !--------------------------------------------------------------------------
1530 ! ssbg_computeSsmisSurfaceType
1531 !--------------------------------------------------------------------------
1532 subroutine ssbg_computeSsmisSurfaceType(obsSpaceData)
1533 !
1534 !:Purpose: Compute surface type element and update obsSpaceData.
1535 !
1536 implicit none
1537
1538 ! Arguments:
1539 type(struct_obs), intent(inout) :: obsSpaceData ! ObsSpaceData object
1540
1541 ! Locals:
1542 real, parameter :: satAzimuthAngle = 210.34
1543 real, parameter :: satZenithAngle = 53.1
1544 integer, allocatable :: landSeaQualifier(:)
1545 integer :: codtyp
1546 integer :: headerIndex
1547 logical :: ssmisDataPresent
1548 real :: obsLatitude(1)
1549 real :: obsLongitude(1)
1550
1551 write(*,*) 'ssbg_computeSsmisSurfaceType: Starting'
1552
1553 ssmisDataPresent=.false.
1554 call obs_set_current_header_list(obsSpaceData,'TO')
1555
1556 HEADER0: do
1557 headerIndex = obs_getHeaderIndex(obsSpaceData)
1558 if (headerIndex < 0) exit HEADER0
1559 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
1560 if ( tvs_isIdBurpInst(codtyp,'ssmis') ) then
1561 ssmisDataPresent = .true.
1562 exit HEADER0
1563 end if
1564 end do HEADER0
1565
1566 if ( .not. ssmisDataPresent ) then
1567 write(*,*) 'WARNING: WILL NOT RUN ssbg_computeSsmisSurfaceType since no SSMIS DATA is found'
1568 return
1569 end if
1570
1571 call obs_set_current_header_list(obsSpaceData,'TO')
1572 HEADER1: do
1573 headerIndex = obs_getHeaderIndex(obsSpaceData)
1574 if (headerIndex < 0) exit HEADER1
1575 obsLatitude(1) = obs_headElem_r( obsSpaceData, OBS_LAT, headerIndex )
1576 obsLongitude(1) = obs_headElem_r( obsSpaceData, OBS_LON, headerIndex )
1577 obsLongitude(1) = obsLongitude(1)*MPC_DEGREES_PER_RADIAN_R8
1578 if( obsLongitude(1) > 180. ) obsLongitude(1) = obsLongitude(1) - 360.
1579 obsLatitude(1) = obsLatitude(1) *MPC_DEGREES_PER_RADIAN_R8
1580 call wentz_sfctype_ssmis(1, obsLatitude, obsLongitude, landSeaQualifier)
1581 call obs_headSet_i(obsSpaceData, OBS_STYP, headerIndex, landSeaQualifier(1))
1582 call obs_headSet_r(obsSpaceData, OBS_SZA, headerIndex,satZenithAngle)
1583 call obs_headSet_r(obsSpaceData, OBS_AZA, headerIndex, satAzimuthAngle)
1584 end do HEADER1
1585
1586 write(*,*) 'ssbg_computeSsmisSurfaceType: Finished'
1587
1588 end subroutine ssbg_computeSsmisSurfaceType
1589
1590 !--------------------------------------------------------------------------
1591 ! ssbg_grossValueCheck
1592 !--------------------------------------------------------------------------
1593 subroutine ssbg_grossValueCheck(numObsToProcess, obsTb, obsTbMin, obsTbMax, grossRej)
1594 !
1595 !:Purpose: Check obsTb for values that are missing or outside physical limits.
1596 !
1597 implicit none
1598
1599 ! Arguments:
1600 integer, intent(in) :: numObsToProcess ! Number of obs points to process
1601 real, intent(in) :: obsTb(:) ! Brightness temperature of observations
1602 real, intent(in) :: obsTbMin ! Min(obsTb) threshold for rejection
1603 real, intent(in) :: obsTbMax ! Max(obsTb) threshold for rejection
1604 logical, allocatable, intent(out) :: grossRej(:) ! Logical array of obs with gross error (obs to reject)
1605
1606 ! Locals:
1607 integer :: hiIndex
1608 integer :: loIndex
1609 integer :: obsIndex
1610
1611 call utl_reAllocate(grossRej, numObsToProcess)
1612
1613 grossRej(1:numObsToProcess) = .true.
1614 loIndex = 1
1615 do obsIndex = 1, numObsToProcess
1616
1617 hiIndex = obsIndex*ssbg_maxNumChan
1618 if ( all( obsTb(loIndex:hiIndex) > obsTbMin ) .and. all( obsTb(loIndex:hiIndex) < obsTbMax ) ) then
1619 grossRej(obsIndex) = .false.
1620 end if
1621 loIndex = hiIndex + 1
1622
1623 end do
1624
1625 end subroutine ssbg_grossValueCheck
1626
1627 !--------------------------------------------------------------------------
1628 ! ssbg_satqcSsmis
1629 !--------------------------------------------------------------------------
1630 subroutine ssbg_satqcSsmis(obsSpaceData, headerIndex, obsToReject)
1631 !
1632 !:Purpose: This program is applied as a first stage of processing to
1633 ! SSMIS data after it is received from UK MetOffice and
1634 ! organized into boxes by a program of Jose Garcia. The
1635 ! processing applied in this program includes:
1636 ! -- interpolate Wentz surface land mask to each obs pt
1637 ! (nearest neighbour) to define land/sea qualifier (008012)
1638 ! -- interpolate binary ice mask to each obs pt (nearest
1639 ! neighbour) to define terrain-type element (013039) where
1640 ! 0 = sea ice and 1 = snow-covered land
1641 ! -- interpolate model MG and LG fields to a grid surrounding each obs
1642 ! pt to identify obs that are over open water, far from coast/ice
1643 ! -- identify those obs for which the UKMO rain marker
1644 ! is ON (ie. 020029 = 1) indicating poor quality
1645 ! -- apply a cloud filter to identify those obs in cloudy regions;
1646 ! write CLW and IWV (over ocean) to output BURP file
1647 ! -- re-write data to output BURP file while modifying flags
1648 ! for those obs which are not over open water, or have been
1649 ! identified in rain/cloud areas, or are of poor quality
1650 ! -- define satellite zenith angle element (007024) and add
1651 ! this and land/sea qualifier and terrain-type elements
1652 ! to the output file
1653 !
1654 implicit none
1655
1656 ! Arguments:
1657 type(struct_obs), intent(inout) :: obsSpaceData ! ObsSpaceData object
1658 integer, intent(in) :: headerIndex ! Current header index
1659 logical, allocatable, intent(out) :: obsToReject(:) ! Observations that will be rejected
1660
1661 ! Locals:
1662 ! arrays to get from obsspacedata
1663 character(len=9) :: burpFileSatId
1664 integer, allocatable :: landSeaQualifier(:)
1665 integer, allocatable :: terrainType(:)
1666 integer, allocatable :: ukRainObs(:)
1667 real, allocatable :: obsLatitude(:)
1668 real, allocatable :: obsLongitude(:)
1669 real, allocatable :: obsTb(:)
1670 real, allocatable :: rclw(:)
1671 real, allocatable :: riwv(:)
1672 real, allocatable :: satZenithAngle(:)
1673 real, allocatable :: scatW(:)
1674 integer, allocatable :: ier(:)
1675 integer :: actualNumChannel
1676 integer :: bodyIndex
1677 integer :: bodyIndexBeg
1678 integer :: bodyIndexEnd
1679 integer :: codtyp
1680 integer :: currentChannelNumber
1681 integer :: headerCompt
1682 integer :: instr
1683 integer :: instrId
1684 integer :: obsIndex
1685 integer :: iSat
1686 integer :: iSatId
1687 integer, save :: numLandObs
1688 integer, save :: numUkBadObs
1689 integer, save :: numGrossObs
1690 integer, save :: numCloudyObs
1691 integer, save :: numBadIWVObs
1692 integer, save :: numPrecipObs
1693 integer, save :: numTotFilteredObs
1694 integer, save :: numLandScatObs
1695 integer, save :: numSeaScatObs
1696 integer, save :: numDryIndexObs
1697 integer, save :: numSeaIceObs
1698 integer, save :: numScatPrecipObs
1699 integer, save :: numCloudsinObs
1700 integer, save :: numWaterObs
1701 integer, save :: numObsF16
1702 integer, save :: numObsF17
1703 integer, save :: numObsF18
1704 integer :: numObsToProcess
1705 integer :: platf
1706 integer :: platfId
1707 integer :: refPosObs
1708 integer :: sensorIndex
1709 logical, allocatable :: cloudObs(:)
1710 logical, allocatable :: grossRej(:)
1711 logical, allocatable :: iwvReject(:)
1712 logical, allocatable :: precipObs(:)
1713 logical, allocatable :: rainDetectionUKMethod(:)
1714 logical, allocatable :: waterObs(:)
1715 logical, save :: ifFirstCall = .true.
1716 logical :: sensorIndexFound
1717 logical :: ssmisDataPresent
1718 real, allocatable :: amsubDrynessIndex(:)
1719 real, allocatable :: scatL(:)
1720 real, allocatable :: ztb91(:)
1721 real, allocatable :: ztb150(:)
1722 real, allocatable :: ztb_amsub3(:)
1723 real, allocatable :: ztb_amsub5(:)
1724
1725 ! Check if its ssmis data:
1726 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
1727 ssmisDataPresent = tvs_isIdBurpInst(codtyp,'ssmis')
1728
1729 if ( .not. ssmisDataPresent ) then
1730 write(*,*) 'WARNING: WILL NOT RUN ssbg_satqcSsmis since no SSMIS DATA is found'
1731 return
1732 end if
1733
1734 ! find tvs_sensor index corresponding to current obs
1735
1736 platf = obs_headElem_i( obsSpaceData, OBS_SAT, headerIndex )
1737 instr = obs_headElem_i( obsSpaceData, OBS_INS, headerIndex )
1738
1739 call tvs_mapSat( platf, platfId, iSat )
1740 call tvs_mapInstrum( instr, instrId )
1741
1742 sensorIndexFound = .false.
1743 HEADER0: do sensorIndex = 1, tvs_nsensors
1744 if ( platfId == tvs_platforms(sensorIndex) .and. &
1745 iSat == tvs_satellites(sensorIndex) .and. &
1746 instrId == tvs_instruments(sensorIndex) ) then
1747 sensorIndexFound = .true.
1748 exit HEADER0
1749 end if
1750 end do HEADER0
1751 if ( .not. sensorIndexFound ) call utl_abort('ssbg_satqcSsmis: sensor Index not found')
1752
1753 ! find actual Number of channels
1754 actualNumChannel = tvs_coefs(sensorIndex)%coef%fmv_ori_nchn
1755
1756 headerCompt = 1
1757 numObsToProcess = 1
1758
1759 ! Allocate intent out arrays
1760
1761 call utl_reAllocate(obsToReject, numObsToProcess*actualNumChannel)
1762
1763 ! Allocate Fortran working arrays
1764
1765 call utl_reAllocate(amsubDrynessIndex, numObsToProcess)
1766 call utl_reAllocate(cloudObs, numObsToProcess)
1767 call utl_reAllocate(grossRej, numObsToProcess)
1768 call utl_reAllocate(ier, numObsToProcess)
1769 call utl_reAllocate(iwvReject, numObsToProcess)
1770 call utl_reAllocate(precipObs, numObsToProcess)
1771 call utl_reAllocate(rainDetectionUKMethod, numObsToProcess)
1772 call utl_reAllocate(scatL, numObsToProcess)
1773 call utl_reAllocate(waterObs, numObsToProcess)
1774 call utl_reAllocate(ztb91, numObsToProcess)
1775 call utl_reAllocate(ztb150, numObsToProcess)
1776 call utl_reAllocate(ztb_amsub3, numObsToProcess)
1777 call utl_reAllocate(ztb_amsub5, numObsToProcess)
1778
1779 ! ELEMENTS FROM OBSSPACEDATA
1780 ! Allocate Header elements
1781 call utl_reAllocate(landSeaQualifier, numObsToProcess)
1782 call utl_reAllocate(obsLatitude, numObsToProcess)
1783 call utl_reAllocate(obsLongitude, numObsToProcess)
1784 call utl_reAllocate(rclw, numObsToProcess)
1785 call utl_reAllocate(riwv, numObsToProcess)
1786 call utl_reAllocate(satZenithAngle, numObsToProcess)
1787 call utl_reAllocate(scatW, numObsToProcess)
1788 call utl_reAllocate(terrainType, numObsToProcess)
1789 call utl_reAllocate(ukRainObs, numObsToProcess)
1790 ! Allocate Body elements
1791 call utl_reAllocate(obsTb, numObsToProcess*actualNumChannel)
1792 !initialization
1793 obsTb(:) = ssbg_realMissing
1794 riwv(:) = ssbg_realMissing
1795
1796 ! Lecture dans obsspacedata
1797
1798 burpFileSatId = obs_elem_c ( obsSpaceData, 'STID' , headerIndex )
1799 rclw(headerCompt) = obs_headElem_r( obsSpaceData, OBS_CLWO, headerIndex )
1800 ukRainObs(headerCompt) = obs_headElem_i( obsSpaceData, OBS_RAIN, headerIndex )
1801 landSeaQualifier(headerCompt) = obs_headElem_i( obsSpaceData, OBS_STYP, headerIndex )
1802 terrainType(headerCompt) = obs_headElem_i( obsSpaceData, OBS_TTYP, headerIndex )
1803 ! If terrain type is missing, set it to -1 for the QC programs
1804 if (terrainType(headerCompt) == 99) terrainType(headerCompt) = -1
1805 obsLatitude(headerCompt) = obs_headElem_r( obsSpaceData, OBS_LAT, headerIndex )
1806 obsLongitude(headerCompt) = obs_headElem_r( obsSpaceData, OBS_LON, headerIndex )
1807 satZenithAngle(headerCompt) = obs_headElem_r( obsSpaceData, OBS_SZA, headerIndex )
1808 ! Convert lat/lon to degrees
1809 obsLatitude(headerCompt) = obsLatitude(headerCompt) *MPC_DEGREES_PER_RADIAN_R8
1810 obsLongitude(headerCompt) = obsLongitude(headerCompt)*MPC_DEGREES_PER_RADIAN_R8
1811 if( obsLongitude(headerCompt) > 180. ) obsLongitude(headerCompt) = obsLongitude(headerCompt) - 360.
1812
1813 ! To read body elements
1814 bodyIndexBeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
1815 bodyIndexEnd = bodyIndexBeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) -1
1816
1817 BODY: do bodyIndex = bodyIndexbeg, bodyIndexEnd
1818 currentChannelNumber = nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex ))-tvs_channelOffset(sensorIndex)
1819 if (obs_bodyElem_r( obsSpaceData, OBS_BCOR, bodyIndex ) /= ssbg_rmisg) then
1820 obsTb(currentChannelNumber) = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex ) - obs_bodyElem_r( obsSpaceData, OBS_BCOR, bodyIndex )
1821 else
1822 obsTb(currentChannelNumber) = obs_bodyElem_r( obsSpaceData, OBS_VAR, bodyIndex )
1823 end if
1824 end do BODY
1825
1826 ! initialization
1827 if (ifFirstCall) then
1828 ifFirstCall = .false.
1829 numBadIWVObs = 0
1830 numCloudyObs = 0
1831 numDryIndexObs = 0
1832 numGrossObs = 0
1833 numLandObs = 0
1834 numLandScatObs = 0
1835 numObsF16 = 0
1836 numObsF17 = 0
1837 numObsF18 = 0
1838 numPrecipObs = 0
1839 numScatPrecipObs = 0
1840 numSeaIceObs = 0
1841 numSeaScatObs = 0
1842 numTotFilteredObs = 0
1843 numUkBadObs = 0
1844 numCloudsinObs = 0
1845 numWaterObs = 0
1846 end if
1847 waterObs(:) = .false.
1848
1849 ! Record the total number of obs pts read for each satellite.
1850 ! Set the satellite ID number.
1851 select case (burpFileSatId(2:7) )
1852 case ( 'DMSP16' )
1853 numObsF16 = numObsF16 + numObsToProcess
1854 iSatId = 1
1855 case ( 'DMSP17' )
1856 numObsF17 = numObsF17 + numObsToProcess
1857 iSatId = 2
1858 case ( 'DMSP18' )
1859 numObsF18 = numObsF18 + numObsToProcess
1860 iSatId = 3
1861 end select
1862
1863 !--------------------------------------------------------------------
1864 ! Determine which obs pts are over open water (i.e NOT near coasts or
1865 ! over/near land/ice.
1866 ! - using model MG field, plus 0.1 deg ice mask and model LG field
1867 ! Also, define the terrain-type qualifier values (itt).
1868 ! Also, ensure land/sea qualifier (ilq) = 1 (sea) for points identified
1869 ! as open water by waterobs array.
1870 !--------------------------------------------------------------------
1871
1872 call land_ice_mask_ssmis(numObsToProcess, obsLatitude, obsLongitude, landSeaQualifier, &
1873 terrainType, waterObs)
1874
1875 !--------------------------------------------------------------------
1876 ! Determine which obs pts have been flagged for rain (using UKMO method of identifying
1877 ! bad quality SSMIS data for DMSP16). First, initialize all obs as good.
1878 !--------------------------------------------------------------------
1879
1880 rainDetectionUKMethod(:) = .false.
1881 where ( ukRainObs == 1 ) rainDetectionUKMethod = .true.
1882
1883 !--------------------------------------------------------------------
1884 ! Check for values of TB that are missing or outside physical limits.
1885 ! **NOTE: REJECT ALL CHANNELS IF ONE IS FOUND TO BE BAD.
1886 !--------------------------------------------------------------------
1887
1888 grossRej(:) = .false.
1889 call ssbg_grossValueCheck(numObsToProcess, obsTb, 50., 400., grossRej)
1890
1891 !--------------------------------------------------------------------
1892 ! Apply a CLW regression technique to determine cloudy obs pts.
1893 ! Apply a IWV regression technique to determine IWV obs out of bounds.
1894 ! Detect Sea_Ice percent from Tb(Ta) and set waterobs() flag to .FALSE. if
1895 ! amount > 70%
1896 ! To begin, assume that all obs are good.
1897 !--------------------------------------------------------------------
1898
1899 cloudobs(:) = .false.
1900 iwvReject(:) = .false.
1901 precipobs(:) = .false.
1902
1903 ! --------------------------------------------------------------------
1904
1905 call cld_filter_fweng(numObsToProcess, obsTb, algOption, waterObs, grossRej, &
1906 cloudObs, iwvReject, precipObs, rclw, riwv, iSatId, &
1907 obsLatitude, numSeaIceObs)
1908 !
1909 ! NOTE: rclw, riwv missing value is zmisg=9.9e09
1910 ! --> if ( (clw > clw_thresh) .or. precipobs(obsIndex) ) cloudobs(obsIndex) = .true.
1911 ! --> if ( clw == zmisg ) cloudobs(obsIndex) = .true.
1912 !
1913 ! Extract Tb for channels 18 (91H GHz) and 8 (150H GHz) for Bennartz SI
1914 ! Extract Tb for channels 11 (AMSU-B 3) and 9 (AMSU-B 5) for Dryness Index (DI)
1915
1916 refPosObs = 1
1917 do obsIndex = 1, numObsToProcess
1918 ztb91(obsIndex) = obsTb(refPosObs+17)
1919 ztb150(obsIndex) = obsTb(refPosObs+7)
1920 ztb_amsub3(obsIndex) = obsTb(refPosObs+10)
1921 ztb_amsub5(obsIndex) = obsTb(refPosObs+8)
1922 refPosObs = obsIndex*ssbg_maxNumChan + 1
1923 end do
1924
1925 !---------------------------------------------------------------------------------
1926 ! Compute the Bennartz scattering index (SI) for each point from channels 8 and 18
1927 ! (detects precipitation affected AMSU-B-like channel radiances)
1928 ! -- using constant satellite zenith angle = 53.1
1929 ! -- SSMIS channel 18 (91.655 GHz) is used for AMSU-B like channel 1 (89 GHz)
1930 ! -- SSMIS channel 8 is used for AMSU-B like channel 2 (both at 150 GHz)
1931 !---------------------------------------------------------------------------------
1932
1933 call bennartz(ier, numObsToProcess, ztb91, ztb150, satZenithAngle, landSeaQualifier, &
1934 scatL, scatW)
1935
1936 !--------------------------------------------------------------------
1937 ! Compute the AMSU-B Dryness Index for all points
1938 !--------------------------------------------------------------------
1939
1940 where ( (ztb_amsub3 /= ssbg_realMissing) .and. (ztb_amsub5 /= ssbg_realMissing) )
1941 amsubDrynessIndex = ztb_amsub3 - ztb_amsub5
1942 elsewhere
1943 amsubDrynessIndex = ssbg_realMissing
1944 end where
1945
1946 !--------------------------------------------------------------------
1947 ! Review all the checks previously made to determine which obs are to be accepted
1948 ! for assimilation and which are to be flagged for exclusion. First, initialize
1949 ! all obs for inclusion.
1950 ! ukBadObs() = .true. if UKMetO QC sets "rain" flag: solar-intrusion or other anomaly
1951 ! grossRej() = .true. if any channel had a gross error at the point
1952 ! cloudObs() = .true. if CLW > 0.01 kg/m**2 or precip (over water)
1953 ! precipObs() = .true. if precip. detected (CLW=missing)
1954 ! waterObs() = .true. if open water point (far from land or sea ice)
1955 ! iwvReject() = .true. if IWV > 80 kg/m**2
1956 !--------------------------------------------------------------------
1957
1958 obsToReject(:) = .false.
1959 HEADER1: do obsIndex = 1, numObsToProcess
1960
1961 ! Reject all channels if UKMet rain flag or gross Tb error detected
1962 if ( rainDetectionUKMethod(obsIndex) .or. grossRej(obsIndex) ) then
1963 ! "BAD" Observations
1964 obsToReject(:) = .true.
1965 numTotFilteredObs = numTotFilteredObs + 1
1966
1967 else ! if ( ukBadObs(obsIndex) .or. grossRej(obsIndex) )
1968 ! "GOOD" Observations
1969 !-----------------------------------------------------------------------------------
1970 ! OVER LAND OR SEA-ICE,
1971 ! -- reject lower tropospheric channels and imager channels
1972 ! -- check Bennartz SI and DI for AMSU-B like channels
1973 !-----------------------------------------------------------------------------------
1974 ! -- CLW/PRECIP not determined over land
1975 ! -- surface emissivity effects lower tropospheric channels
1976 if ( .not. waterObs(obsIndex) ) then
1977 obsToReject(1:ipc) = .true. ! AMSU-A 3-5(6)
1978 obsToReject(12:18) = .true. ! SSMI-like imager 1-7
1979 obsToReject(8:9) = .true. ! AMSU-B 2,5
1980 ! Check BSI and DI for AMSU-B 3-4
1981 ! Bennartz Scattering Index
1982 ! Land point
1983 if ( scatL(obsIndex) > 0.0 .and. scatL(obsIndex) /= ssbg_realMissing ) obsToReject(10:11) = .true.
1984 ! Sea-ice point
1985 if ( scatW(obsIndex) > 40.0 .and. scatW(obsIndex) /= ssbg_realMissing ) obsToReject(10:11) = .true.
1986 ! Missing scattering index
1987 if ( scatW(obsIndex) == ssbg_realMissing .and. scatL(obsIndex) == ssbg_realMissing ) obsToReject(10:11) = .true.
1988 if ( any(obsToReject(10:11))) then
1989 numLandScatObs = numLandScatObs + 1
1990 end if
1991 ! Dryness index
1992 if ( amsubDrynessIndex(obsIndex) > 0.0 ) then
1993 obsToReject(11) = .true.
1994 end if
1995 if ( amsubDrynessIndex(obsIndex) > -10.0 ) then
1996 obsToReject(10) = .true.
1997 end if
1998 if ( amsubDrynessIndex(obsIndex) > -20.0 ) then
1999 obsToReject(9) = .true.
2000 end if
2001 if ( amsubDrynessIndex(obsIndex) > -10.0 ) numDryIndexObs = numDryIndexObs + 1
2002 numTotFilteredObs = numTotFilteredObs + 1
2003 end if
2004 !-----------------------------------------------------------------------------------
2005 ! OVER WATER,
2006 ! -- reject tropospheric channels and imager channels if cloudy,
2007 ! precip, or excessive IWV
2008 ! -- check Bennartz SI for AMSU-B like channels
2009 !-----------------------------------------------------------------------------------
2010 if ( waterObs(obsIndex) ) then
2011 if (cloudObs(obsIndex) .or. iwvReject(obsIndex)) then
2012 if ( iwvReject(obsIndex) ) then
2013 ! SSMI-like imager channels
2014 obsToReject(12:18) = .true.
2015 ! AMSU-A like channels 3-5(6)
2016 obsToReject(1:ipc) = .true.
2017 else ! ----- CLOUDY OBSERVATION (OR MISSING CLOUD) ------
2018 ! SSMI-like imager channels
2019 obsToReject(12:18) = .true.
2020 ! AMSU-A like channels 3-5(6)
2021 if ( rclw(obsIndex) > clw_amsu_rej ) then
2022 obsToReject(2:ipc) = .true.
2023 end if
2024 if ( rclw(obsIndex) > clw_amsu_rej_ch3 ) then
2025 obsToReject(1) = .true.
2026 end if
2027 end if
2028 numTotFilteredObs = numTotFilteredObs + 1
2029 end if
2030 ! Check BSI for AMSU B channels 2-5 for all water obs
2031 ! Bennartz Scattering Index ( AMSU-B 2-5)
2032 ! Open water point
2033 if ( scatW(obsIndex) > 15.0 .and. scatW(obsIndex) /= ssbg_realMissing ) then
2034 obsToReject(8:11) = .true.
2035 numSeaScatObs = numSeaScatObs + 1
2036 if ( precipObs(obsIndex) ) numScatPrecipObs = numScatPrecipObs + 1
2037 end if
2038
2039 end if ! if waterObs
2040
2041 end if ! if ( ukBadObs(obsIndex) .or. grossRej(obsIndex) [end else] )
2042
2043 if ( .not. waterObs(obsIndex) ) numLandObs = numLandObs + 1
2044 if ( rainDetectionUKMethod(obsIndex) ) numUkBadObs = numUkBadObs + 1
2045 if ( grossRej(obsIndex) ) numGrossObs = numGrossObs + 1
2046 if ( cloudObs(obsIndex) ) numCloudsinObs = numCloudsinObs + 1
2047 if ( waterObs(obsIndex) ) numWaterObs = numWaterObs + 1
2048 if ( cloudObs(obsIndex) .and. waterObs(obsIndex) ) numCloudyObs = numCloudyObs + 1
2049 if ( iwvReject(obsIndex) .and. waterObs(obsIndex) .and. &
2050 & (.not. cloudObs(obsIndex)) ) numBadIWVObs = numBadIWVObs + 1
2051 if ( precipObs(obsIndex) .and. waterObs(obsIndex) ) then
2052 numPrecipObs = numPrecipObs + 1
2053 end if
2054 end do HEADER1
2055
2056 !-------------------------------------------------------------------------------
2057 ! Update ObsspaceData with the computed values of: ZLQ, ZTT, IWV and scatW
2058 !-------------------------------------------------------------------------------
2059 call obs_headSet_i(obsSpaceData, OBS_STYP, headerIndex, landSeaQualifier(1))
2060 call obs_headSet_i(obsSpaceData, OBS_TTYP, headerIndex, terrainType(1))
2061 if (scatW(1) /= ssbg_realMissing) then
2062 call obs_headSet_r(obsSpaceData, OBS_SIO, headerIndex, scatW(1))
2063 else
2064 call obs_headSet_r(obsSpaceData, OBS_SIO, headerIndex, MPC_missingValue_R4)
2065 end if
2066 call obs_headSet_r(obsSpaceData, OBS_CLWO, headerIndex, rclw(1))
2067 call obs_headSet_r(obsSpaceData, OBS_IWV, headerIndex, riwv(1))
2068
2069 if (headerIndex == obs_numHeader(obsSpaceData)) then
2070 write(*,*) '*******************************************************************'
2071 write(*,*) '******************* SATQC PROGRAM STATS****************************'
2072 write(*,*) '*******************************************************************'
2073 write(*,*)
2074 write(*,*) 'F16 obs = ', numObsF16
2075 write(*,*) 'F17 obs = ', numObsF17
2076 write(*,*) 'F18 obs = ', numObsF18
2077 write(*,*)
2078 write(*,*) 'Land obs = ', numLandObs
2079 write(*,*) 'UKMO bad obs = ', numUkBadObs
2080 write(*,*) 'Gross rej obs = ', numGrossObs
2081 write(*,*) 'Cloudy obs = ', numCloudyObs
2082 write(*,*) 'Bad IWV obs = ', numBadIWVObs
2083 write(*,*) 'Precip obs = ', numPrecipObs
2084 write(*,*) 'Clouds in obs = ', numCloudsinObs
2085 write(*,*) 'Obs with water = ', numWaterObs
2086 write(*,*) '-------------------------------------------------------------------'
2087 write(*,*) 'Total filtered obs = ', numTotFilteredObs
2088 write(*,*) '-------------------------------------------------------------------'
2089 write(*,*) ' AMSU-B-Like Channel Filtering Report'
2090 write(*,*) 'Land/ice points with data flagged for Bennartz SI = ', numLandScatObs
2091 write(*,*) 'Water points with data flagged for Bennartz SI = ', numSeaScatObs
2092 write(*,*) 'Dry obs = ', numDryIndexObs
2093 write(*,*) 'Sea Ice obs = ', numSeaIceObs
2094 write(*,*) '*******************************************************************'
2095 end if
2096
2097 end subroutine ssbg_satqcSsmis
2098
2099 !--------------------------------------------------------------------------
2100 ! ssbg_updateObsSpaceAfterSatQc
2101 !--------------------------------------------------------------------------
2102 subroutine ssbg_updateObsSpaceAfterSatQc(obsSpaceData, headerIndex, obsToReject)
2103 !
2104 !:Purpose: Update obspacedata variables (obstTB and obs flags) after QC
2105 !
2106 implicit none
2107
2108 ! Arguments:
2109 type(struct_obs), intent(inout) :: obsSpaceData ! ObsSpaceData object
2110 integer, intent(in) :: headerIndex ! Current header index
2111 logical, intent(in) :: obsToReject(:) ! Observations that will be rejected
2112
2113 ! Locals:
2114 integer, allocatable :: obsFlags(:)
2115 integer, allocatable :: obsGlobalFlag(:)
2116 integer, allocatable :: satScanPosition(:)
2117 integer :: actualNumChannel
2118 integer :: bodyIndex
2119 integer :: bodyIndexBeg
2120 integer :: bodyIndexEnd
2121 integer :: channelIndex
2122 integer :: currentChannelNumber
2123 integer :: dataIndex
2124 integer :: headerCompt
2125 integer :: instr
2126 integer :: instrId
2127 integer :: iSat
2128 integer :: numObsToProcess
2129 integer :: obsIndex
2130 integer :: platf
2131 integer :: platfId
2132 integer :: sensorIndex
2133 logical :: sensorIndexFound
2134
2135 ! find tvs_sensor index corresponding to current obs
2136
2137 platf = obs_headElem_i( obsSpaceData, OBS_SAT, headerIndex )
2138 instr = obs_headElem_i( obsSpaceData, OBS_INS, headerIndex )
2139
2140 call tvs_mapSat( platf, platfId, iSat )
2141 call tvs_mapInstrum( instr, instrId )
2142
2143 sensorIndexFound = .false.
2144 HEADER0: do sensorIndex = 1, tvs_nsensors
2145 if ( platfId == tvs_platforms(sensorIndex) .and. &
2146 iSat == tvs_satellites(sensorIndex) .and. &
2147 instrId == tvs_instruments(sensorIndex) ) then
2148 sensorIndexFound = .true.
2149 exit HEADER0
2150 end if
2151 end do HEADER0
2152 if ( .not. sensorIndexFound ) call utl_abort('ssbg_updateObsSpaceAfterSatQc: sensor Index not found')
2153
2154 ! find actual Number of channels
2155 actualNumChannel = tvs_coefs(sensorIndex)%coef%fmv_ori_nchn
2156
2157 headerCompt = 1
2158 numObsToProcess = 1
2159 ! Allocate Header elements
2160 call utl_reAllocate(obsGlobalFlag, numObsToProcess)
2161 call utl_reAllocate(satScanPosition, numObsToProcess)
2162
2163 ! Allocation
2164 call utl_reAllocate(obsFlags, numObsToProcess*actualNumChannel)
2165
2166 ! Read elements in obsspace
2167
2168 obsGlobalFlag(headerCompt) = obs_headElem_i( obsSpaceData, OBS_ST1, headerIndex )
2169 satScanPosition(headerCompt) = obs_headElem_i( obsSpaceData, OBS_FOV, headerIndex )
2170
2171 bodyIndexBeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
2172 bodyIndexEnd = bodyIndexBeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) -1
2173
2174 BODY2: do bodyIndex = bodyIndexbeg, bodyIndexEnd
2175 currentChannelNumber = nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex ))-tvs_channelOffset(sensorIndex)
2176 obsFlags(currentChannelNumber) = obs_bodyElem_i( obsSpaceData, OBS_FLG, bodyIndex )
2177 end do BODY2
2178
2179 !------------------------------------------------------------------
2180 ! 1 - Mark global flags if any obsToReject
2181 !------------------------------------------------------------------
2182 do obsIndex = 1, numObsToProcess
2183 if ( any(obsToReject) ) obsGlobalFlag(obsIndex) = ibset(obsGlobalFlag(obsIndex), 6)
2184 end do
2185
2186 !------------------------------------------------------------------
2187 ! 2 - Mark obs flags for each value of obsToReject
2188 !------------------------------------------------------------------
2189 dataIndex = 0
2190 do obsIndex = 1, numObsToProcess
2191 do channelIndex = 1, actualNumChannel
2192 dataIndex = dataIndex+1
2193 if (resetQc) obsFlags(dataIndex) = 0
2194 if (obsToReject(dataIndex)) obsFlags(dataIndex) = ibset(obsFlags(dataIndex),7)
2195 end do
2196 end do
2197
2198 !-----------------------------------------------------------------
2199 ! Subtract 270 from FOV values (element 005043).
2200 !-----------------------------------------------------------------
2201
2202 do obsIndex = 1, numObsToProcess
2203 if (satScanPosition(obsIndex) > 270) satScanPosition(obsIndex) = satScanPosition(obsIndex) - 270
2204 end do
2205
2206 ! write elements in obsspace
2207 call obs_headSet_i(obsSpaceData, OBS_FOV, headerIndex, satScanPosition(1))
2208 call obs_headSet_i(obsSpaceData, OBS_ST1, headerIndex, obsGlobalFlag(1))
2209
2210 bodyIndexBeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
2211 bodyIndexEnd = bodyIndexBeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) -1
2212
2213 BODY: do bodyIndex = bodyIndexBeg, bodyIndexEnd
2214 currentChannelNumber=nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex ))-tvs_channelOffset(sensorIndex)
2215 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyIndex, obsFlags(currentChannelNumber))
2216 end do BODY
2217
2218 end subroutine ssbg_updateObsSpaceAfterSatQc
2219
2220 !--------------------------------------------------------------------------
2221 ! ssbg_inovqcSsmis
2222 !--------------------------------------------------------------------------
2223 subroutine ssbg_inovqcSsmis(obsSpaceData, headerIndex, flagsInovQc)
2224 !
2225 !:Purpose: Identify those observations in SSMIS data that have O-P
2226 ! values greater than a threshold proportional to known
2227 ! standard deviations (computed when the bias correction
2228 ! coefficients were derived). The flags of these observations
2229 ! are adjusted accordingly (ie bit 9 switched ON).
2230 ! Also,
2231 ! -- flag channels for systematic rejection based on UTIL
2232 ! value in stats_*_errtot file or because flag bit 6 OFF
2233 ! (uncorrected data)
2234 ! -- reject sets of AMSU-like channels based on O-P for
2235 ! a single channel
2236 ! -- reject selected AMSU-like channels over land when
2237 ! model surface height exceeds a specified limit
2238 ! (topography check)
2239 !
2240 implicit none
2241
2242 ! Arguments:
2243 type(struct_obs), intent(inout) :: obsSpaceData ! ObsSpaceData object
2244 integer, intent(in) :: headerIndex ! Current header index
2245 integer, allocatable, intent(out) :: flagsInovQc(:) ! Flags for assimilation/rejection of obs
2246
2247 ! Locals:
2248 character(len=9) :: burpFileSatId ! Satellite ID
2249 integer, allocatable :: obsChannels(:) ! channel numbers
2250 integer, allocatable :: obsFlags(:) ! data flags
2251 integer :: actualNumChannel ! actual Num channel
2252 integer :: bodyIndex
2253 integer :: bodyIndexBeg
2254 integer :: bodyIndexEnd
2255 integer :: codtyp ! code type
2256 integer :: channelIndex
2257 integer :: currentChannelNumber
2258 integer :: headerCompt
2259 integer :: instr
2260 integer :: instrId
2261 integer :: iSat
2262 integer :: numObsToProcess ! Number of obs points to process
2263 integer :: platf
2264 integer :: platfId
2265 integer :: sensorIndex ! find tvs_sensor index corresponding to current obs
2266 logical :: sensorIndexFound
2267 logical :: ssmisDataPresent
2268 real , allocatable :: modelInterpTer(:) ! topo in standard file interpolated to obs point
2269 real , allocatable :: obsLatitude(:) ! obs. point latitude
2270 real , allocatable :: obsLongitude(:) ! obs. point longitude
2271 real , allocatable :: ompTb(:) ! OMP values
2272
2273 !-------------------------------------------------------------------------
2274 ! 1) INOVQC begins
2275 !-------------------------------------------------------------------------
2276
2277 ! Check if its ssmis data:
2278 ssmisDataPresent = .false.
2279 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
2280 if ( tvs_isIdBurpInst(codtyp,'ssmis') ) then
2281 ssmisDataPresent = .true.
2282 end if
2283
2284 if ( .not. ssmisDataPresent ) then
2285 write(*,*) 'WARNING: WILL NOT RUN ssbg_inovqcSsmis since no SSMIS DATA were found'
2286 return
2287 end if
2288
2289 ! find tvs_sensor index corresponding to current obs
2290
2291 platf = obs_headElem_i( obsSpaceData, OBS_SAT, headerIndex )
2292 instr = obs_headElem_i( obsSpaceData, OBS_INS, headerIndex )
2293
2294 call tvs_mapSat( platf, platfId, iSat )
2295 call tvs_mapInstrum( instr, instrId )
2296
2297 sensorIndexFound = .false.
2298 HEADER1: do sensorIndex = 1, tvs_nsensors
2299 if ( platfId == tvs_platforms(sensorIndex) .and. &
2300 iSat == tvs_satellites(sensorIndex) .and. &
2301 instrId == tvs_instruments(sensorIndex) ) then
2302 sensorIndexFound = .true.
2303 exit HEADER1
2304 end if
2305 end do HEADER1
2306 if ( .not. sensorIndexFound ) call utl_abort('ssbg_inovqcSsmis: sensor Index not found')
2307
2308 !--------------------------------------------------------------------
2309 ! 2) Allocating arrays
2310 !--------------------------------------------------------------------
2311
2312 ! find actual Number of channels
2313 actualNumChannel = tvs_coefs(sensorIndex)%coef%fmv_ori_nchn
2314
2315 headerCompt = 1
2316 numObsToProcess = 1
2317
2318 ! ELEMENTS FROM OBSSPACEDATA
2319 ! Allocate header elements
2320 call utl_reAllocate(obsLatitude, numObsToProcess)
2321 call utl_reAllocate(obsLongitude, numObsToProcess)
2322 ! Allocate Body elements
2323 call utl_reAllocate(obsChannels, numObsToProcess*actualNumChannel)
2324 call utl_reAllocate(obsFlags, numObsToProcess*actualNumChannel)
2325 call utl_reAllocate(ompTb, numObsToProcess*actualNumChannel)
2326 ! Allocate intent out array
2327 call utl_reAllocate(flagsInovQc, numObsToProcess*actualNumChannel)
2328
2329 ! Lecture dans obsspacedata
2330 burpFileSatId = obs_elem_c ( obsSpaceData, 'STID' , headerIndex )
2331 obsLatitude(headerCompt) = obs_headElem_r( obsSpaceData, OBS_LAT, headerIndex )
2332 obsLongitude(headerCompt) = obs_headElem_r( obsSpaceData, OBS_LON, headerIndex )
2333
2334 bodyIndexBeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
2335 bodyIndexEnd = bodyIndexBeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) -1
2336
2337 BODY: do bodyIndex = bodyIndexbeg, bodyIndexEnd
2338 currentChannelNumber = nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex ))-tvs_channelOffset(sensorIndex)
2339 ompTb(currentChannelNumber) = obs_bodyElem_r( obsSpaceData, OBS_OMP, bodyIndex )
2340 obsFlags(currentChannelNumber) = obs_bodyElem_i( obsSpaceData, OBS_FLG, bodyIndex )
2341 end do BODY
2342
2343 do channelIndex=1, actualNumChannel
2344 obsChannels(channelIndex) = channelIndex+tvs_channelOffset(sensorIndex)
2345 end do
2346
2347 ! Convert lat/lon to degrees
2348 obsLatitude(headerCompt) = obsLatitude(headerCompt) *MPC_DEGREES_PER_RADIAN_R8
2349 obsLongitude(headerCompt) = obsLongitude(headerCompt)*MPC_DEGREES_PER_RADIAN_R8
2350 if( obsLongitude(headerCompt) > 180. ) obsLongitude(headerCompt) = obsLongitude(headerCompt) - 360.
2351
2352 !--------------------------------------------------------------------
2353 ! 3) Extract model topography from the input GEO file and interpolate
2354 ! to all observation points in the box
2355 !--------------------------------------------------------------------
2356
2357 call ssbg_readGeophysicFieldsAndInterpolate(obsLatitude, obsLongitude, modelInterpTer)
2358
2359 !--------------------------------------------------------------------
2360 ! 4) Perform quality control on the data, checking O-P and topography
2361 !--------------------------------------------------------------------
2362
2363 call check_stddev(obsChannels, ompTb, flagsInovQc, actualNumChannel, numObsToProcess, &
2364 & sensorIndex, burpFileSatId, obsFlags)
2365
2366 call check_topo(modelInterpTer, flagsInovQc, actualNumChannel, numObsToProcess)
2367
2368 end subroutine ssbg_inovqcSsmis
2369
2370 !--------------------------------------------------------------------------
2371 ! check_stddev
2372 !--------------------------------------------------------------------------
2373 subroutine check_stddev(obsChannels, ompTb, flagsInovQc, actualNumChannel, numObsToProcess, &
2374 sensorIndex, burpFileSatId, obsFlags)
2375 !
2376 !:Purpose: Perform quality control on the radiances by analysing the
2377 ! magnitude of the residuals.
2378
2379 !------------------------------------------------------------------
2380 ! Variable Definitions:
2381 ! ---------------------
2382 ! obsChannels - input - channel numbers (residuals)
2383 ! ompTb - input - radiance residuals (o-p)
2384 ! obsFlags - input - radiance data flags (bit 7 on for data that
2385 ! will not be assimilated)
2386 ! flagsInovQc - output - quality contol indicator for each channel of each
2387 ! observation point
2388 ! =0 ok
2389 ! =1 not checked because FLAG bit 7 ON,
2390 ! =2 reject by UTIL value or FLAG bit 6 OFF (data not bias corrected),
2391 ! =3 reject by rogue check,
2392 ! =4 reject by both UTIL value and rogue check.
2393 ! ==> individual channels may be rejected from each obs pt
2394 ! actualNumChannel - input - number of residual channels
2395 ! numObsToProcess - input - number of groups of NVAL*NELE
2396 ! sensorIndex - input - number of satellite (index # --> 1-nsat)
2397 ! burpFileSatId - input - identificateur du satellite
2398 ! rogueFac -internal- constant which determines the severity of the
2399 ! quality control applied to the O-P magnitude
2400 ! for each channel
2401 ! productRogueSTD -internal- product of rogueFac and standard deviation for
2402 ! each channel
2403 ! oer_tovutil -external- UTIL values read from the total error statistics file
2404 ! 0 = blacklisted channel
2405 ! 1 = assimilated channel
2406 ! oer_toverrst -external- standard deviation statistics read from the total
2407 ! error statistics file
2408 !
2409 implicit none
2410
2411 ! Arguments:
2412 integer, intent(in) :: obsChannels(:) ! Channel numbers
2413 real , intent(in) :: ompTb(:) ! Radiance residuals
2414 integer, intent(out) :: flagsInovQc(:) ! Flags for assimilation/rejection of obs
2415 integer, intent(in) :: actualNumChannel ! Number of channels
2416 integer, intent(in) :: numObsToProcess ! Number of obs points to process
2417 integer, intent(in) :: sensorIndex ! Identification number of satellite
2418 character(len=9), intent(in) :: burpFileSatId ! Satellite identification in BURP file
2419 integer, intent(in) :: obsFlags(:) ! Radiance data flags
2420
2421 ! Locals:
2422 real, parameter :: factorCh1 = 2.0 ! factor for channel 1 O-P for rejection of channels 1-4
2423 real, parameter :: maxOmpCh8 = 5.0 ! max Abs(O-P) for channel 8 for rejection of channels 8-11 (units = K)
2424 integer :: chanIndex
2425 integer :: chanIndex1
2426 integer :: chanIndex4
2427 integer :: chanIndex8
2428 integer :: chanIndex11
2429 integer :: chanIndexLast
2430 integer :: channelNumber
2431 integer :: obsChanIndex
2432 integer :: obsIndex
2433 real :: productRogueSTD
2434 real :: rogueFac(ssbg_maxNumChan)
2435
2436 ! Initialization. Assume all observations are to be kept.
2437
2438 ! rogueFac(1:7) should be consistent with bgck.satqc_amsua.f (for AMSUA ch. 3-10)
2439 ! rogueFac(8:11) should be consistent with bgck.satqc_amsub.f (for AMSUB ch. 2-5)
2440 ! rogueFac(12:18) should be consistent with ssmi_inovqc.ftn90 (for SSMI ch. 1-7)
2441 ! rogueFac(19:24) are for upper-atmospheric channels (strato/meso-sphere)
2442 rogueFac = (/2.0, 3.0, 4.0, 4.0, 4.0, 4.0, 4.0, 2.0, &
2443 & 4.0, 4.0, 4.0, 2.0, 2.0, 2.0, 2.0, 2.0, &
2444 & 2.0, 2.0, 3.0, 3.0, 3.0, 3.0, 3.0, 3.0/)
2445
2446 flagsInovQc(:) = 0
2447
2448 ! Loop through all observation points in the current record and check
2449 ! first whether a channel is to be systematically removed and second
2450 ! whether the O-P value for each channel falls within the acceptable
2451 ! range defined by rogueFac and oer_toverrst (sensorIndex identifies the satellite).
2452
2453 ! *** Only do check for data that could be assimilated ***
2454
2455 ! Also, apply multi-channel rejections for AMSU-like channels using
2456 ! O-P values from channels 1 and 8.
2457 ! If Abs(O-P) chan. 1 (AMSU-A 3) > 2*errtot, reject ch. 1-4 (AMSU-A 3-6)
2458 ! If Abs(O-P) chan. 8 (AMSU-B 2) > 5K, reject ch. 8-11 (AMSU-B 2-5)
2459
2460 !--------------------------------------------------------------------------
2461 ! 1. Rogue check O-P for all data flagged for assimilation
2462 !--------------------------------------------------------------------------
2463
2464 HEADER0: do obsChanIndex = 1, numObsToProcess*actualNumChannel
2465
2466 channelNumber = obsChannels(obsChanIndex)
2467 productRogueSTD = rogueFac(channelNumber) * oer_toverrst(channelNumber,sensorIndex)
2468
2469 if ( .not. btest(obsFlags(obsChanIndex),7) ) then
2470
2471 if ( oer_tovutil(channelNumber,sensorIndex) == 0 .or. (.not. btest(obsFlags(obsChanIndex),6)) ) then
2472
2473 ! systematic rejection of this channel
2474 flagsInovQc(obsChanIndex) = 2
2475 if ( abs( ompTb(obsChanIndex) ) >= productRogueSTD ) then
2476 ! rogue check failure as well
2477 flagsInovQc(obsChanIndex) = 4
2478 end if
2479
2480 else
2481
2482 ! keep channel, now perform rogue check
2483 if ( abs( ompTb(obsChanIndex) ) >= productRogueSTD ) then
2484 flagsInovQc(obsChanIndex) = 3
2485 end if
2486
2487 end if
2488
2489 else
2490
2491 flagsInovQc(obsChanIndex) = 1
2492
2493 end if
2494
2495 ! Keep statistics of obs rejected by rogue check.
2496 if ( (flagsInovQc(obsChanIndex) > 2) .and. (ssbg_debug) ) then
2497 write(*,*) 'check_stddev: '
2498 write(*,*) burpFileSatId(2:9),' Rogue check reject: ', &
2499 & ' Obs = ',obsChanIndex,' Channel = ',channelNumber, &
2500 & ' Check value = ',productRogueSTD, &
2501 & ' O-P = ',ompTb(obsChanIndex)
2502 end if
2503
2504 end do HEADER0
2505
2506 !----------------------------------------------------------------------------------
2507 ! 2. Apply additional criteria for rejections of multiple AMSU-like channels using
2508 ! Channels 1 and 8 Abs(O-P) (for data that could be assimilated [bit 7 off])
2509 !----------------------------------------------------------------------------------
2510
2511 ! Loop over observation points (actualNumChannel = num. channels [24])
2512 HEADER1: do obsIndex = 1, numObsToProcess
2513
2514 chanIndexLast = obsIndex*actualNumChannel ! index of ch.24 obs
2515 chanIndex1 = chanIndexLast - (actualNumChannel-1) ! index of ch.1 obs
2516 chanIndex4 = chanIndex1 + 3 ! index of ch.4 obs
2517 chanIndex8 = chanIndex4 + 4 ! index of ch.8 obs
2518 chanIndex11 = chanIndex8 + 3 ! index of ch.11 obs
2519
2520 ! AMSU-A like ch. 1-4
2521 if ( sum(flagsInovQc(chanIndex1:chanIndex4)) /= 4 ) then
2522 productRogueSTD = factorCh1 * oer_toverrst(1,sensorIndex)
2523 if ( abs(ompTb(chanIndex1)) >= productRogueSTD ) then
2524 do chanIndex = chanIndex1, chanIndex4
2525 channelNumber = obsChannels(chanIndex)
2526 if (flagsInovQc(chanIndex) /= 1) flagsInovQc(chanIndex) = max(flagsInovQc(chanIndex),3)
2527 end do
2528 end if
2529 end if
2530
2531 ! AMSU-B like ch. 8-11
2532 if ( sum(flagsInovQc(chanIndex8:chanIndex11)) /= 4 ) then
2533 productRogueSTD = maxOmpCh8
2534 if ( abs(ompTb(chanIndex8)) >= productRogueSTD ) then
2535 do chanIndex = chanIndex8, chanIndex11
2536 channelNumber = obsChannels(chanIndex)
2537 if (flagsInovQc(chanIndex) /= 1) flagsInovQc(chanIndex) = max(flagsInovQc(chanIndex),3)
2538 end do
2539 end if
2540 end if
2541
2542 end do HEADER1
2543
2544 end subroutine check_stddev
2545
2546 !--------------------------------------------------------------------------
2547 ! check_topo
2548 !--------------------------------------------------------------------------
2549 subroutine check_topo(modelInterpTer, flagsInovQc, actualNumChannel, numObsToProcess)
2550
2551 !:Purpose: Perform rejection of observations for selected channels based
2552 ! on model surface height (for channels assimilated over land)
2553
2554 ! -- for a single satellite (burpFileSatId,sensorIndex)
2555 ! -- for a single box (nt observations)
2556
2557 !------------------------------------------------------------------
2558 ! Variable Definitions:
2559 ! ---------------------
2560 ! modelInterpTer - input - model surface height (m) at each observation point
2561 ! flagsInovQc - in/out - quality control indicator for each channel of each
2562 ! observation point
2563 ! -- on INPUT
2564 ! =0 ok
2565 ! =1 not checked because FLAG bit 7 ON,
2566 ! =2 reject by UTIL value or bit 6 OFF (data not bias corrected),
2567 ! =3 reject by rogue check,
2568 ! =4 reject by both UTIL value and rogue check.
2569 ! -- on OUTPUT,
2570 ! =0 ok
2571 ! =1 not checked because FLAG bit 7 ON,
2572 ! =2 reject by UTIL value,
2573 ! =3 reject by rogue check,
2574 ! =4 reject by both UTIL value and rogue check.
2575 ! =5 rejection due to topography
2576 ! =6 rejection due to topography, reject by UTIL value
2577 ! =7 rejection due to topography, reject by rogue check
2578 ! =8 rejection due to topography, reject by UTIL value,
2579 ! reject by rogue check
2580 ! ==> individual channels may be rejected from each obs pt
2581 ! actualNumChannel - input - number of residual channels
2582 ! numObsToProcess - input - number of groups of NVAL*NELE
2583 !
2584 implicit none
2585
2586 ! Arguments:
2587 real, intent(in) :: modelInterpTer(:) ! Model surface height (m) for each obs
2588 integer, intent(inout) :: flagsInovQc(:) ! Flags for assimilation/rejection of obs
2589 integer, intent(in) :: actualNumChannel ! Number of channels
2590 integer, intent(in) :: numObsToProcess ! Number of obs points to process
2591
2592 ! Locals:
2593 integer, parameter :: nChanCheck=4 ! number of channels to check
2594 integer :: chanIndex
2595 integer :: chanIndex1
2596 integer :: checkedChan(nChanCheck)
2597 integer :: checkedChanIndex
2598 integer :: obsIndex
2599 integer :: topoReject
2600 logical :: debugSupp
2601 real :: heightLimit
2602 real :: topoHeight
2603 real :: topoLimit(nChanCheck)
2604
2605 !------------------------------------------------------------------
2606 ! Define channels to check and height limits (m) for rejection
2607 !------------------------------------------------------------------
2608 checkedChan = (/ 4, 9, 10, 11 /)
2609 topoLimit = (/ 250., 1000., 2000., 2500. /)
2610
2611 if (ssbg_debug) then
2612 write(*,*) 'check_topo: Maximum input GZ = ', maxval(modelInterpTer)
2613 write(*,*) 'check_topo: numObsToProcess = ', numObsToProcess
2614 end if
2615
2616 !------------------------------------------------------------------
2617 ! Perform the check for the selected channels
2618 !------------------------------------------------------------------
2619
2620 ! Loop over all numObsToProcess observation locations
2621
2622 topoReject = 0
2623 HEADER: do obsIndex = 1, numObsToProcess
2624 debugSupp = .false.
2625 chanIndex1 = (obsIndex-1)*actualNumChannel + 1 ! index of ch.1 obs
2626 topoHeight = modelInterpTer(obsIndex) ! model topography height [m]
2627
2628 if (ssbg_debug) then
2629 if (topoHeight == maxval(modelInterpTer) .and. topoHeight > minval(topoLimit)) then
2630 write(*,*) 'check_topo: ****** Max height point! topoHeight = ', topoHeight
2631 debugSupp = .true.
2632 end if
2633 end if
2634
2635 SUBHEADER: do checkedChanIndex = 1, nChanCheck ! loop over channels to check (checkedChan)
2636 heightLimit = topoLimit(checkedChanIndex) ! height limit [m] for channel checkedChan(checkedChanIndex)
2637 chanIndex = chanIndex1 + (checkedChan(checkedChanIndex)-1) ! channel checkedChan(checkedChanIndex) index
2638 if ( flagsInovQc(chanIndex) /= 1 ) then
2639 if ( topoHeight > heightLimit ) then
2640 flagsInovQc(chanIndex) = max(1,flagsInovQc(chanIndex)) + 4
2641 if (ssbg_debug .and. debugSupp) write(*,*) 'check_topo: Incrementing topoReject for max topoHeight point for ch.= ', checkedChan(checkedChanIndex)
2642 topoReject = topoReject + 1
2643 if (ssbg_debug) then
2644 if ( topoReject <= nChanCheck ) then
2645 write(*,*) 'check_topo:'
2646 write(*,*) ' Channel = ', checkedChan(checkedChanIndex), &
2647 & ' Height limit (m) = ', heightLimit, &
2648 & ' Model height (m) = ', topoHeight
2649 end if
2650 end if
2651 end if
2652 end if
2653 end do SUBHEADER
2654
2655 end do HEADER
2656
2657 if (ssbg_debug .and. (topoReject > 0) ) then
2658 write(*,*) 'check_topo: Number of topography rejections and observations for this box = ', topoReject, numObsToProcess*actualNumChannel
2659 end if
2660
2661 end subroutine check_topo
2662
2663 !--------------------------------------------------------------------------
2664 ! ssbg_updateObsSpaceAfterInovQc
2665 !--------------------------------------------------------------------------
2666 subroutine ssbg_updateObsSpaceAfterInovQc(obsSpaceData, headerIndex, flagsInovQc)
2667 !
2668 !:Purpose: Update obspacedata variables (obstTB and obs flags) after QC
2669 !
2670 implicit none
2671
2672 ! Arguments:
2673 type(struct_obs), intent(inout) :: obsSpaceData ! ObsSpaceData object
2674 integer, intent(in) :: headerIndex ! Current header index
2675 integer, intent(in) :: flagsInovQc(:) ! Flags for assimilation/rejection of obs
2676
2677 ! Locals:
2678 integer, allocatable :: obsFlags(:)
2679 integer, allocatable :: obsGlobalFlag(:)
2680 integer, allocatable :: satScanPosition(:)
2681 logical :: sensorIndexFound
2682 integer :: actualNumChannel
2683 integer :: bodyIndex
2684 integer :: bodyIndexBeg
2685 integer :: bodyIndexEnd
2686 integer :: channelIndex
2687 integer :: currentChannelNumber
2688 integer :: dataIndex
2689 integer :: headerCompt
2690 integer :: instr
2691 integer :: instrId
2692 integer :: iSat
2693 integer :: numObsToProcess
2694 integer :: obsIndex
2695 integer :: platf
2696 integer :: platfId
2697 integer :: sensorIndex
2698
2699 ! find tvs_sensor index corresponding to current obs
2700
2701 platf = obs_headElem_i( obsSpaceData, OBS_SAT, headerIndex )
2702 instr = obs_headElem_i( obsSpaceData, OBS_INS, headerIndex )
2703
2704 call tvs_mapSat( platf, platfId, iSat )
2705 call tvs_mapInstrum( instr, instrId )
2706
2707 sensorIndexFound = .false.
2708 HEADER: do sensorIndex = 1, tvs_nsensors
2709 if ( platfId == tvs_platforms(sensorIndex) .and. &
2710 iSat == tvs_satellites(sensorIndex) .and. &
2711 instrId == tvs_instruments(sensorIndex) ) then
2712 sensorIndexFound = .true.
2713 exit HEADER
2714 end if
2715 end do HEADER
2716 if ( .not. sensorIndexFound ) call utl_abort('ssbg_updateObsSpaceAfterInovQc: sensor Index not found')
2717
2718 ! find actual Number of channels
2719 actualNumChannel = tvs_coefs(sensorIndex)%coef%fmv_ori_nchn
2720
2721 headerCompt = 1
2722 numObsToProcess = 1
2723 ! Allocate Header elements
2724 call utl_reAllocate(obsGlobalFlag, numObsToProcess)
2725 call utl_reAllocate(satScanPosition, numObsToProcess)
2726
2727 ! Allocation
2728 call utl_reAllocate(obsFlags,numObsToProcess*actualNumChannel)
2729
2730 ! Read elements in obsspace
2731
2732 obsGlobalFlag(headerCompt) = obs_headElem_i( obsSpaceData, OBS_ST1, headerIndex )
2733 satScanPosition(headerCompt) = obs_headElem_i( obsSpaceData, OBS_FOV, headerIndex )
2734 bodyIndexBeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
2735 bodyIndexEnd = bodyIndexBeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) -1
2736
2737 BODY2: do bodyIndex = bodyIndexbeg, bodyIndexEnd
2738 currentChannelNumber = nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex ))-tvs_channelOffset(sensorIndex)
2739 obsFlags(currentChannelNumber) = obs_bodyElem_i( obsSpaceData, OBS_FLG, bodyIndex )
2740 end do BODY2
2741 ! Modify flags
2742
2743 !------------------------------------------------------------------
2744 ! 1 - Mark global flags if any flagsInovQc > 0
2745 !------------------------------------------------------------------
2746 do obsIndex = 1 , numObsToProcess
2747 if ( maxval(flagsInovQc) > 0 ) obsGlobalFlag(obsIndex) = ibset(obsGlobalFlag(obsIndex), 6)
2748 end do
2749
2750 !------------------------------------------------------------------
2751 ! 2 - Mark obs flags for each value of flagsInovQc
2752 !------------------------------------------------------------------
2753
2754 dataIndex = 0
2755 do obsIndex = 1 , numObsToProcess
2756 do channelIndex = 1, actualNumChannel
2757 dataIndex = dataIndex+1
2758 if (resetQc) obsFlags(dataIndex) = 0
2759
2760 select case (flagsInovQc(dataIndex))
2761 case(1)
2762 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),9)
2763 case(2)
2764 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),11)
2765 case(3)
2766 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),9)
2767 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),16)
2768 case(4)
2769 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),9)
2770 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),11)
2771 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),16)
2772 case(5)
2773 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),9)
2774 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),18)
2775 case(6)
2776 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),9)
2777 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),11)
2778 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),18)
2779 case(7)
2780 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),9)
2781 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),16)
2782 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),18)
2783 case(8)
2784 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),9)
2785 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),11)
2786 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),16)
2787 obsFlags(dataIndex) = ibset(obsFlags(dataIndex),18)
2788 end select
2789
2790 end do
2791 end do
2792
2793 !-----------------------------------------------------------------
2794 ! Subtract 270 from FOV values (element 005043).
2795 !-----------------------------------------------------------------
2796 do obsIndex = 1 , numObsToProcess
2797 if (satScanPosition(obsIndex) > 270) satScanPosition(obsIndex) = satScanPosition(obsIndex) - 270
2798 end do
2799
2800 ! write elements in obsspace
2801 call obs_headSet_i(obsSpaceData, OBS_FOV, headerIndex, satScanPosition(1))
2802 call obs_headSet_i(obsSpaceData, OBS_ST1, headerIndex, obsGlobalFlag(1))
2803
2804 bodyIndexBeg = obs_headElem_i( obsSpaceData, OBS_RLN, headerIndex )
2805 bodyIndexEnd = bodyIndexBeg + obs_headElem_i( obsSpaceData, OBS_NLV, headerIndex ) -1
2806 BODY: do bodyIndex = bodyIndexBeg, bodyIndexEnd
2807 currentChannelNumber=nint(obs_bodyElem_r( obsSpaceData, OBS_PPP, bodyIndex ))-tvs_channelOffset(sensorIndex)
2808 call obs_bodySet_i(obsSpaceData, OBS_FLG, bodyIndex, obsFlags(currentChannelNumber))
2809 end do BODY
2810
2811 end subroutine ssbg_updateObsSpaceAfterInovQc
2812
2813 !--------------------------------------------------------------------------
2814 ! ssbg_bgCheckSSMIS
2815 !--------------------------------------------------------------------------
2816 subroutine ssbg_bgCheckSSMIS(obsSpaceData)
2817 !
2818 !:Purpose: Do the background check for SSMIS data (satQC and inovQC).
2819 !
2820 implicit none
2821
2822 ! Arguments:
2823 type(struct_obs), intent(inout) :: obsSpaceData ! ObsSpaceData object
2824
2825 ! Locals:
2826 integer, allocatable :: flagsInovQc(:)
2827 logical, allocatable :: obsToReject(:)
2828 integer :: codtyp
2829 integer :: dataIndex
2830 integer :: dataIndex1
2831 integer :: headerIndex
2832 integer :: indexFlags
2833 integer :: inovQcSize
2834 integer :: statsInovQcFlags(10)
2835 logical :: otherDataPresent
2836 logical :: ssmisDataPresent
2837 real :: percentInovQcFlags(9)
2838
2839 write(*,*) 'ssbg_bgCheckSSMIS: Starting'
2840
2841 call utl_tmg_start(119,'--BgckSSMIS')
2842 otherDataPresent = .false.
2843 ssmisDataPresent = .false.
2844 call obs_set_current_header_list(obsSpaceData,'TO')
2845 HEADER0: do
2846 headerIndex = obs_getHeaderIndex(obsSpaceData)
2847 if (headerIndex < 0) exit HEADER0
2848 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
2849 if ( tvs_isIdBurpInst(codtyp,'ssmis' ) ) then
2850 ssmisDataPresent = .true.
2851 else
2852 otherDataPresent = .true.
2853 end if
2854 end do HEADER0
2855
2856 if ( ssmisDataPresent .and. otherDataPresent ) then
2857 call utl_abort ('ssbg_bgCheckSSMIS: Other data than SSMIS also included in obsSpaceData')
2858 endif
2859
2860 if ( .not. ssmisDataPresent ) then
2861 write(*,*) 'WARNING: WILL NOT RUN ssbg_bgCheckSSMIS since no SSMIS'
2862 return
2863 end if
2864
2865 statsInovQcFlags(:) = 0
2866 percentInovQcFlags(:) = 0.0
2867
2868 ! read nambgck
2869 call ssbg_init()
2870 !Quality Control loop over all observations
2871 !
2872 ! loop over all header indices of the specified family with surface obs
2873
2874 call obs_set_current_header_list(obsSpaceData,'TO')
2875 HEADER: do
2876 headerIndex = obs_getHeaderIndex(obsSpaceData)
2877 if (headerIndex < 0) exit HEADER
2878 codtyp = obs_headElem_i(obsSpaceData, OBS_ITY, headerIndex)
2879 if ( .not. (tvs_isIdBurpInst(codtyp,'ssmis')) ) then
2880 write(*,*) 'WARNING in ssbg_bgCheckSSMIS: Observation with codtyp = ', codtyp, ' is not SSMIS'
2881 cycle HEADER
2882 end if
2883
2884 !###############################################################################
2885 ! STEP 1) call satQC SSMIS program !
2886 !###############################################################################
2887 if (ssbg_debug) write(*,*) 'ssbg_bgCheckSSMIS: STEP 1) call satQC SSMIS program'
2888 call ssbg_satqcSsmis(obsSpaceData, headerIndex, obsToReject)
2889
2890 !###############################################################################
2891 ! STEP 2) update Flags after satQC SSMIS program !
2892 !###############################################################################
2893 if (ssbg_debug) write(*,*) 'ssbg_bgCheckSSMIS: STEP 2) update Flags after satQC SSMIS program'
2894 call ssbg_updateObsSpaceAfterSatQc(obsSpaceData, headerIndex, obsToReject)
2895
2896 !###############################################################################
2897 ! STEP 3) call inovQC SSMIS program !
2898 !###############################################################################
2899 if (ssbg_debug) write(*,*) 'ssbg_bgCheckSSMIS: STEP 3) call inovQC SSMIS program'
2900 call ssbg_inovqcSsmis(obsSpaceData, headerIndex, flagsInovQc)
2901
2902 !###############################################################################
2903 ! STEP 4) update Flags after inovQC SSMIS program !
2904 !###############################################################################
2905 if (ssbg_debug) write(*,*) 'ssbg_bgCheckSSMIS: STEP 4) update Flags after inovQC SSMIS program'
2906 call ssbg_updateObsSpaceAfterInovQc(obsSpaceData, headerIndex, flagsInovQc)
2907
2908 !###############################################################################
2909 ! STEP 5) compute statistics of different inovQc flags types !
2910 !###############################################################################
2911 inovQcSize = size(flagsInovQc)
2912 if (maxval(flagsInovQc) > 8) call utl_abort('ssbg_bgCheckSSMIS: Problem with flagsInovQc, value greater than 8.')
2913 do dataIndex = 1,inovQcSize
2914 dataIndex1 = flagsInovQc(dataIndex)+1
2915 ! Counting number of flags with value flagsInovQc(dataIndex)
2916 statsInovQcFlags(dataIndex1) = statsInovQcFlags(dataIndex1) + 1
2917 ! Counting total number of flags (observations)
2918 statsInovQcFlags(10) = statsInovQcFlags(10) + 1
2919 end do
2920
2921 end do HEADER
2922
2923 !###############################################################################
2924 ! STEP 6) displaying statistics of inovQc flags !
2925 !###############################################################################
2926
2927 ! If data were not checked because all had the flag bit 7 ON, then set all percentage statistics to 0.
2928 if ( statsInovQcFlags(10) == statsInovQcFlags(2) ) then
2929 percentInovQcFlags(:) = 0.0
2930 else
2931 do indexFlags = 1,9
2932 percentInovQcFlags(indexFlags) = float(statsInovQcFlags(indexFlags))/float(statsInovQcFlags(10)-statsInovQcFlags(2))*100
2933 end do
2934 end if
2935
2936 write(*,*) '------------------- Innovation Rejection Statistics -----------------------'
2937 write(*,*) ' Flag description Nm. Obs. Prct (%) '
2938 write(*,256) ' Total number of observations : ', statsInovQcFlags(10)
2939 write(*,256) ' (1) Not checked because flag bit 7 ON : ', statsInovQcFlags(2)
2940 write(*,256) ' Remaining number of observations : ', statsInovQcFlags(10)-statsInovQcFlags(2)
2941 write(*,257) ' (0) Observations that are OK : ', statsInovQcFlags(1), percentInovQcFlags(1)
2942 write(*,257) ' (2) Rejected by UTIL value or bit 6 OFF : ', statsInovQcFlags(3), percentInovQcFlags(3)
2943 write(*,257) ' (3) Rejected by rogue check (O-P) : ', statsInovQcFlags(4), percentInovQcFlags(4)
2944 write(*,257) ' (4) Rejected by both UTIL and rogue check : ', statsInovQcFlags(5), percentInovQcFlags(5)
2945 write(*,257) ' (5) Topography rejection : ', statsInovQcFlags(6), percentInovQcFlags(6)
2946 write(*,257) ' (6) Topography rejection and by UTIL value : ', statsInovQcFlags(7), percentInovQcFlags(7)
2947 write(*,257) ' (7) Topography rejection and by rogue check : ', statsInovQcFlags(8), percentInovQcFlags(8)
2948 write(*,257) ' (8) Topography rejection, by UTIL and rogue check : ', statsInovQcFlags(9), percentInovQcFlags(9)
2949 write(*,*) '---------------------------------------------------------------------------'
2950
2951256 format(A55,i9)
2952257 format(A55,i9,f7.2,' %')
2953
2954 call utl_tmg_stop(119)
2955
2956 write(*,*) 'ssbg_bgCheckSSMIS: Finished'
2957
2958 end subroutine ssbg_bgCheckSSMIS
2959
2960end module bgckSSMIS_mod
