var/da/da_test/da_check_vp_errors.inc

   1 subroutine da_check_vp_errors(vp1, vp2, ne, &
   2                                its,ite, jts,jte, kts,kte)
   3
   4    !---------------------------------------------------------------------------
   5    ! Purpose: Test invertibility of transform to/from Vp or Vv
   6    !
   7    ! Method:  Perform statistics on differences in initial and final Vv or Vp
   8    !---------------------------------------------------------------------------
   9
  10    implicit none
  11
  12    type (vp_type), intent(in)     :: vp1         ! Test input
  13    type (vp_type), intent(in)     :: vp2         ! Test output.
  14    integer, intent(in)            :: ne          ! Ensemble size.
  15    integer, intent(in)            :: its,ite, jts,jte, kts,kte ! tile   dims.
  16
  17    real                           :: inv_size    ! 1/size of array.
  18    real                           :: rms_fild    ! RMS of field.
  19    real                           :: rms_diff    ! RMS of differnce.
  20
  21    real, dimension(its:ite, jts:jte, kts:kte) :: diff ! Difference
  22 !  real                           :: diff_alpha(its:ite,jts:jte,kts:kte,1:ne)
  23    real                           :: diff_alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:ne)
  24    real                           :: inv_size_ens    ! 1/size of ensemble array.
  25
  26    if (trace_use) call da_trace_entry("da_check_vp_errors")
  27
  28    inv_size = 1.0 / real((ite-its+1) * (jte-jts+1) * (kte-kts+1))
  29    inv_size_ens = 1.0 / real((ite_int-its_int+1) * (jte_int-jts_int+1) * (kte_int-kts_int+1)*ne)
  30
  31    !-------------------------------------------------------------------------
  32    ! [1.0]: Check v1 differences:
  33    !-------------------------------------------------------------------------
  34
  35    diff(its:ite,jts:jte, kts:kte) = vp2 % v1(its:ite,jts:jte,kts:kte) - &
  36                                     vp1 % v1(its:ite,jts:jte,kts:kte)
  37
  38    rms_fild = sqrt(sum(vp1 % v1(its:ite, jts:jte,kts:kte) &
  39                        * vp1 % v1(its:ite, jts:jte,kts:kte)) * inv_size)
  40    rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
  41                        * diff(its:ite, jts:jte,kts:kte)) * inv_size)
  42
  43    if (rms_fild == 0.0) then
  44       write(unit=stdout, fmt='(a)') ' v1 is zero '
  45    else
  46       write(unit=stdout, fmt='(a,1pe10.4)') &
  47            ' v1 RMS error/RMS field = ', rms_diff/rms_fild
  48    end if
  49
  50    !-------------------------------------------------------------------------
  51    ! [2.0]: Check v2 differences:
  52    !-------------------------------------------------------------------------
  53
  54    diff(its:ite,jts:jte, kts:kte) = vp2 % v2(its:ite,jts:jte,kts:kte) - &
  55                                     vp1 % v2(its:ite,jts:jte,kts:kte)
  56
  57    rms_fild = sqrt(sum(vp1 % v2(its:ite, jts:jte,kts:kte) &
  58                        * vp1 % v2(its:ite, jts:jte,kts:kte)) * inv_size)
  59    rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
  60                        * diff(its:ite, jts:jte,kts:kte)) * inv_size)
  61
  62    if (rms_fild == 0.0) then
  63       write(unit=stdout, fmt='(a)') ' v2 is zero '
  64    else
  65       write(unit=stdout, fmt='(a,1pe10.4)') &
  66            ' v2 RMS error/RMS field = ', rms_diff/rms_fild
  67    end if
  68
  69    !-------------------------------------------------------------------------
  70    ! [3.0]: Check v3 differences:
  71    !-------------------------------------------------------------------------
  72
  73    diff(its:ite,jts:jte, kts:kte) = vp2 % v3(its:ite,jts:jte,kts:kte) - &
  74                                     vp1 % v3(its:ite,jts:jte,kts:kte)
  75
  76    rms_fild = sqrt(sum(vp1 % v3(its:ite, jts:jte,kts:kte) &
  77                        * vp1 % v3(its:ite, jts:jte,kts:kte)) * inv_size)
  78    rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
  79                        * diff(its:ite, jts:jte,kts:kte)) * inv_size)
  80
  81    if (rms_fild == 0.0) then
  82       write(unit=stdout, fmt='(a)') ' v3 is zero '
  83    else
  84       write(unit=stdout, fmt='(a,1pe10.4)') &
  85          ' v3 RMS error/RMS field = ', rms_diff/rms_fild
  86    end if
  87
  88    !-------------------------------------------------------------------------
  89    ! [4.0]: Check v4 differences:
  90    !-------------------------------------------------------------------------
  91
  92    diff(its:ite,jts:jte, kts:kte) = vp2 % v4(its:ite,jts:jte,kts:kte) - &
  93                                     vp1 % v4(its:ite,jts:jte,kts:kte)
  94
  95    rms_fild = sqrt(sum(vp1 % v4(its:ite, jts:jte,kts:kte) &
  96                        * vp1 % v4(its:ite, jts:jte,kts:kte)) * inv_size)
  97    rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
  98                        * diff(its:ite, jts:jte,kts:kte)) * inv_size)
  99
 100    if (rms_fild == 0.0) then
 101       write(unit=stdout, fmt='(a)') ' v4 is zero '
 102    else
 103       write(unit=stdout, fmt='(a,1pe10.4)') &
 104          ' v4 RMS error/RMS field = ', rms_diff/rms_fild
 105    end if
 106
 107    !-------------------------------------------------------------------------
 108    ! [5.0]: Check v5 differences:
 109    !-------------------------------------------------------------------------
 110
 111    inv_size = 1.0 / real((ite-its+1) * (jte-jts+1))
 112
 113    diff(its:ite, jts:jte,kts:kte) = vp2 % v5(its:ite, jts:jte,kts:kte) - &
 114       vp1 % v5(its:ite, jts:jte,kts:kte)
 115
 116    rms_fild = sqrt(sum(vp1 % v5(its:ite, jts:jte,kts:kte) * &
 117       vp1 % v5(its:ite, jts:jte,kts:kte)) * inv_size)
 118    rms_diff = sqrt(sum(diff(its:ite, jts:jte,1) * &
 119       diff(its:ite, jts:jte,1)) * inv_size)
 120
 121    if (rms_fild == 0.0) then
 122       write(unit=stdout, fmt='(a)') ' v5 is zero '
 123    else
 124       write(unit=stdout, fmt='(a,1pe10.4)') &
 125          ' v5 RMS error/RMS field = ', rms_diff/rms_fild
 126    end if
 127
 128    !-------------------------------------------------------------------------
 129    ! [6.0]: Check alpha differences:
 130    !-------------------------------------------------------------------------
 131
 132    ! changed all dimensions below to reflect ensemble tiles
 133    if (ne > 0) then
 134 !     inv_size = 1.0 / real((ite-its+1) * (jte-jts+1) * ne)
 135       diff_alpha(its:ite,jts:jte,kts:kte,1:ne) = vp2 % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:ne) - &
 136                                          vp1 % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:ne)
 137       rms_fild = sqrt(sum(vp1 % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:ne) &
 138                           * vp1 % alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:ne)) * inv_size_ens)
 139       rms_diff = sqrt(sum(diff_alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:ne) &
 140                           * diff_alpha(its_int:ite_int,jts_int:jte_int,kts_int:kte_int,1:ne)) * inv_size_ens)
 141
 142       if (rms_fild /= 0.0) then
 143          write(unit=stdout, fmt='(a,1pe10.4)') ' alpha RMS error/RMS field = ',&
 144             rms_diff/rms_fild
 145       end if
 146    end if
 147
 148    if (trace_use) call da_trace_exit("da_check_vp_errors")
 149
 150 end subroutine da_check_vp_errors
 151
 152