var/da/da_biascorr_airmass/rad_bias.f90

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