MODULE module_thin_ob
USE module_type
USE module_func
USE module_mm5
USE module_map
USE map_utils
USE module_namelist
CONTAINS
SUBROUTINE THIN_OB (max_number_of_obs, obs, number_of_obs, mix, mjx, &
missing_r, Ob_name, Ob_fm, Delta_P)
! --------------------------------------------------------------------------
! Purpose: To thin the sateliite (SATOB, SSMI Retrieval and Tb) data.
! Method : To select only ONE data nearest to the center of the (I,J,K)
! grid-box if there are more than one OBS within that gridbox.
!
! The vertical P_thickness is defined by Delta_P.
! Currently for SATOB the maximum of 100 level in vertical
! and maximum of 20 data points within a grid box are allowed.
! For the single level SSMI the maximum of 400 data points
! within a gridbox is allowed.
!
! Yong-Run Guo 05/01/2001
! --------------------------------------------------------------------------
!
implicit none
INTEGER, INTENT (in) :: max_number_of_obs
TYPE (report), DIMENSION (max_number_of_obs), INTENT (inout) :: obs
INTEGER, INTENT (in) :: number_of_obs
CHARACTER (LEN = *), INTENT (in) :: Ob_name
INTEGER, INTENT (in) :: Ob_fm
REAL, INTENT (in) :: missing_r
REAL , OPTIONAL, INTENT (in) :: Delta_P
INTEGER :: n, fm, i, j, k, kk, num_box, &
max_num, mix, mjx, num, nn, total, &
max_box, LV, N_selected, N1_selected
INTEGER, ALLOCATABLE :: Num_ob(:,:,:), Ob_index(:,:,:,:)
REAL :: X, Y, dx, dy, RR, Rmin, R1min
TYPE ( meas_data ) :: new
INTEGER :: nvalids_fm
!------------------------------------------------------------------
!
! 1.0 Allocate the working arrays
!
! The default size: LV, max_box can be changed to meet the needs.
!
! For SATOB
if (Ob_fm == 88) then
LV = 100
max_box = 50
! For SSMI
else if (Ob_fm == 125 .or. Ob_fm == 126 .or. Ob_fm == 281) then
LV = 1
max_box = 400
else
WRITE(0,'(A,I3,A)') ' FM=',Ob_fm,' NOT SATELLITE OBS,', &
' NO THINING APPLIED TO IT.'
return
endif
allocate (Num_ob (mix,mjx,LV))
allocate (Ob_index(mix,mjx,LV,max_box))
WRITE(0,'(/"Thining OBS ==> ",A,1X,"FM=",I3,2X, &
& "mix, mjx, LV, max_box:",4I6)') &
& Ob_name, Ob_fm, mix, mjx, LV, max_box
!
! 2.0 To count the number of OBS within the grid boxes
! ------------------------------------------------
Num_ob = 0
nvalids_fm = 0
! 2.1 Loop over stations
! ------------------
stations: &
DO n = 1, number_of_obs
stations_valid: &
IF (obs (n)%info%discard ) THEN
CYCLE stations
ELSE stations_valid
READ (obs (n) % info % platform (4:6), '(I3)') fm
if (fm /= Ob_fm) CYCLE stations
nvalids_fm = nvalids_fm + 1
SELECT CASE (Ob_fm)
! 2.2 count the number of SATOB within the grid box for thining
!
CASE (88)
new = obs (n) % surface % meas
if (eps_equal (new%pressure%data, missing_r, 1.)) then
CYCLE stations
else
! 2.2.1 Calculate the model horizonatl coordinates: I, J, and the
! vertical index K
if (fg_format == 'MM5') then
call LLXY (obs (n) % location % latitude, &
obs (n) % location % longitude, X, Y)
else if (fg_format == 'WRF') then
call latlon_to_ij(map_info, obs(n)%location%latitude, &
obs (n)%location%longitude, X, Y)
X = X + .5
Y = Y + .5
endif
J = INT(X + .5)
I = INT(Y + .5)
K = INT(new % pressure % data / Delta_p)
! 2.2.2 Grid-box counter
Num_ob (i,j,k) = Num_ob (i,j,k) + 1
! 2.2.3 Keep the OBS indices for each of the grid-box
Ob_index(i,j,k, Num_ob(i,j,k)) = n
endif
! 2.3 count the number of SSMI within the grid box for thining
!
CASE (125, 126)
! 2.3.1 Calculate the model horizonatl coordinates: I, J
!
if (fg_format == 'MM5') then
call LLXY (obs (n) % location % latitude, &
obs (n) % location % longitude, X, Y)
else if (fg_format == 'WRF') then
call latlon_to_ij(map_info, obs(n)%location%latitude, &
obs (n)%location%longitude, X, Y)
X = X + .5
Y = Y + .5
endif
J = INT(X + .5)
I = INT(Y + .5)
K = 1
! 2.3.2 Grid-box counter
Num_ob (i,j,k) = Num_ob (i,j,k) + 1
! 2.3.3 Keep the OBS indices for each of the grid-box
Ob_index(i,j,k, Num_ob(i,j,k)) = n
! 2.4 count the number of QSCAT within the grid box for thining
CASE ( 281)
! 2.4.1 Calculate the model horizonatl coordinates: I, J
!
if (fg_format == 'MM5') then
call LLXY (obs (n) % location % latitude, &
obs (n) % location % longitude, X, Y)
else if (fg_format == 'WRF') then
call latlon_to_ij(map_info, obs(n)%location%latitude, &
obs (n)%location%longitude, X, Y)
X = X + .5
Y = Y + .5
endif
J = INT(X + .5)
I = INT(Y + .5)
K = 1
! 2.4.2 Grid-box counter
Num_ob (i,j,k) = Num_ob (i,j,k) + 1
! 2.4.3 Keep the OBS indices for each of the grid-box
Ob_index(i,j,k, Num_ob(i,j,k)) = n
! Others
CASE DEFAULT;
! Should never come here
END SELECT
ENDIF stations_valid
ENDDO stations
! 2.5 Print the total number of boxes and columns data available
! ----------------------------------------------------------
!
WRITE(0,'("VALID NO.=",I6)') nvalids_fm
! 2.5.1 Print the total number of boxes data available
num_box = 0
DO K = 1,LV
kk = 0
max_num = 0
DO I = 1,MIX
DO J = 1,MJX
IF (Num_ob (i,j,k) > 0) THEN
max_num = max(max_num, Num_ob (i,j,k))
kk = kk + 1
num_box = num_box + 1
ENDIF
ENDDO
ENDDO
if (kk >0) &
WRITE(0,'(A,I5,A,I5,A,I3)') "LEVEL ", K, " NUM=", KK, " Max=", max_num
ENDDO
write(0,'(10X,"Total Number of Boxes ",I6)') num_box
! 2.5.2 Print the total number of columes data available
num_box = 0
DO I = 1,MIX
DO J = 1,MJX
kk = 0
DO K = 1,LV
KK = KK + Num_ob (i,j,k)
ENDDO
if (kk > 0) num_box = num_box + 1
enddo
enddo
write(0,'(10X,"Number of columns =", I5)') num_box
! 3.0 Select the ONE data within one grid box
! ---------------------------------------
DO K = 1,LV
DO I = 1,MIX
NEXT_box:DO J = 1,MJX
N_selected = -99
N1_selected = -99
num = Num_ob (i,j,k)
if (num <= 1) cycle NEXT_box
! 3.1 The data nearest to the center of grid box (I,J) is selected
Rmin = 100.0
R1min = 100.0
Select_ob:do nn = 1,num
n = Ob_index(i,j,k,nn)
SELECT CASE (Ob_fm)
! 3.1.1 Select one SATOB within a grid-box
CASE (88)
new = obs (n) % surface % meas
if (new % pressure % qc == 0 .and. &
new % speed % qc == 0 .and. &
new % direction % qc == 0) then
if (fg_format == 'MM5') then
call LLXY (obs (n) % location % latitude, &
obs (n) % location % longitude, X, Y)
else if (fg_format == 'WRF') then
call latlon_to_ij(map_info, obs(n)%location%latitude, &
obs (n)%location%longitude, X, Y)
X = X + .5
Y = Y + .5
endif
X = X + .5
Y = Y + .5
dx = X - float(J)
dy = Y - float(I)
RR = dx*dx + dy*dy
if (RR < Rmin) then
N_selected = n
Rmin = RR
endif
endif
! 3.1.2 Select one good SSMI Retrieval Speed and one PW within a grid-box
! This may be in same obs(n): N_selected = N1_selected = n
!
CASE (125)
if (fg_format == 'MM5') then
call LLXY (obs (n) % location % latitude, &
obs (n) % location % longitude, X, Y)
else if (fg_format == 'WRF') then
call latlon_to_ij(map_info, obs(n)%location%latitude, &
obs (n)%location%longitude, X, Y)
X = X + .5
Y = Y + .5
endif
X = X + .5
Y = Y + .5
dx = X - float(J)
dy = Y - float(I)
RR = dx*dx + dy*dy
! Speed
IF (ASSOCIATED (obs (n) % surface)) then
new = obs (n) % surface % meas
if (new % speed % qc == 0) then
if (RR < Rmin) then
N_selected = n
Rmin = RR
endif
endif
endif
! PW
if (obs(n)%ground% pw % qc == 0) then
if (RR < R1min) then
N1_selected = n
R1min = RR
endif
endif
! 3.1.3 Select one good SSMI Tb within a grid-box
! Only keep the data when the data for all channels are good
CASE (126)
if (obs (n) % ground % tb19v % qc == 0 .and. &
obs (n) % ground % tb19h % qc == 0 .and. &
obs (n) % ground % tb22v % qc == 0 .and. &
obs (n) % ground % tb37v % qc == 0 .and. &
obs (n) % ground % tb37h % qc == 0 .and. &
obs (n) % ground % tb85v % qc == 0 .and. &
obs (n) % ground % tb85h % qc == 0) then
if (fg_format == 'MM5') then
call LLXY (obs (n) % location % latitude, &
obs (n) % location % longitude, X, Y)
else if (fg_format == 'WRF') then
call latlon_to_ij(map_info, obs(n)%location%latitude, &
obs (n)%location%longitude, X, Y)
X = X + .5
Y = Y + .5
endif
X = X + .5
Y = Y + .5
dx = X - float(J)
dy = Y - float(I)
RR = dx*dx + dy*dy
if (RR < Rmin) then
N_selected = n
Rmin = RR
endif
endif
! 3.1.2 Select one good QUIKSCAT within a grid-box
!
CASE (281)
if (fg_format == 'MM5') then
call LLXY (obs (n) % location % latitude, &
obs (n) % location % longitude, X, Y)
else if (fg_format == 'WRF') then
call latlon_to_ij(map_info, obs(n)%location%latitude, &
obs (n)%location%longitude, X, Y)
X = X + .5
Y = Y + .5
endif
X = X + .5
Y = Y + .5
dx = X - float(J)
dy = Y - float(I)
RR = dx*dx + dy*dy
! if the speed is OK, keep this report.
IF (ASSOCIATED (obs (n) % surface)) then
new = obs (n) % surface % meas
if (new % speed % qc == 0) then
if (RR < Rmin) then
N_selected = n
Rmin = RR
endif
endif
endif
! Others
CASE DEFAULT;
! Should never come here
END SELECT
enddo Select_ob
! 3.2 Discard the unselected data
Discard:do nn = 1,num
n = Ob_index(i,j,k,nn)
SELECT CASE (Ob_fm)
! 3.2.1 Thining SATOB data by setting discard = .TRUE.
CASE (88)
if (n /= N_selected) obs(n) % info % discard = .TRUE.
! 3.2.3 Thining SSMI Retrieval data by setting discard = .TRUE.
CASE (125)
if (n /= N_selected .and. &
n /= N1_selected) obs(n) % info % discard = .TRUE.
! 3.2.4 Thining SSMI Tb data by setting discard = .TRUE.
CASE (126)
if (n /= N_selected) obs(n) % info % discard = .TRUE.
! 3.2.5 Thining QUIKSCAT data by setting discard = .TRUE.
CASE (281)
if (n /= N_selected) obs(n) % info % discard = .TRUE.
! Others
CASE DEFAULT;
! Should never come here
END SELECT
enddo Discard
ENDDO NEXT_box
ENDDO
ENDDO
! 4.0 deallocate the working arrays
deallocate (Num_ob)
deallocate (Ob_index)
END SUBROUTINE THIN_OB
END MODULE module_thin_ob