var/da/da_biascorr_airmass/da_bias_verif.f90

   1   program da_bias_verif
   2
   3   use rad_bias, only : bias, long, jpchan, jpscan,jband, jpnx, &
   4      da_read_biasprep, read_biascoef, get_scorr, qc_amsua, qc_amsub, mask, &
   5      write_biascoef, jpny, qc_ssmis
   6
   7 ! PURPOSE.
   8 ! --------
   9 ! Apply bias correction to radiance data
  10 ! to statistically verify algorithm
  11
  12   implicit none
  13
  14   type (bias) :: tovs
  15
  16   real (kind=long) :: vmean_dep(JPCHAN), vmean_abs(JPCHAN)
  17   real (kind=long) :: vec_dep(JPCHAN), vec_abs(JPCHAN)
  18   real (kind=long) :: vstd_dep(JPCHAN), vstd_abs(JPCHAN)
  19   real (kind=long) :: vmean_dep1(JPCHAN), vmean_abs1(JPCHAN)
  20   real (kind=long) :: vstd_dep1(JPCHAN), vstd_abs1(JPCHAN)
  21
  22   real (kind=long) :: airbias(JPCHAN)
  23
  24   real (kind=long) :: xcoef0(JPCHAN), xcoef(JPCHAN,JPNX)
  25
  26   real (kind=long) :: vmnrl(JPCHAN,JPSCAN) = 0.0
  27   real (kind=long) :: vmnrlb(JPCHAN,JPSCAN,JBAND) = 0.0
  28
  29   real (kind=long) :: SCORR(JPCHAN)
  30
  31   integer :: nsel(10)
  32   integer :: nobs(JPCHAN),nobs1(JPCHAN)
  33
  34   logical :: LMASK
  35
  36   logical :: lscan = .FALSE.
  37
  38   integer :: I, iband, ierr
  39   integer :: J, JSCAN, jv
  40
  41   integer :: kscanx=90
  42   ! integer :: jbandx
  43   logical :: check_limb=.false., check_mask=.false., global
  44   real    :: FAC = 3.0      ! Number of SD' for QC
  45
  46   integer :: nchan,nscan,nband,npred
  47
  48   namelist /INPUTS/ global,lscan, kscanx, check_limb, check_mask, FAC
  49
  50 !------------------------------------------------------------------------------
  51 !        1.   SETUP.
  52 !             -----
  53   read (5,INPUTS,end=100)
  54   100 continue
  55   write (6,INPUTS)
  56
  57 ! 1.0 read bias correction coefficients
  58 !   -------------------------------------------------------------
  59      open(unit=12,file='scor.asc',form='formatted')
  60      read (12,'(4i6)') nchan,nscan,nband,npred
  61      close(12)
  62
  63      call read_biascoef(nchan,nscan,nband,npred,global, &
  64                       VMNRL(1:nchan,1:nscan),VMNRLB(1:nchan,1:nscan,1:nband), &
  65                       xcoef(1:nchan,1:npred),xcoef0(1:nchan), &
  66                       nobs(1:nchan),vmean_abs(1:nchan), &
  67                       vstd_abs(1:nchan),vmean_dep(1:nchan), &
  68                       vstd_dep(1:nchan) )
  69
  70 !----------------------------------------------------------------------
  71 !        2.   READ IN DATA, Q.C., CALC MEANS AND VARIANCES.
  72 !             ---- -- ----- ----- ---- ----- --- ---------
  73
  74   nobs1(:)  = 0
  75   nsel(:)  = 0
  76   vmean_dep1(:) = 0.0
  77   vmean_abs1(:) = 0.0
  78   vstd_dep1(:)  = 0.0
  79   vstd_abs1(:)  = 0.0
  80
  81   read (unit=10,end=300)  tovs%nchan, tovs%npred    ! Read in data
  82   rewind(unit=10)
  83
  84   allocate(tovs%tb(tovs%nchan))
  85   allocate(tovs%omb(tovs%nchan))
  86   allocate(tovs%bias(tovs%nchan))
  87   allocate(tovs%qc_flag(tovs%nchan))
  88   allocate(tovs%cloud_flag(tovs%nchan))
  89   allocate(tovs%pred(tovs%npred))
  90
  91   300 continue
  92
  93 loop2:&
  94   do
  95     call da_read_biasprep(tovs,10,ierr)
  96     if (ierr == 0) then      ! not end
  97          continue
  98     elseif (ierr == 1) then  ! end
  99          exit
 100     else                     ! error
 101          stop 'read error in da_veri_bias'
 102     end if
 103
 104 ! latitude band
 105     iband = INT(tovs%lat/30.000001) + 3
 106     if (lscan) then
 107         JSCAN = tovs%scanpos
 108       if (global) then
 109         call GET_SCORR(JPCHAN,SCORR(1:JPCHAN),tovs%lat,vmnrlb,JSCAN)
 110       else
 111         SCORR(1:nchan) = vmnrl(1:nchan,JSCAN)
 112       end if
 113     ELSE
 114         JSCAN = KSCANX/2
 115         SCORR(1:JPCHAN) = 0.0
 116     end if
 117
 118     vec_dep(1:nchan) = tovs%omb(1:nchan) - SCORR(1:nchan)
 119     vec_abs(1:nchan) =  tovs%tb(1:nchan) - SCORR(1:nchan)
 120
 121 ! 3.2 QC:
 122       if (tovs%sensor_id == 3) then
 123            call qc_amsua(tovs)
 124       elseif(tovs%sensor_id == 4) then
 125            call qc_amsub(tovs)
 126       elseif(tovs%sensor_id == 15) then  ! mhs
 127            call qc_amsub(tovs)
 128       elseif(tovs%sensor_id == 10) then  ! ssmis
 129            call qc_ssmis(tovs)
 130       end if
 131
 132 ! 3.3 limb scan check
 133 !---------------------
 134    if (check_limb) then
 135     if ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCANX-1))) then
 136       nsel(2) = nsel(2) + 1
 137       CYCLE loop2
 138     end if
 139    end if
 140
 141 ! 3.4 Reject data outside radiosonde mask
 142 !------------------------------------------
 143     if (check_mask) then
 144       call MASK(tovs%lat,tovs%lon,LMASK)
 145       if (.NOT. LMASK) then
 146         nsel(8) = nsel(8) + 1
 147         CYCLE loop2
 148       end if
 149     end if
 150
 151 ! 3.5 Reject outliers : facx*sigma,
 152 !--------------------------------------------------------------------------
 153     do j=1, tovs%nchan
 154       if ( (abs(vec_dep(j)-vmean_dep(j)) > (vstd_dep(j)*FAC)) ) then
 155         tovs%qc_flag(j) = -1
 156       end if
 157     end do
 158
 159     do i=1,nchan
 160       airbias(i) = xcoef0(i)
 161       do j=1,npred
 162         airbias(i) = airbias(i) + tovs%pred(j)*xcoef(i,j)
 163       end do
 164       vec_dep(i) = vec_dep(i) - airbias(i)
 165       vec_abs(i) = vec_abs(i) - airbias(i)
 166       if ( tovs%omb(i) == -888888.00 ) vec_dep(i)=tovs%omb(i)
 167     end do
 168
 169 ! mean/std statistics for scan/airmass-bias corrected values
 170     do j=1, nchan
 171       if ( tovs%qc_flag(j) == 1 ) then
 172         jv = j
 173         nobs1(j) = nobs1(j) + 1
 174         vmean_dep1(j) = vmean_dep1(j) + vec_dep(j)
 175          vstd_dep1(j) = vstd_dep1(j) + vec_dep(j)*vec_dep(j)
 176         vmean_abs1(j) = vmean_abs1(j) + vec_abs(j)
 177          vstd_abs1(j) = vstd_abs1(j) + vec_abs(j)*vec_abs(j)
 178       end if
 179     end do
 180
 181     write (11,'(30i15)') tovs%qc_flag(1:nchan)
 182     write (11,'(30f15.3)') tovs%omb(1:nchan)  ! omb no-biascorrection
 183     write (11,'(30f15.3)') vec_dep(1:nchan)   ! omb with bias correction
 184
 185   end do loop2
 186
 187 ! Calculate means, standard deviations and covariances
 188
 189   where (nobs1(:) /= 0)
 190     vmean_dep1(:) = vmean_dep1(:)/nobs1(:)
 191     vstd_dep1(:)  = vstd_dep1(:)/nobs1(:) - vmean_dep1(:)**2
 192     vstd_dep1(:)  = SQRT(MAX(0.0_8,vstd_dep1(:)))
 193     vmean_abs1(:) = vmean_abs1(:)/nobs1(:)
 194     vstd_abs1(:)  = vstd_abs1(:)/nobs1(:) - vmean_abs1(:)**2
 195     vstd_abs1(:)  = SQRT(MAX(0.0_8,vstd_abs1(:)))
 196   end where
 197
 198 !----------------------------------------------------------------------------
 199
 200    deallocate(tovs%tb)
 201    deallocate(tovs%omb)
 202    deallocate(tovs%qc_flag)
 203    deallocate(tovs%cloud_flag)
 204    deallocate(tovs%pred)
 205
 206 ! out coefs to ASCII file bcor.asc
 207   call write_biascoef(nchan,nscan,nband,npred,global, &
 208                       VMNRL(1:nchan,1:nscan),VMNRLB(1:nchan,1:nscan,1:nband), &
 209                       xcoef(1:nchan,1:npred),xcoef0(1:nchan), &
 210                       nobs1(1:nchan),vmean_abs1(1:nchan), &
 211                       vstd_abs1(1:nchan),vmean_dep1(1:nchan), &
 212                       vstd_dep1(1:nchan) )
 213
 214 end program da_bias_verif