var/da/da_gpspw/da_ao_stats_gpspw.inc

   1 subroutine da_ao_stats_gpspw (stats_unit, iv, re)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: TBD
   5    !-----------------------------------------------------------------------
   6
   7    implicit none
   8
   9    integer,        intent(in)    :: stats_unit    ! Output unit for stats.
  10    type (iv_type), intent(inout) :: iv            ! iv
  11    type (y_type),  intent(in)    :: re            ! A - O
  12
  13    type (stats_gpspw_type) :: stats
  14    integer                 :: ntpw
  15    integer                 :: n
  16    real                    :: o_minus_b, o_minus_a, sigma_o, sigma_b
  17    real                    :: o_minus_b_0, o_minus_a_0, sigma_o_0, sigma_b_0
  18
  19    if (trace_use_dull) call da_trace_entry("da_ao_stats_gpspw")
  20
  21    ntpw = 0
  22    o_minus_b = 0.0; o_minus_a = 0.0; sigma_o = 0.0; sigma_b = 0.0
  23
  24    stats%maximum%tpw = maxmin_type (missing_r, 0, 0)
  25    stats%minimum%tpw = maxmin_type(-missing_r, 0, 0)
  26
  27    stats%average = residual_gpspw1_type(0.0)
  28    stats%rms_err = stats%average
  29
  30    do n=1, iv%info(gpspw)%nlocal
  31       if (iv%info(gpspw)%proc_domain(1,n)) then
  32          call da_stats_calculate (n, 0, iv%gpspw(n)%tpw%qc, &
  33             re%gpspw(n)%tpw, ntpw,   &
  34             stats%minimum%tpw,  stats%maximum%tpw, &
  35             stats%average%tpw,  stats%rms_err%tpw)
  36
  37          if ( pseudo_tpw .or. pseudo_ztd ) then
  38             o_minus_b = iv%GPSpw(n)%tpw%inv
  39             o_minus_a = re%gpspw(n)%tpw
  40             sigma_o   = iv%gpspw(n)%tpw%error
  41
  42             ! Calculate equivalent sigma_b using
  43             ! O-A=(O-B)*sigma_o/(sigma_o+sigma_b)
  44             sigma_b = sqrt ((o_minus_b - o_minus_a) / o_minus_a) * sigma_o
  45          end if
  46       end if    ! end if (iv%info(gpspw)%proc_domain(1,n))
  47    end do
  48
  49    ! Do inter-processor communication to gather statistics.
  50    call da_proc_sum_int (ntpw)
  51    iv%nstats(gpspw) = ntpw
  52
  53    call da_proc_stats_combine(stats%average%tpw, stats%rms_err%tpw, &
  54       stats%minimum%tpw%value, stats%maximum%tpw%value, &
  55       stats%minimum%tpw%n, stats%maximum%tpw%n, &
  56       stats%minimum%tpw%l, stats%maximum%tpw%l)
  57
  58    if ( pseudo_tpw .or. pseudo_ztd ) then
  59       o_minus_b_0 = wrf_dm_sum_real (o_minus_b)
  60       o_minus_a_0 = wrf_dm_sum_real (o_minus_a)
  61       sigma_o_0 = wrf_dm_sum_real (sigma_o)
  62       sigma_b_0 = wrf_dm_sum_real (sigma_b)
  63       if (rootproc) then
  64          write(stats_unit,'(/A,A3,A,f12.3)')  &
  65             ' Pseudo ', pseudo_var, '     O-B: ', o_minus_b_0
  66          write(stats_unit,' (A,A3,A,f12.3)')  &
  67             ' Pseudo ', pseudo_var, '     O-A: ', o_minus_a_0
  68          write(stats_unit,' (A,A3,A,f12.3)')  &
  69             ' Pseudo ', pseudo_var, ' sigma_o: ', sigma_o_0
  70          write(stats_unit,'(A,A3,A,f12.3)')  &
  71             ' Pseudo ', pseudo_var, ' sigma_b: ', sigma_b_0
  72       end if
  73    end if
  74
  75    if (rootproc) then
  76       if (ntpw /= 0) then
  77          if (use_gpspwObs) then
  78            write(unit=stats_unit, fmt='(/a/)') ' Diagnostics of AO for gpspw'
  79          else if (use_gpsztdObs) then
  80            write(unit=stats_unit, fmt='(/a/)') ' Diagnostics of AO for gpsztd'
  81          endif
  82
  83          call da_print_stats_gpspw(stats_unit, ntpw, stats)
  84        end if
  85    end if
  86
  87    if (trace_use_dull) call da_trace_exit("da_ao_stats_gpspw")
  88
  89 end subroutine da_ao_stats_gpspw
  90
  91