var/da/da_test/da_check_xtovptox_errors.inc

   1 subroutine da_check_xtovptox_errors(xa, xa2_u, xa2_v, xa2_w, xa2_t, &
   2                                      xa2_p, xa2_q, xa2_rho, &
   3                                      xa2_qt, xa2_qcw, xa2_qrn)
   4
   5    !---------------------------------------------------------------------------
   6    ! Purpose: Test invertibility of v = U^{-1} x followed by x = Uv.
   7    !
   8    !  Method:  Perform statistics on differences in initial and final x.
   9    !---------------------------------------------------------------------------
  10
  11    implicit none
  12
  13    type (x_type), intent(in)      :: xa          ! Test input
  14
  15    real, dimension(ims:ime, jms:jme, kms:kme), &
  16                  intent(in)      :: xa2_u, xa2_v, xa2_t, &
  17                                     xa2_p, xa2_q, xa2_rho, &
  18                                     xa2_qt, xa2_qcw, xa2_qrn
  19    real, dimension(ims:ime, jms:jme, kms:kme), &
  20                  intent(in)      :: xa2_w    !xiao
  21
  22
  23    real                           :: rms_fild    ! RMS of field.
  24    real                           :: rms_diff    ! RMS of differnce.
  25
  26    real, dimension(ims:ime, jms:jme, kms:kme) :: diff ! Difference
  27
  28    if (trace_use) call da_trace_entry("da_check_xtovpx_errors")
  29
  30    !----------------------------------------------------------------------
  31    ! [1.0]: Check u differences:
  32    !----------------------------------------------------------------------
  33
  34    diff(its:ite, jts:jte, kts:kte) = xa2_u(its:ite, jts:jte, kts:kte) &
  35                                    - xa% u(its:ite, jts:jte, kts:kte)
  36
  37    rms_fild = sqrt(sum(xa % u(its:ite, jts:jte, kts:kte) &
  38                        * xa % u(its:ite, jts:jte, kts:kte)))
  39    rms_diff = sqrt(sum(diff(its:ite, jts:jte, kts:kte) &
  40                        * diff(its:ite, jts:jte, kts:kte)))
  41
  42    if (rms_fild == 0.0) then
  43       write(unit=stdout, fmt='(a)') ' u is zero '
  44    else
  45       write(unit=stdout, fmt='(a,1pe10.4)') ' u RMS error = ', rms_diff
  46       write(unit=stdout, fmt='(a,1pe10.4)') ' u RMS field = ', rms_fild
  47       write(unit=stdout, fmt='(a,1pe10.4)') ' u RMS error/RMS field = ', &
  48          rms_diff/rms_fild
  49    end if
  50
  51    !----------------------------------------------------------------------
  52    ! [2.0]: Check v differences:
  53    !----------------------------------------------------------------------
  54
  55    diff(its:ite, jts:jte, kts:kte) = xa2_v(its:ite, jts:jte, kts:kte) &
  56                                    - xa% v(its:ite, jts:jte, kts:kte)
  57
  58    rms_fild = sqrt(sum(xa % v(its:ite, jts:jte, kts:kte) &
  59                        * xa % v(its:ite, jts:jte, kts:kte)))
  60    rms_diff = sqrt(sum(diff(its:ite, jts:jte, kts:kte) &
  61                        * diff(its:ite, jts:jte, kts:kte)))
  62
  63    if (rms_fild == 0.0) then
  64       write(unit=stdout, fmt='(a)') ' v is zero '
  65    else
  66       write(unit=stdout, fmt='(a,1pe10.4)') ' v RMS error = ', rms_diff
  67       write(unit=stdout, fmt='(a,1pe10.4)') ' v RMS field = ', rms_fild
  68       write(unit=stdout, fmt='(a,1pe10.4)') ' v RMS error/RMS field = ', &
  69       rms_diff/rms_fild
  70    end if
  71
  72    !----------------------------------------------------------------------
  73    ! [3.0]: Check t differences:
  74    !----------------------------------------------------------------------
  75
  76    diff(its:ite, jts:jte, kts:kte) = xa2_t(its:ite, jts:jte, kts:kte) &
  77                                    - xa% t(its:ite, jts:jte, kts:kte)
  78
  79    rms_fild = sqrt(sum(xa % t(its:ite, jts:jte, kts:kte) &
  80                        * xa % t(its:ite, jts:jte, kts:kte)))
  81    rms_diff = sqrt(sum(diff(its:ite, jts:jte, kts:kte) &
  82                        * diff(its:ite, jts:jte, kts:kte)))
  83
  84    if (rms_fild == 0.0) then
  85       write(unit=stdout, fmt='(a)') ' t is zero '
  86    else
  87       write(unit=stdout, fmt='(a,1pe10.4)') ' t RMS error = ', rms_diff
  88       write(unit=stdout, fmt='(a,1pe10.4)') ' t RMS field = ', rms_fild
  89       write(unit=stdout, fmt='(a,1pe10.4)') ' t RMS error/RMS field = ', &
  90          rms_diff/rms_fild
  91    end if
  92
  93    !----------------------------------------------------------------------
  94    ! [4.0]: Check p differences:
  95    !----------------------------------------------------------------------
  96
  97    diff(its:ite, jts:jte, kts:kte) = xa2_p(its:ite, jts:jte, kts:kte) &
  98                                    - xa% p(its:ite, jts:jte, kts:kte)
  99
 100    rms_fild = sqrt(sum(xa % p(its:ite, jts:jte, kts:kte) &
 101                        * xa % p(its:ite, jts:jte, kts:kte)))
 102    rms_diff = sqrt(sum(diff(its:ite, jts:jte, kts:kte) &
 103                        * diff(its:ite, jts:jte, kts:kte)))
 104
 105    if (rms_fild == 0.0) then
 106       write(unit=stdout, fmt='(a)') ' p is zero '
 107    else
 108       write(unit=stdout, fmt='(a,1pe10.4)') ' p RMS error = ', rms_diff
 109       write(unit=stdout, fmt='(a,1pe10.4)') ' p RMS field = ', rms_fild
 110       write(unit=stdout, fmt='(a,1pe10.4)') ' p RMS error/RMS field = ', &
 111          rms_diff/rms_fild
 112    end if
 113
 114    !----------------------------------------------------------------------
 115    ! [5.0]: Check q differences:
 116    !----------------------------------------------------------------------
 117
 118    diff(its:ite, jts:jte, kts:kte) = xa2_q(its:ite, jts:jte, kts:kte) &
 119                                    - xa% q(its:ite, jts:jte, kts:kte)
 120
 121    rms_fild = sqrt(sum(xa % q(its:ite, jts:jte, kts:kte) &
 122                        * xa % q(its:ite, jts:jte, kts:kte)))
 123    rms_diff = sqrt(sum(diff(its:ite, jts:jte, kts:kte) &
 124                        * diff(its:ite, jts:jte, kts:kte)))
 125
 126    if (rms_fild == 0.0) then
 127       write(unit=stdout, fmt='(a)') ' q is zero '
 128    else
 129       write(unit=stdout, fmt='(a,1pe10.4)') ' q RMS error = ', rms_diff
 130       write(unit=stdout, fmt='(a,1pe10.4)') ' q RMS field = ', rms_fild
 131       write(unit=stdout, fmt='(a,1pe10.4)') ' q RMS error/RMS field = ', &
 132          rms_diff/rms_fild
 133    end if
 134
 135    !----------------------------------------------------------------------
 136    ! [6.0]: Check rho differences:
 137    !----------------------------------------------------------------------
 138
 139    diff(its:ite, jts:jte, kts:kte) = xa2_rho(its:ite, jts:jte, kts:kte) &
 140                                    - xa% rho(its:ite, jts:jte, kts:kte)
 141
 142    rms_fild = sqrt(sum(xa % rho(its:ite, jts:jte, kts:kte) &
 143                        * xa % rho(its:ite, jts:jte, kts:kte)))
 144    rms_diff = sqrt(sum(diff(its:ite, jts:jte, kts:kte) &
 145                        * diff(its:ite, jts:jte, kts:kte)))
 146
 147    if (rms_fild == 0.0) then
 148       write(unit=stdout, fmt='(a)') ' rho is zero '
 149    else
 150       write(unit=stdout, fmt='(a,1pe10.4)') ' rho RMS error = ', rms_diff
 151       write(unit=stdout, fmt='(a,1pe10.4)') ' rho RMS field = ', rms_fild
 152       write(unit=stdout, fmt='(a,1pe10.4)') ' rho RMS error/RMS field = ', &
 153          rms_diff/rms_fild
 154    end if
 155
 156    !----------------------------------------------------------------------
 157    ! [7.0]: Check w differences:
 158    !----------------------------------------------------------------------
 159
 160    diff(its:ite, jts:jte, kts:kte+1) = xa2_w(its:ite, jts:jte, kts:kte+1) &
 161                                      - xa% w(its:ite, jts:jte, kts:kte+1)
 162
 163    rms_fild = sqrt(sum(xa % w(its:ite, jts:jte, kts:kte+1) &
 164                        * xa % w(its:ite, jts:jte, kts:kte+1)))
 165    rms_diff = sqrt(sum(diff(its:ite, jts:jte, kts:kte+1) &
 166                        * diff(its:ite, jts:jte, kts:kte+1)))
 167
 168    if (rms_fild == 0.0) then
 169       write(unit=stdout, fmt='(a)') ' w is zero '
 170    else
 171       write(unit=stdout, fmt='(a,1pe10.4)') ' w RMS error = ', rms_diff
 172       write(unit=stdout, fmt='(a,1pe10.4)') ' w RMS field = ', rms_fild
 173       write(unit=stdout, fmt='(a,1pe10.4)') ' w RMS error/RMS field = ', &
 174          rms_diff/rms_fild
 175    end if
 176
 177    if (trace_use) call da_trace_exit("da_check_xtovpx_errors")
 178
 179 end subroutine da_check_xtovptox_errors
 180
 181