var/obsproc/src/module_thin_ob.F90

   1 MODULE module_thin_ob
   2
   3   USE module_type
   4   USE module_func
   5   USE module_mm5
   6   USE module_map
   7   USE map_utils
   8   USE module_namelist
   9
  10 CONTAINS
  11
  12 SUBROUTINE THIN_OB (max_number_of_obs, obs, number_of_obs, mix, mjx, &
  13                                  missing_r, Ob_name, Ob_fm, Delta_P)
  14
  15 ! --------------------------------------------------------------------------
  16 ! Purpose: To thin the sateliite (SATOB, SSMI Retrieval and Tb) data.
  17 ! Method : To select only ONE data nearest to the center of the (I,J,K)
  18 !          grid-box if there are more than one OBS within that gridbox.
  19 !
  20 !          The vertical P_thickness is defined by Delta_P.
  21 !          Currently for SATOB the maximum of 100 level in vertical
  22 !          and maximum of 20 data points within a grid box are allowed.
  23 !          For the single level SSMI the maximum  of 400 data points
  24 !          within a gridbox is allowed.
  25 !
  26 !                                          Yong-Run Guo 05/01/2001
  27 ! --------------------------------------------------------------------------
  28 !
  29
  30   implicit none
  31
  32   INTEGER,                                  INTENT (in) :: max_number_of_obs
  33   TYPE (report), DIMENSION (max_number_of_obs), INTENT (inout) :: obs
  34   INTEGER,                                      INTENT (in) :: number_of_obs
  35   CHARACTER (LEN = *),                          INTENT (in) :: Ob_name
  36   INTEGER,                                      INTENT (in) :: Ob_fm
  37   REAL,                                         INTENT (in) :: missing_r
  38   REAL   ,                     OPTIONAL,        INTENT (in) :: Delta_P
  39
  40   INTEGER                       :: n, fm, i, j, k, kk, num_box, &
  41                                    max_num, mix, mjx, num,  nn, total, &
  42                                    max_box, LV, N_selected, N1_selected
  43   INTEGER, ALLOCATABLE          :: Num_ob(:,:,:), Ob_index(:,:,:,:)
  44   REAL                          :: X, Y, dx, dy, RR, Rmin, R1min
  45   TYPE ( meas_data )            :: new
  46   INTEGER                       :: nvalids_fm
  47
  48 !------------------------------------------------------------------
  49 !
  50 ! 1.0 Allocate the working arrays
  51 !
  52 !     The default size: LV, max_box can be changed to meet the needs.
  53 !
  54
  55    ! For SATOB
  56
  57        if (Ob_fm == 88) then
  58          LV      = 100
  59          max_box =  50
  60
  61    ! For SSMI
  62
  63        else if (Ob_fm == 125 .or. Ob_fm == 126 .or. Ob_fm == 281) then
  64          LV      =   1
  65          max_box = 400
  66
  67        else
  68          WRITE(0,'(A,I3,A)') ' FM=',Ob_fm,' NOT SATELLITE OBS,', &
  69                              ' NO THINING APPLIED TO IT.'
  70          return
  71        endif
  72
  73        allocate (Num_ob  (mix,mjx,LV))
  74        allocate (Ob_index(mix,mjx,LV,max_box))
  75
  76        WRITE(0,'(/"Thining OBS ==> ",A,1X,"FM=",I3,2X,      &
  77        &                   "mix, mjx, LV, max_box:",4I6)') &
  78        &    Ob_name, Ob_fm, mix, mjx, LV, max_box
  79 !
  80 ! 2.0  To count the number of OBS within the grid boxes
  81 !      ------------------------------------------------
  82
  83        Num_ob = 0
  84        nvalids_fm = 0
  85
  86 ! 2.1 Loop over stations
  87 !     ------------------
  88
  89 stations: &
  90        DO n = 1, number_of_obs
  91
  92 stations_valid: &
  93          IF (obs (n)%info%discard ) THEN
  94
  95            CYCLE  stations
  96
  97          ELSE stations_valid
  98
  99            READ (obs (n) % info % platform (4:6), '(I3)') fm
 100
 101            if (fm /= Ob_fm)  CYCLE stations
 102
 103            nvalids_fm = nvalids_fm + 1
 104
 105        SELECT CASE (Ob_fm)
 106
 107 ! 2.2 count the number of SATOB within the grid box for thining
 108 !
 109        CASE (88)
 110
 111            new = obs (n) % surface % meas
 112
 113            if (eps_equal (new%pressure%data, missing_r, 1.)) then
 114               CYCLE stations
 115            else
 116
 117 ! 2.2.1 Calculate the model horizonatl coordinates: I, J, and the
 118 !       vertical index K
 119
 120            if (fg_format == 'MM5') then
 121              call LLXY (obs (n) % location % latitude, &
 122                         obs (n) % location % longitude, X, Y)
 123            else if (fg_format == 'WRF') then
 124              call latlon_to_ij(map_info, obs(n)%location%latitude, &
 125                                obs (n)%location%longitude, X, Y)
 126              X = X + .5
 127              Y = Y + .5
 128            endif
 129              J = INT(X + .5)
 130              I = INT(Y + .5)
 131              K = INT(new % pressure % data / Delta_p)
 132
 133 ! 2.2.2 Grid-box counter
 134
 135              Num_ob (i,j,k) = Num_ob (i,j,k) + 1
 136
 137 ! 2.2.3 Keep the OBS indices for each of the grid-box
 138
 139              Ob_index(i,j,k, Num_ob(i,j,k)) = n
 140
 141            endif
 142
 143 ! 2.3 count the number of SSMI within the grid box for thining
 144 !
 145
 146        CASE (125, 126)
 147
 148 ! 2.3.1 Calculate the model horizonatl coordinates: I, J
 149 !
 150            if (fg_format == 'MM5') then
 151              call LLXY (obs (n) % location % latitude, &
 152                         obs (n) % location % longitude, X, Y)
 153            else if (fg_format == 'WRF') then
 154              call latlon_to_ij(map_info, obs(n)%location%latitude, &
 155                                obs (n)%location%longitude, X, Y)
 156              X = X + .5
 157              Y = Y + .5
 158            endif
 159              J = INT(X + .5)
 160              I = INT(Y + .5)
 161              K = 1
 162 ! 2.3.2 Grid-box counter
 163
 164              Num_ob (i,j,k) = Num_ob (i,j,k) + 1
 165
 166 ! 2.3.3 Keep the OBS indices for each of the grid-box
 167
 168              Ob_index(i,j,k, Num_ob(i,j,k)) = n
 169
 170 ! 2.4 count the number of QSCAT within the grid box for thining
 171
 172        CASE ( 281)
 173
 174 ! 2.4.1 Calculate the model horizonatl coordinates: I, J
 175 !
 176            if (fg_format == 'MM5') then
 177              call LLXY (obs (n) % location % latitude, &
 178                         obs (n) % location % longitude, X, Y)
 179            else if (fg_format == 'WRF') then
 180              call latlon_to_ij(map_info, obs(n)%location%latitude, &
 181                                obs (n)%location%longitude, X, Y)
 182              X = X + .5
 183              Y = Y + .5
 184            endif
 185              J = INT(X + .5)
 186              I = INT(Y + .5)
 187              K = 1
 188 ! 2.4.2 Grid-box counter
 189
 190              Num_ob (i,j,k) = Num_ob (i,j,k) + 1
 191
 192 ! 2.4.3 Keep the OBS indices for each of the grid-box
 193
 194              Ob_index(i,j,k, Num_ob(i,j,k)) = n
 195
 196        !  Others
 197
 198        CASE DEFAULT;
 199
 200            !  Should never come here
 201
 202        END SELECT
 203          ENDIF stations_valid
 204
 205       ENDDO  stations
 206
 207 ! 2.5 Print the total number of boxes and columns data available
 208 !     ----------------------------------------------------------
 209 !
 210       WRITE(0,'("VALID NO.=",I6)') nvalids_fm
 211
 212 ! 2.5.1 Print the total number of boxes data available
 213
 214       num_box = 0
 215       DO K = 1,LV
 216       kk = 0
 217       max_num = 0
 218         DO I = 1,MIX
 219         DO J = 1,MJX
 220           IF (Num_ob (i,j,k) > 0) THEN
 221             max_num = max(max_num, Num_ob (i,j,k))
 222             kk = kk + 1
 223             num_box = num_box + 1
 224           ENDIF
 225         ENDDO
 226         ENDDO
 227         if (kk >0) &
 228         WRITE(0,'(A,I5,A,I5,A,I3)') "LEVEL ", K, " NUM=", KK, " Max=", max_num
 229       ENDDO
 230
 231       write(0,'(10X,"Total Number of Boxes ",I6)') num_box
 232
 233 ! 2.5.2 Print the total number of columes data available
 234
 235         num_box = 0
 236         DO I = 1,MIX
 237         DO J = 1,MJX
 238           kk = 0
 239           DO K = 1,LV
 240             KK = KK + Num_ob (i,j,k)
 241           ENDDO
 242           if (kk > 0) num_box = num_box + 1
 243         enddo
 244         enddo
 245         write(0,'(10X,"Number of columns =", I5)') num_box
 246
 247 ! 3.0 Select the ONE data within one grid box
 248 !     ---------------------------------------
 249
 250       DO K = 1,LV
 251
 252           DO I = 1,MIX
 253
 254   NEXT_box:DO J = 1,MJX
 255
 256           N_selected  = -99
 257           N1_selected = -99
 258
 259           num = Num_ob (i,j,k)
 260           if (num <= 1) cycle NEXT_box
 261
 262 ! 3.1 The data nearest to the center of grid box (I,J) is selected
 263
 264           Rmin  = 100.0
 265           R1min = 100.0
 266
 267    Select_ob:do nn = 1,num
 268             n = Ob_index(i,j,k,nn)
 269
 270        SELECT CASE (Ob_fm)
 271
 272 ! 3.1.1 Select one SATOB within a grid-box
 273
 274        CASE (88)
 275
 276             new = obs (n) % surface % meas
 277             if (new % pressure  % qc == 0 .and. &
 278                 new % speed     % qc == 0 .and. &
 279                 new % direction % qc == 0) then
 280                 if (fg_format == 'MM5') then
 281                   call LLXY (obs (n) % location % latitude, &
 282                              obs (n) % location % longitude, X, Y)
 283                 else if (fg_format == 'WRF') then
 284                   call latlon_to_ij(map_info, obs(n)%location%latitude, &
 285                                     obs (n)%location%longitude, X, Y)
 286                   X = X + .5
 287                   Y = Y + .5
 288                 endif
 289                 X = X + .5
 290                 Y = Y + .5
 291
 292                 dx = X - float(J)
 293                 dy = Y - float(I)
 294                 RR = dx*dx + dy*dy
 295                 if (RR < Rmin) then
 296                   N_selected = n
 297                   Rmin = RR
 298                 endif
 299             endif
 300
 301 ! 3.1.2 Select one good SSMI Retrieval Speed and one PW within a grid-box
 302 !       This may be in same obs(n): N_selected = N1_selected = n
 303 !
 304        CASE (125)
 305                 if (fg_format == 'MM5') then
 306                   call LLXY (obs (n) % location % latitude, &
 307                              obs (n) % location % longitude, X, Y)
 308                 else if (fg_format == 'WRF') then
 309                   call latlon_to_ij(map_info, obs(n)%location%latitude, &
 310                                     obs (n)%location%longitude, X, Y)
 311                   X = X + .5
 312                   Y = Y + .5
 313                 endif
 314              X = X + .5
 315              Y = Y + .5
 316
 317              dx = X - float(J)
 318              dy = Y - float(I)
 319              RR = dx*dx + dy*dy
 320
 321     ! Speed
 322              IF (ASSOCIATED (obs (n) % surface)) then
 323                 new = obs (n) % surface % meas
 324                 if (new % speed     % qc == 0) then
 325                    if (RR < Rmin) then
 326                       N_selected = n
 327                       Rmin = RR
 328                    endif
 329                 endif
 330              endif
 331     ! PW
 332              if (obs(n)%ground% pw  % qc == 0) then
 333                 if (RR < R1min) then
 334                    N1_selected = n
 335                    R1min = RR
 336                 endif
 337              endif
 338
 339 ! 3.1.3 Select one good SSMI Tb within a grid-box
 340 !       Only keep the data when the data for all channels are good
 341
 342        CASE (126)
 343
 344             if (obs (n) % ground % tb19v % qc == 0 .and. &
 345                 obs (n) % ground % tb19h % qc == 0 .and. &
 346                 obs (n) % ground % tb22v % qc == 0 .and. &
 347                 obs (n) % ground % tb37v % qc == 0 .and. &
 348                 obs (n) % ground % tb37h % qc == 0 .and. &
 349                 obs (n) % ground % tb85v % qc == 0 .and. &
 350                 obs (n) % ground % tb85h % qc == 0) then
 351                 if (fg_format == 'MM5') then
 352                   call LLXY (obs (n) % location % latitude, &
 353                              obs (n) % location % longitude, X, Y)
 354                 else if (fg_format == 'WRF') then
 355                   call latlon_to_ij(map_info, obs(n)%location%latitude, &
 356                                     obs (n)%location%longitude, X, Y)
 357                   X = X + .5
 358                   Y = Y + .5
 359                 endif
 360                 X = X + .5
 361                 Y = Y + .5
 362
 363                 dx = X - float(J)
 364                 dy = Y - float(I)
 365                 RR = dx*dx + dy*dy
 366                 if (RR < Rmin) then
 367                   N_selected = n
 368                   Rmin = RR
 369                 endif
 370              endif
 371
 372 ! 3.1.2 Select one good QUIKSCAT within a grid-box
 373 !
 374        CASE (281)
 375                 if (fg_format == 'MM5') then
 376                   call LLXY (obs (n) % location % latitude, &
 377                              obs (n) % location % longitude, X, Y)
 378                 else if (fg_format == 'WRF') then
 379                   call latlon_to_ij(map_info, obs(n)%location%latitude, &
 380                                     obs (n)%location%longitude, X, Y)
 381                   X = X + .5
 382                   Y = Y + .5
 383                 endif
 384                 X = X + .5
 385                 Y = Y + .5
 386
 387              dx = X - float(J)
 388              dy = Y - float(I)
 389              RR = dx*dx + dy*dy
 390
 391     ! if the speed is OK, keep this report.
 392
 393              IF (ASSOCIATED (obs (n) % surface)) then
 394                 new = obs (n) % surface % meas
 395                 if (new % speed     % qc == 0) then
 396                    if (RR < Rmin) then
 397                       N_selected = n
 398                       Rmin = RR
 399                    endif
 400                 endif
 401              endif
 402
 403        !  Others
 404
 405        CASE DEFAULT;
 406
 407            !  Should never come here
 408
 409        END SELECT
 410
 411           enddo Select_ob
 412
 413 ! 3.2 Discard the unselected data
 414
 415    Discard:do nn = 1,num
 416             n = Ob_index(i,j,k,nn)
 417
 418        SELECT CASE (Ob_fm)
 419
 420 ! 3.2.1 Thining SATOB data by setting discard = .TRUE.
 421
 422        CASE (88)
 423
 424             if (n /= N_selected) obs(n) % info % discard = .TRUE.
 425
 426 ! 3.2.3 Thining SSMI Retrieval data by setting discard = .TRUE.
 427
 428        CASE (125)
 429
 430             if (n /= N_selected .and. &
 431                 n /= N1_selected) obs(n) % info % discard = .TRUE.
 432
 433
 434 ! 3.2.4 Thining SSMI Tb data by setting discard = .TRUE.
 435
 436        CASE (126)
 437
 438             if (n /= N_selected) obs(n) % info % discard = .TRUE.
 439
 440 ! 3.2.5 Thining QUIKSCAT data by setting discard = .TRUE.
 441
 442        CASE (281)
 443
 444             if (n /= N_selected) obs(n) % info % discard = .TRUE.
 445
 446        !  Others
 447
 448        CASE DEFAULT;
 449
 450            !  Should never come here
 451
 452        END SELECT
 453
 454           enddo Discard
 455
 456      ENDDO NEXT_box
 457
 458      ENDDO
 459
 460      ENDDO
 461
 462 ! 4.0 deallocate the working arrays
 463
 464      deallocate (Num_ob)
 465      deallocate (Ob_index)
 466
 467 END SUBROUTINE THIN_OB
 468
 469 END MODULE module_thin_ob