var/da/da_biascorr_airmass/da_bias_airmass.f90

   1 program da_bias_airmass
   2
   3    use rad_bias
   4
   5   ! PURPOSE.
   6   ! --------
   7   ! To calculate spatially varying bias correction COEFFICIENTS
   8   ! for TOVS using model bias predictors.
   9
  10   ! EXTERNALS.
  11   ! ----------
  12   ! REGRESS_ONE
  13
  14   ! REFERENCE.
  15   ! ----------
  16   ! Harris and Kelly 1998.
  17
  18   implicit none
  19
  20   type (bias) :: tovs
  21
  22   real (kind=long) :: vmean_dep(JPCHAN), vmean_abs(JPCHAN)
  23   real (kind=long) :: vec_dep(JPCHAN), vec_abs(JPCHAN)
  24   real (kind=long) :: vstd_dep(JPCHAN), vstd_abs(JPCHAN)
  25   real (kind=long) :: vmn(JPCHAN), vstd(JPCHAN)
  26
  27   real :: pred(JPNX)
  28
  29   real (kind=long) :: xbar(JPNY,JPNX), ybar(JPNY)
  30   real (kind=long) :: xcov(JPNY,JPNX,JPNX), ycov(JPNY), xycov(JPNY,JPNX)
  31
  32   real (kind=long) :: coef(JPNY,JPNX), coef0(JPNY)
  33   real (kind=long) :: reserr(JPNY), rescov(JPNY,JPNY)
  34   real (kind=long) :: xvec(JPNY,JPNX,JPNX)
  35   real (kind=long) :: xreser(JPCHAN), xcoef0(JPCHAN), xcoef(JPCHAN,JPNX)
  36
  37   real (kind=long) :: vmnbd(JPNY,6), vstdbd(JPNY,6)
  38   real (kind=long) :: vmnbdl(JPNY,6,0:1), vstdbdl(JPNY,6,0:1)
  39   real (kind=long) :: vmnb0(JPNY,6), vstdb0(JPNY,6)
  40   real (kind=long) :: vmnbl0(JPNY,6,0:1), vstdbl0(JPNY,6,0:1)
  41   real (kind=long) :: dvmnbd(JPNY,6), dvstdb(JPNY,6)
  42
  43   real (kind=long) :: vmnrl(JPCHAN,JPSCAN) = 0.0
  44   real (kind=long) :: vmnrlb(JPCHAN,JPSCAN,JBAND) = 0.0
  45
  46   real (kind=long) :: SCORR(JPCHAN)
  47
  48   integer :: nsel(10)
  49   integer :: nobs(JPCHAN)
  50   integer :: nobsy(JPNY),lchan
  51
  52   integer :: nband(6), nbandl(6,0:1), n_band_ch(6,JPNY,0:1) = 0
  53   integer :: nband0(JPNY,6), nbandd(JPNY,6)
  54   integer :: nbandl0(JPNY,6,0:1), nbanddl(JPNY,6,0:1)
  55
  56   logical :: LMASK
  57
  58   logical :: lscan = .FALSE., global = .FALSE.
  59
  60   integer :: IR, ibin, I, iband, II, IV, ib, ierr
  61   integer :: JS, J, JJ, JCOUNT, JBOX, JMINI, JSCAN, jv, IIV, JJV, sband
  62
  63   real (kind=long) :: xcorr(JPCHAN)
  64   real (kind=long) :: coef_sat, coef_year, coef_month, coef_day, coef_time
  65   integer :: CDATE = 1998010105                           ! Current Date
  66
  67   integer :: icount = 0
  68   integer :: KSCAN = JPSCAN
  69   logical :: check_limb=.false., check_mask=.false.
  70   real    :: FAC = 3.0      ! Number of SD' for QC
  71
  72
  73   namelist /INPUTS/ global, lscan, kscan, check_limb, check_mask, &
  74                     FAC,CDATE
  75
  76 !------------------------------------------------------------------------------
  77 !        1.   SETUP.
  78 !             -----
  79
  80   read (5,INPUTS,end=100)
  81   100 continue
  82   write (6,INPUTS)
  83
  84   if (lscan) OPEN(12,FORM='UNformatTED')  ! Scan Biases with lat bands
  85   OPEN(unit=10,FORM='UNformatTED')        ! Input data from SELECT
  86
  87 ! Read scan biases
  88
  89   if (lscan) then
  90      write (6,112) (JS,JS=1,KSCAN)
  91      112 format(/1X,'SCAN DEPendENT BIASES APPLIED'/1X,3X,18I7/(4X,18I7))
  92   ELSE
  93      write (6,*) 'NO SCAN CORRECTION APPLIED'
  94   end if
  95
  96 !----------------------------------------------------------------------
  97 !        2.   READ IN DATA, Q.C., CALC MEANS AND VARIANCES.
  98 !             ---- -- ----- ----- ---- ----- --- ---------
  99
 100   nobs(:)  = 0
 101   nsel(:)  = 0
 102   vmean_dep(:) = 0.0
 103   vmean_abs(:) = 0.0
 104   vstd_dep(:)  = 0.0
 105   vstd_abs(:)  = 0.0
 106
 107   200  continue
 108
 109   read (unit=10,end=265)  tovs%nchan, tovs%npred    ! Read in data
 110   rewind(unit=10)
 111
 112   allocate(tovs%tb(tovs%nchan))
 113   allocate(tovs%omb(tovs%nchan))
 114   allocate(tovs%bias(tovs%nchan))
 115   allocate(tovs%qc_flag(tovs%nchan))
 116   allocate(tovs%cloud_flag(tovs%nchan))
 117   allocate(tovs%pred(tovs%npred))
 118
 119   if (lscan) then
 120    if (global) then
 121     do J=1, tovs%nchan   ! get scan biases (latitude bands)
 122     do sband=1,JBAND
 123         read (12) JJ, vmnrlb(J,1:KSCAN,sband)
 124     end do
 125     end do
 126    else
 127     do J=1, tovs%nchan   ! get scan biases (no band dependence)
 128         read (12) JJ, vmnrl(J,1:KSCAN)
 129     end do
 130    end if
 131   end if
 132
 133 loop1:&
 134   do
 135     icount=icount+1
 136
 137     call da_read_biasprep(tovs,10,ierr)
 138     if (ierr == 0) then      ! not end
 139          continue
 140     elseif (ierr == 1) then  ! end
 141          exit
 142     else                     ! error
 143          stop 'read error in da_airmass_bias'
 144     end if
 145
 146 ! select lat band and get scan bias
 147       iband = INT(tovs%lat/30.000001) + 3
 148       if (lscan) then
 149         JSCAN = tovs%scanpos
 150        if (global) then
 151         call GET_SCORR(JPCHAN,SCORR(1:JPCHAN),tovs%lat,vmnrlb,JSCAN)
 152        else
 153         SCORR(1:tovs%nchan) = vmnrl(1:tovs%nchan,JSCAN)
 154        end if
 155       ELSE
 156         JSCAN = KSCAN/2
 157         SCORR(1:JPCHAN) = 0.0
 158       end if
 159
 160 ! apply scan bias to the departure (from FG or analysis) and the TB
 161       vec_dep(1:tovs%nchan) = tovs%omb(1:tovs%nchan) - SCORR(1:tovs%nchan)
 162       vec_abs(1:tovs%nchan) =  tovs%tb(1:tovs%nchan) - SCORR(1:tovs%nchan)
 163
 164 !
 165 !   2.2 QC: extrme values/extrme departure/window channel
 166 !
 167       if (tovs%sensor_id == 3) then
 168            call qc_amsua(tovs)
 169       elseif(tovs%sensor_id == 4) then
 170            call qc_amsub(tovs)
 171       elseif(tovs%sensor_id == 15) then ! mhs
 172            call qc_amsub(tovs)
 173       elseif(tovs%sensor_id == 10) then ! ssmis
 174            call qc_ssmis(tovs)
 175       end if
 176
 177 !-------------------------
 178 ! 2.3 limb sounding check
 179 !-------------------------
 180   if (check_limb) then
 181     if ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCAN-1))) then  ! Reject edge of scan
 182       nsel(2) = nsel(2) + 1
 183       CYCLE loop1
 184     end if
 185   end if
 186
 187 !-----------------------------------------
 188 ! 2.4 Reject data outside radiosonde mask
 189 !-----------------------------------------
 190     if (check_mask) then
 191       call MASK(tovs%lat,tovs%lon,LMASK)
 192       if (.NOT. LMASK) then
 193         nsel(8) = nsel(8) + 1
 194         CYCLE loop1
 195       end if
 196     end if
 197
 198 ! Good data : count and use in the statistics
 199
 200     do j=1, tovs%nchan
 201       if ( tovs%qc_flag(j) == 1 ) then
 202 !  statistics for channel
 203 !--------------------------
 204         nobs(j) = nobs(j) + 1
 205         vmean_dep(j) = vmean_dep(j) + vec_dep(j)
 206          vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
 207         vmean_abs(j) = vmean_abs(j) + vec_abs(j)
 208          vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
 209       end if
 210     end do
 211
 212   end do loop1
 213
 214   265 continue
 215
 216 !---------------------------------
 217 !  2.8 mean and std for channels
 218 !---------------------------------
 219   where (nobs(:) /= 0)
 220     vmean_dep(:) = vmean_dep(:)/nobs(:)
 221     vstd_dep(:)  = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
 222     vstd_dep(:)  = SQRT(MAX(0.0_8,vstd_dep(:)))
 223     vmean_abs(:) = vmean_abs(:)/nobs(:)
 224     vstd_abs(:)  = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
 225     vstd_abs(:)  = SQRT(MAX(0.0_8,vstd_abs(:)))
 226   end where
 227
 228   write (6,270) icount,nobs(1:tovs%nchan)
 229   270 format (/1X,'NUMBER OF DATA Total/ACCEPTED'/1X,i10/1X,15I10)
 230
 231   write (6,288)
 232   288 format (/1X,'FIRST PASS: MEANS AND STANDARD DEVIATIONS')
 233
 234   do j=1, tovs%nchan
 235     jv = j
 236     write (6,289) jv, nobs(j), vmean_abs(j), vstd_abs(j), vmean_dep(j), vstd_dep(j)
 237     289 format (1X,I5,I10,4F15.2)
 238   end do
 239
 240 !-----------------------------------------------------------------------------
 241 !     3.   SECOND PASS THROUGH DATA, EXTRA Q.C. (sigma-elimination)
 242 !          ------ ---- ------- ----- ----- ----
 243
 244   300 continue
 245
 246   vmn(:) = vmean_dep(:)
 247   vstd(:)= vstd_dep(:)
 248
 249   rewind(unit=10)
 250
 251   vmean_dep(:) = 0.0                                    ! Clear matrices
 252   vmean_abs(:) = 0.0                                    ! Clear matrices
 253   nobs(:) = 0
 254   nsel(:) = 0
 255   nband(:) = 0
 256   nbandl(:,:) = 0
 257   nobsy(:) = 0
 258
 259   xbar(:,:)   = 0.0
 260   ybar(:)     = 0.0
 261   xcov(:,:,:) = 0.0
 262   ycov(:)     = 0.0
 263   xycov(:,:)  = 0.0
 264
 265   icount = 0
 266 loop2:&
 267   do
 268     icount = icount + 1
 269
 270     call da_read_biasprep(tovs,10,ierr)
 271     if (ierr == 0) then      ! not end
 272          continue
 273     elseif (ierr == 1) then  ! end
 274          exit
 275     else                     ! error
 276          stop 'read error in da_airmass_bias'
 277     end if
 278
 279 ! latitude band
 280     iband = INT(tovs%lat/30.000001) + 3
 281     if (lscan) then
 282        JSCAN = tovs%scanpos
 283        if (global) then
 284         call GET_SCORR(JPCHAN,SCORR(1:JPCHAN),tovs%lat,vmnrlb,JSCAN)
 285        else
 286         SCORR(1:tovs%nchan) = vmnrl(1:tovs%nchan,JSCAN)
 287        end if
 288     ELSE
 289        JSCAN = KSCAN/2
 290        SCORR(1:JPCHAN) = 0.0
 291     end if
 292
 293     vec_dep(1:JPCHAN) = tovs%omb(1:JPCHAN) - SCORR(1:JPCHAN)
 294     vec_abs(1:JPCHAN) = tovs%tb(1:JPCHAN) - SCORR(1:JPCHAN)
 295
 296 ! 3.2 QC:
 297       if (tovs%sensor_id == 3) then
 298            call qc_amsua(tovs)
 299       elseif(tovs%sensor_id == 4) then
 300            call qc_amsub(tovs)
 301       elseif(tovs%sensor_id == 15) then  ! mhs
 302            call qc_amsub(tovs)
 303       elseif(tovs%sensor_id == 10) then  ! ssmis
 304            call qc_ssmis(tovs)
 305       end if
 306
 307 ! 3.3 limb scan check
 308 !---------------------
 309    if (check_limb) then
 310     if ((tovs%scanpos <= 2) .OR. (tovs%scanpos >= (KSCAN-1))) then
 311       nsel(2) = nsel(2) + 1
 312       CYCLE loop2
 313     end if
 314    end if
 315
 316 ! 3.4 Reject data outside radiosonde mask
 317 !------------------------------------------
 318     if (check_mask) then
 319       call MASK(tovs%lat,tovs%lon,LMASK)
 320       if (.NOT. LMASK) then
 321         nsel(8) = nsel(8) + 1
 322         CYCLE loop2
 323       end if
 324     end if
 325
 326 ! 3.5 Reject outliers : facx*sigma, sigma calculated in first pass : loop1
 327 !--------------------------------------------------------------------------
 328     do j=1, tovs%nchan
 329       if ( (abs(vec_dep(j)-vmn(j)) > (vstd(j)*FAC)) ) then
 330         tovs%qc_flag(j) = -1
 331       end if
 332     end do
 333
 334 ! mean/std statistics for relative scan-bias corrected values
 335     do j=1, tovs%nchan
 336       if ( tovs%qc_flag(j) == 1 ) then
 337         jv = j
 338         nobs(j) = nobs(j) + 1
 339         vmean_dep(j) = vmean_dep(j) + vec_dep(j)
 340          vstd_dep(j) = vstd_dep(j) + vec_dep(j)*vec_dep(j)
 341         vmean_abs(j) = vmean_abs(j) + vec_abs(j)
 342          vstd_abs(j) = vstd_abs(j) + vec_abs(j)*vec_abs(j)
 343       end if
 344     end do
 345
 346       PRED(1:tovs%npred) = tovs%pred(1:tovs%npred)
 347
 348 ! compute regression variables mean/var/cov: y:departure; x:predictors
 349     do j=1, tovs%nchan
 350      if ( tovs%qc_flag(j) == 1 ) then
 351         jv = j
 352
 353         ybar(j) = ybar(j) + vec_dep(j)             ! mean of y
 354         ycov(j) = ycov(j) + vec_dep(j)*vec_dep(j)  ! variance of y
 355
 356         do i=1, tovs%npred                         ! Covariances for regression
 357            xbar(j,i) = xbar(j,i) + pred(i)         ! mean of x
 358           xycov(j,i) = xycov(j,i) + vec_dep(j)*pred(i) ! cov of x and y
 359           do ii=1, tovs%npred
 360             xcov(j,i,ii) = xcov(j,i,ii) + pred(i)*pred(ii) ! cov of x and x
 361           end do
 362         end do
 363
 364      end if
 365     end do
 366
 367   end do loop2
 368
 369   365 continue
 370
 371 ! Calculate means, standard deviations and covariances
 372
 373   where (nobs(:) /= 0)
 374     vmean_dep(:) = vmean_dep(:)/nobs(:)
 375     vstd_dep(:)  = vstd_dep(:)/nobs(:) - vmean_dep(:)**2
 376     vstd_dep(:)  = SQRT(MAX(0.0_8,vstd_dep(:)))
 377     vmean_abs(:) = vmean_abs(:)/nobs(:)
 378     vstd_abs(:)  = vstd_abs(:)/nobs(:) - vmean_abs(:)**2
 379     vstd_abs(:)  = SQRT(MAX(0.0_8,vstd_abs(:)))
 380   end where
 381
 382   do j=1, tovs%nchan
 383     if (nobs(j) /= 0) then
 384       ybar(j) = ybar(j)/nobs(j)
 385       ycov(j) = ycov(j)/nobs(j) - ybar(j)*ybar(j)
 386       do i=1, tovs%npred
 387          xbar(j,i) = xbar(j,i)/nobs(j)
 388         xycov(j,i) = xycov(j,i)/nobs(j) - xbar(j,i)*ybar(j)
 389       end do
 390       do i=1, tovs%npred
 391         xcov(j,i,1:tovs%npred) = xcov(j,i,1:tovs%npred)/nobs(j) - xbar(j,i)*xbar(j,1:tovs%npred)
 392       end do
 393     end if
 394   end do
 395
 396   write (6,270) icount,nobs(1:tovs%nchan)
 397
 398   write (6,388)
 399   388 format (/1X,'SECOND PASS: MEANS AND STANDARD DEVIATIONS')
 400
 401   do j=1, tovs%nchan
 402      jv = j
 403     write (6,289) jv, nobs(j), vmean_abs(j), vstd_abs(j), vmean_dep(j), vstd_dep(j)
 404   end do
 405
 406   print *, ' '
 407   print *, 'PREDICTOR MEANS AND STANDARD DEVIATIONS'
 408   do j=1, tovs%nchan
 409      jv = j
 410     print *, ' '
 411     print *, 'CHANNEL ', jv, ' NOBS = ', nobs(j)
 412     do i=1, tovs%npred
 413       write (6,390) i, xbar(j,i), SQRT(xcov(j,i,i))
 414     end do
 415   end do
 416   390 format (1X,I5,4F15.2)
 417
 418 !----------------------------------------------------------------------------
 419 !       4.   CALCULATE REGRESSION COEFFICIENTS.
 420 !            --------- ---------- ------------
 421
 422     do j=1, tovs%nchan
 423        jv = j
 424       if ( nobs(j) >= 10 ) then
 425        print *, 'REGRESSION : CHANNEL ', jv
 426        call REGRESS_ONE(tovs%npred,xbar(j,1:tovs%npred),ybar(j), &
 427                         xcov(j,1:tovs%npred,1:tovs%npred), &
 428                         ycov(j),xycov(j,1:tovs%npred), &
 429                         tovs%npred,coef(j,1:tovs%npred),coef0(j),reserr(j),&
 430                         rescov(j,j),xvec(j,1:tovs%npred,1:tovs%npred))
 431       else
 432        print *, 'nobs < 10: ignoring REGRESSION : CHANNEL ', jv
 433       end if
 434     end do
 435
 436     print *, 'PREDICTOR EIGENVECTORS'
 437     do j=1, tovs%nchan
 438        jv = j
 439       print *, ' '
 440       print *,  'CHANNEL ', jv
 441       do i=1, tovs%npred
 442         print 888, xvec(j,1:tovs%npred,i)
 443       end do
 444     end do
 445     888 format(1X,6F12.4)
 446
 447     xcoef0 = 0.0
 448     xreser = 0.0
 449     xcoef = 0.0
 450
 451     do JJ=1,tovs%nchan
 452       J = JJ
 453       xcoef0(J) = coef0(JJ)
 454       xreser(J) = reserr(JJ)
 455       do I=1, tovs%npred
 456         xcoef(J,I) = coef(JJ,I)
 457       end do
 458     end do
 459
 460 !        5.   OUTPUT RESULTS.
 461 !             ------ -------
 462
 463     coef_year  = real(CDATE/1000000)
 464     coef_month = real(MOD(CDATE,1000000)/10000)
 465     coef_day   = real(MOD(CDATE,10000)/100)                  ! Set date
 466     coef_time  = real(MOD(CDATE,100))
 467
 468     write (6,*)
 469     write (6,502) &
 470                   coef_year, coef_month, coef_day,coef_time
 471     502  format (/1X, '   YEAR',F5.0,'   MONTH',F5.0,'   DAY',F5.0,' TIME',F6.0,//&
 472                   1X,3X,' CH    MEAN  STDDEV  RES.SD  COEFFICIENTS')
 473
 474     do I=1, tovs%nchan
 475       write (6,505) I,vmean_dep(I),vstd_dep(I),xreser(I),(xcoef(I,J),J=1,tovs%npred),xcoef0(I)
 476       505 format (1X,I3,3F8.2,8F12.5)
 477     end do
 478
 479     write (6,508) (I,I=1,tovs%nchan)
 480     508 format (/1X,'RESIDUAL ERROR COVARIANCE'/1X,21I6)
 481
 482     do J=1, tovs%nchan
 483       write (6,510) (rescov(I,J),I=1,tovs%nchan)
 484     end do
 485     510 format (1X,21F6.2)
 486
 487     write (6,511) (I,I=1,tovs%nchan)
 488     511 format (/1X,'RESIDUAL ERROR CORRELATION'/1X,21I6)
 489
 490     do J=1, tovs%nchan
 491       do I=1, tovs%nchan
 492         rescov(I,J) = rescov(I,J)/(reserr(I)*reserr(J))
 493       end do
 494       write (6,510) (rescov(I,J),I=1,tovs%nchan)
 495     end do
 496
 497 !----------------------------------------------------------------------------
 498
 499    deallocate(tovs%tb)
 500    deallocate(tovs%omb)
 501    deallocate(tovs%qc_flag)
 502    deallocate(tovs%cloud_flag)
 503    deallocate(tovs%pred)
 504
 505   close(unit=10)
 506   if (lscan) close(unit=12)
 507
 508 ! out coefs to ASCII file bcor.asc
 509   call write_biascoef(tovs%nchan,kscan,jband,tovs%npred,global, &
 510                       VMNRL(1:tovs%nchan,1:kscan),VMNRLB(1:tovs%nchan,1:kscan,1:jband), &
 511                       xcoef(1:tovs%nchan,1:tovs%npred),xcoef0(1:tovs%nchan), &
 512                       nobs(1:tovs%nchan),vmean_abs(1:tovs%nchan), &
 513                       vstd_abs(1:tovs%nchan),vmean_dep(1:tovs%nchan), &
 514                       vstd_dep(1:tovs%nchan) )
 515
 516   end program da_bias_airmass