var/da/da_varbc/da_varbc_init.inc

   1   subroutine da_varbc_init (iv, be)
   2
   3    !---------------------------------------------------------------------------
   4    !  PURPOSE: Initialize Variational bias correction from VARBC.in file
   5    !
   6    !  Called from da_radiance_init routine
   7    !
   8    !  HISTORY: 10/26/2007 - Creation                     Tom Auligne
   9    !---------------------------------------------------------------------------
  10
  11    implicit none
  12
  13    type (iv_type), intent(inout)  :: iv       ! O-B structure.
  14
  15    type (be_type), intent(inout) :: be        ! background error.
  16
  17    !
  18    !  local arguments
  19    !-------------------
  20    integer   :: n, i, ii, j, k, nchan, npred, npred_ws, ipred, ichan, np, nsensor
  21    integer   :: iunit, iost
  22    integer   :: size_jp, jp_start
  23    integer   :: platform_id, satellite_id, sensor_id, nchanl, npredmax, num_rad
  24    integer, allocatable :: pred_use(:), ipred_ws(:)
  25    real, allocatable    :: pred_ws(:)
  26    character(len=filename_len) :: filename
  27    character(len=120) :: cline
  28    logical   :: lvarbc_read, limatch, lmatch
  29
  30    if (trace_use) call da_trace_entry("da_varbc_init")
  31
  32    !--------------------------------------------------------------
  33    ! [1] Initializations
  34    !--------------------------------------------------------------
  35    size_jp  = 0                            !! Jp Control Variable size
  36    jp_start = be % cv % size_jb + be % cv % size_je
  37
  38    do n = 1, iv % num_inst
  39       nchan = iv%instid(n)%nchan
  40       iv%instid(n)%varbc_info%npredmax = 0
  41       allocate ( iv%instid(n)%varbc(nchan) )
  42       do k = 1, nchan
  43          iv%instid(n)%varbc(k)%npred = 0   ! by default, do not use VarBC
  44          iv%instid(n)%varbc(k)%nobs  = 0
  45       end do
  46    end do
  47
  48    !--------------------------------------------------------------
  49    ! [2] Read VarBC.in file and match with obs data set
  50    !--------------------------------------------------------------
  51    filename = 'VARBC.in'
  52    inquire(file=trim(adjustl(filename)), exist=lvarbc_read)
  53
  54    if (lvarbc_read) then
  55     write(unit=stdout,fmt='(A)') 'VARBC: Reading VARBC.in file'
  56     call da_get_unit(iunit)
  57     open(unit=iunit,file=filename, form='formatted', status='old')
  58
  59     read(iunit, *)
  60     read(iunit, *) nsensor
  61
  62     do j = 1, nsensor
  63       read(iunit, *) ; read(iunit, *) ;read(iunit, *)
  64       read(iunit, *) platform_id, satellite_id, sensor_id, nchanl, npredmax
  65       limatch = .false.
  66        do n = 1, iv % num_inst
  67          if ( (platform_id  == iv%instid(n)%platform_id) .and. &   !! found match !!
  68               (satellite_id == iv%instid(n)%satellite_id) .and. &  !! found match !!
  69               (sensor_id    == iv%instid(n)%sensor_id) ) then      !! found match !!
  70             limatch = .true.                                       !! found match !!
  71
  72 !            write(unit=stdout,fmt='(A,A)') &
  73 !              'VARBC: Matching with obs for ', iv%instid(n)%rttovid_string
  74
  75             num_rad  = iv%instid(n)%num_rad
  76             allocate ( iv%instid(n)%varbc_info%pred(npredmax, num_rad) )
  77
  78             iv%instid(n)%varbc_info%platform_id  = platform_id
  79             iv%instid(n)%varbc_info%satellite_id = satellite_id
  80             iv%instid(n)%varbc_info%sensor_id    = sensor_id
  81             iv%instid(n)%varbc_info%nchanl       = nchanl
  82             iv%instid(n)%varbc_info%npredmax     = npredmax
  83             iv%instid(n)%varbc_info%gammapred    = 9    ! CRTM Gamma Correction
  84             read(iunit, *)
  85             allocate ( iv%instid(n)%varbc_info%pred_mean(npredmax) )
  86             allocate ( iv%instid(n)%varbc_info%pred_std (npredmax) )
  87             allocate ( iv%instid(n)%varbc_info%nbgerr   (npredmax) )
  88             read(iunit, *) iv%instid(n)%varbc_info%pred_mean
  89             read(iunit, *) iv%instid(n)%varbc_info%pred_std
  90             read(iunit, *) iv%instid(n)%varbc_info%nbgerr
  91             read(iunit, *)
  92
  93             allocate ( pred_use(npredmax), ipred_ws(npredmax), pred_ws(npredmax) )
  94             do i = 1, nchanl
  95                lmatch = .false.
  96                read(iunit, '(A)',iostat=iost) cline
  97                read(cline, *) ii, ichan, pred_use
  98                npred    = COUNT (pred_use >= 0)
  99                if (npred <= 0) cycle                !! VarBC channels only
 100
 101                chan_loop: do k = 1, iv%instid(n)%nchan
 102                   if (iv%instid(n)%ichan(k) == ichan) then         !! found match !!
 103                      lmatch = .true.                               !! found match !!
 104                      iv%instid(n)%varbc(k)%ichanl = ichan
 105                      iv%instid(n)%varbc(k)%npred  = npred
 106                      allocate ( iv%instid(n)%varbc(k)%pred_use (npredmax) )
 107                      allocate ( iv%instid(n)%varbc(k)%ipred (npred) )
 108                      allocate ( iv%instid(n)%varbc(k)%index (npred) )
 109                      allocate ( iv%instid(n)%varbc(k)%param (npred) )
 110                      allocate ( iv%instid(n)%varbc(k)%bgerr (npred) )
 111                      allocate ( iv%instid(n)%varbc(k)%vtox    (npred,npred) )
 112                      iv%instid(n)%varbc(k)%pred_use = pred_use
 113                      iv%instid(n)%varbc(k)%vtox(1:npred, 1:npred)     = 0.0
 114                      do ipred = 1, npred
 115                         iv%instid(n)%varbc(k)%vtox(ipred, ipred)      = 1.0
 116                      end do
 117                      iv%instid(n)%varbc(k)%param(1:npred)             = 0.0
 118                      iv%instid(n)%varbc(k)%bgerr(1:npred)             = 0.0
 119                      iv%instid(n)%varbc(k)%ipred(1:npred)             = &
 120                            PACK((/ (ipred, ipred=1,npredmax) /), mask = (pred_use >= 0))
 121
 122                     ! VarBC warm-start parameters
 123                     !-----------------------------
 124                      npred_ws = COUNT (pred_use > 0)
 125                      if (npred_ws > 0) then
 126                         read(cline, *) ii, ichan, pred_use, pred_ws(1:npred_ws)
 127                         ipred_ws(1:npred_ws) = PACK((/ (ipred, ipred=1,npred) /), &
 128                                                mask = (pred_use(iv%instid(n)%varbc(k)%ipred) >  0))
 129                         iv%instid(n)%varbc(k)%param(ipred_ws(1:npred_ws))= pred_ws(1:npred_ws)
 130                      end if
 131
 132                     ! Jp Control Variable size
 133                     !--------------------------
 134                      do ipred = 1, npred
 135                          size_jp =  size_jp + 1
 136                          iv%instid(n)%varbc(k)%index(ipred) = jp_start + size_jp
 137                      end do
 138
 139                      exit chan_loop
 140                   end if
 141                end do chan_loop
 142
 143                if (.not. lmatch) write(unit=stdout,fmt='(A,I4)') &
 144                                        'VARBC: no matching for channel:',ichan
 145             end do
 146             deallocate(pred_use, ipred_ws, pred_ws)
 147          end if
 148        end do
 149        if (.not. limatch) then
 150           read(iunit, *) ; read(iunit, *) ; read(iunit, *) ; read(iunit, *) ; read(iunit, *)
 151           do i = 1, nchanl
 152              read(iunit, *)
 153           end do
 154           write(unit=stdout,fmt='(A,3I4)') &
 155                 'VARBC: no matching for platform/satid/sensor',platform_id, satellite_id, sensor_id
 156        end if
 157     end do
 158     close(iunit)
 159     call da_free_unit(iunit)
 160    else
 161       write(unit=stdout,fmt='(A)') 'VARBC: could not find VARBC.in file --> VARBC switched off'
 162       use_varbc    = .false.
 163       freeze_varbc = .false.
 164    end if
 165
 166    !--------------------------------------------------------------
 167    ! [3] Define VarBC control variable size:
 168    !--------------------------------------------------------------
 169    use_varbc = use_varbc.and.(.not.freeze_varbc)
 170    if (use_varbc) then
 171       be % cv % size_jp = size_jp
 172       cv_size_domain_jp = size_jp
 173    end if
 174
 175    if (trace_use) call da_trace_exit("da_varbc_init")
 176
 177  end subroutine da_varbc_init