var/da/da_gpsref/da_ao_stats_gpsref.inc

   1 subroutine da_ao_stats_gpsref (stats_unit, iv, re)
   2
   3    !--------------------------------------------------------------------
   4    ! Purpose: Called by da_minimisation/da_write_diagnostics.inc.
   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_gpsref_type)         :: stats
  14    integer                          :: ngpsref
  15    integer                          :: n, k
  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_gpsref")
  20
  21    ngpsref = 0
  22    o_minus_b = 0.0; o_minus_a = 0.0; sigma_o = 0.0; sigma_b = 0.0
  23
  24    stats%maximum%ref = maxmin_type (missing_r, 0, 0)
  25    stats%minimum%ref = maxmin_type(-missing_r, 0, 0)
  26
  27    stats%average = residual_gpsref1_type(0.0,0.0,0.0,0.0)
  28    stats%rms_err = stats%average
  29
  30    do n=1, iv%info(gpsref)%nlocal
  31       if (iv%info(gpsref)%proc_domain(1,n)) then
  32          do k=1, iv%info(gpsref)%levels(n)
  33             call da_stats_calculate (n, k, iv%gpsref(n)%ref(k)%qc, &
  34                re%gpsref(n)%ref(k), ngpsref, &
  35                stats%minimum%ref, stats%maximum%ref, &
  36                stats%average%ref, stats%rms_err%ref)
  37
  38             if ( pseudo_ref ) then
  39                o_minus_b = iv%GPSRef(n)%ref(k)%inv
  40                o_minus_a = re%gpsref(n)%ref(k)
  41                sigma_o   = iv%gpsref(n)%ref(k)%error
  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) &
  45                          / o_minus_a) * sigma_o
  46             end if
  47          end do
  48
  49       end if    ! end if (iv%info(gpsref)%proc_domain(1,n))
  50    end do
  51
  52    ! Do inter-processor communication to gather statistics.
  53
  54    call da_proc_sum_int (ngpsref)
  55    iv%nstats(gpsref) = ngpsref
  56
  57    call da_proc_stats_combine(stats%average%ref, stats%rms_err%ref, &
  58       stats%minimum%ref%value, stats%maximum%ref%value, &
  59       stats%minimum%ref%n, stats%maximum%ref%n, &
  60       stats%minimum%ref%l, stats%maximum%ref%l)
  61
  62    if ( pseudo_ref ) then
  63       o_minus_b_0 = wrf_dm_sum_real (o_minus_b)
  64       o_minus_a_0 = wrf_dm_sum_real (o_minus_a)
  65       sigma_o_0 = wrf_dm_sum_real (sigma_o)
  66       sigma_b_0 = wrf_dm_sum_real (sigma_b)
  67       if (rootproc) then
  68          write(stats_unit,'(/A,A3,A,f12.3)')  &
  69             ' Pseudo ', pseudo_var, '     O-B: ', o_minus_b_0
  70          write(stats_unit,' (A,A3,A,f12.3)')  &
  71             ' Pseudo ', pseudo_var, '     O-A: ', o_minus_a_0
  72          write(stats_unit,' (A,A3,A,f12.3)')  &
  73             ' Pseudo ', pseudo_var, ' sigma_o: ', sigma_o_0
  74          write(stats_unit,'(A,A3,A,f12.3)')  &
  75             ' Pseudo ', pseudo_var, ' sigma_b: ', sigma_b_0
  76       end if
  77    end if
  78
  79    if (rootproc) then
  80       if (ngpsref > 0) then
  81          write(unit=stats_unit, fmt='(/a/)') ' Diagnostics of AO for gpsref'
  82             call da_print_stats_gpsref(stats_unit, ngpsref, stats)
  83       end if
  84    end if
  85
  86    if (trace_use_dull) call da_trace_exit("da_ao_stats_gpsref")
  87
  88 end subroutine da_ao_stats_gpsref
  89
  90