var/da/da_varbc/da_varbc_update.inc

   1   subroutine da_varbc_update (it, cv_size, cv, iv)
   2
   3    !---------------------------------------------------------------------------
   4    !  PURPOSE: Update VarBC parameters and write into file
   5    !
   6    ! [1] Update of VarBC parameters (including change of variable)
   7    !
   8    ! [2] Write VARBC.out file with VarBC information
   9    !
  10    !  Called from da_solve
  11    !
  12    !  HISTORY: 10/29/2007 - Creation                          Tom Auligne
  13    !           01/20/2010 - Update for multiple outer-loops   Tom Auligne
  14    !---------------------------------------------------------------------------
  15
  16    implicit none
  17
  18    integer,        intent(in)    :: it          !outer loop counting
  19    integer,        intent(in)    :: cv_size
  20    real,           intent(in)    :: cv(cv_size) ! Control variable structure.
  21    type (iv_type), intent(inout) :: iv          ! Obs. increment structure.
  22
  23    !
  24    !  local arguments
  25    !-------------------
  26    integer              :: inst, ichan, npred, i, npredmax, id
  27    integer              :: iunit, iost
  28    character(len=filename_len) :: filename
  29    character(len=120)   :: fparam
  30    real, allocatable    :: varbc_param_tl(:)
  31    logical, allocatable :: lvarbc_inst(:)
  32
  33    if (trace_use) call da_trace_entry("da_varbc_update")
  34
  35    write(unit=stdout,fmt='(A)') 'VARBC: Updating bias parameters'
  36
  37    allocate( lvarbc_inst(iv%num_inst) )
  38    do inst = 1, iv % num_inst
  39       lvarbc_inst(inst) = ANY(iv%instid(inst)%varbc(1:iv%instid(inst)%nchan)%npred>0)
  40    end do
  41
  42    do inst = 1, iv % num_inst
  43       if (.not. lvarbc_inst(inst)) cycle             !! VarBC instruments only
  44       allocate(varbc_param_tl(iv%instid(inst)%varbc_info%npredmax))
  45
  46       do ichan = 1, iv%instid(inst)%nchan
  47          npred    = iv%instid(inst)%varbc(ichan)%npred
  48          if (npred <= 0) cycle               !! VarBC channels only
  49          if (iv%instid(inst)%varbc(ichan)%nobs >= varbc_nobsmin) then
  50             where (iv%instid(inst)%varbc(ichan)%pred_use == 0) &
  51                    iv%instid(inst)%varbc(ichan)%pred_use = 1
  52          else
  53             if ( satinfo(inst)%iuse(ichan) == 1) then
  54                if (count(iv%instid(inst)%varbc(ichan)%pred_use == 0) > 0) &
  55                   write(unit=stdout,fmt='(A,A,I5)') &
  56                   'VARBC: Not enough data to keep statistics for ', &
  57                   trim(iv%instid(inst)%rttovid_string),             &
  58                   iv%instid(inst)%ichan(ichan)
  59             end if
  60          end if
  61
  62          if (.not. use_varbc) cycle
  63
  64         !---------------------------------------------------------------
  65         ! Change of variable (preconditioning) for parameters increments
  66         !---------------------------------------------------------------
  67
  68          do i = 1, npred
  69             id = iv%instid(inst)%varbc(ichan)%index(i)
  70             varbc_param_tl(i) = &
  71                SUM(cv(id) * iv%instid(inst)%varbc(ichan)%vtox(i,1:npred))
  72          end do
  73
  74         !---------------------------------------------------------------
  75         ! Update VarBC parameters
  76         !---------------------------------------------------------------
  77          iv%instid(inst)%varbc(ichan)%param = &
  78          iv%instid(inst)%varbc(ichan)%param + varbc_param_tl(1:npred)
  79       end do
  80
  81       deallocate(varbc_param_tl)
  82    end do
  83
  84    if (.not. rootproc) then
  85       if (trace_use) call da_trace_exit("da_varbc_update")
  86       return
  87    end if
  88
  89    !---------------------------------------------------------------
  90    ! Write VARBC.out file
  91    !---------------------------------------------------------------
  92
  93    write(unit=stdout,fmt='(A)') 'VARBC: Writing information in VARBC.out file'
  94
  95
  96    call da_get_unit(iunit)
  97    if ( it == max_ext_its ) then
  98       filename = 'VARBC.out'
  99    else
 100       write(unit=filename, fmt='(a,i2.2)') 'VARBC.out_',it
 101    end if
 102    open(unit=iunit,file=filename,form='formatted',iostat = iost,status='replace')
 103
 104    if (iost /= 0) then
 105       message(1)="Cannot open bias correction file "//adjustl(filename)
 106       call da_error(__FILE__,__LINE__,message(1:1))
 107    end if
 108
 109    write(iunit, *) ' VARBC version 1.0 - Number of instruments: '
 110    write(iunit, *) COUNT(lvarbc_inst)
 111
 112    do inst = 1, iv % num_inst
 113       if (.not. lvarbc_inst(inst)) cycle             !! VarBC instruments only
 114       npredmax = iv%instid(inst)%varbc_info%npredmax
 115
 116       write(iunit, *) '------------------------------------------------'
 117       write(iunit, *) 'Platform_id  Sat_id  Sensor_id  Nchanl  Npredmax'
 118       write(iunit, *) '------------------------------------------------'
 119       write(iunit, *) iv%instid(inst)%varbc_info%platform_id, &
 120                       iv%instid(inst)%varbc_info%satellite_id, &
 121                       iv%instid(inst)%varbc_info%sensor_id, &
 122                       iv%instid(inst)%varbc_info%nchanl,&
 123                       iv%instid(inst)%varbc_info%npredmax
 124
 125       write(iunit, *) ' -----> Bias predictor statistics:  Mean & Std & Nbgerr'
 126       write(iunit,'(30F10.1)') iv%instid(inst)%varbc_info%pred_mean(1:npredmax)
 127       write(iunit,'(30F10.1)') iv%instid(inst)%varbc_info%pred_std (1:npredmax)
 128       write(iunit,'(30I10)')   iv%instid(inst)%varbc_info%nbgerr   (1:npredmax)
 129
 130       write(iunit, *) ' -----> Chanl_id Chanl_nb  Pred_use(-1/0/1)  Param'
 131
 132       do ichan = 1, iv%instid(inst)%nchan
 133          npred    = iv%instid(inst)%varbc(ichan)%npred
 134          if (npred <= 0) cycle               !! VarBC channels only
 135          write(fparam,*) '(I4,I6,',npredmax,'I3,',npred,'F8.3)'
 136          write(iunit, fmt=trim(adjustl(fparam))) &
 137              ichan, iv%instid(inst)%varbc(ichan)%ichanl,   &
 138                     iv%instid(inst)%varbc(ichan)%pred_use, &
 139                     iv%instid(inst)%varbc(ichan)%param
 140       end do
 141    end do
 142
 143    deallocate(lvarbc_inst)
 144    close(iunit)
 145    call da_free_unit(iunit)
 146
 147    if (trace_use) call da_trace_exit("da_varbc_update")
 148
 149
 150  end subroutine da_varbc_update