da_check_vp_errors.inc
References to this file elsewhere.
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,1:ne)
23
24 if (trace_use) call da_trace_entry("da_check_vp_errors")
25
26 inv_size = 1.0 / real((ite-its+1) * (jte-jts+1) * (kte-kts+1))
27
28 !-------------------------------------------------------------------------
29 ! [1.0]: Check v1 differences:
30 !-------------------------------------------------------------------------
31
32 diff(its:ite,jts:jte, kts:kte) = vp2 % v1(its:ite,jts:jte,kts:kte) - &
33 vp1 % v1(its:ite,jts:jte,kts:kte)
34
35 rms_fild = sqrt(sum(vp1 % v1(its:ite, jts:jte,kts:kte) &
36 * vp1 % v1(its:ite, jts:jte,kts:kte)) * inv_size)
37 rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
38 * diff(its:ite, jts:jte,kts:kte)) * inv_size)
39
40 if (rms_fild == 0.0) then
41 write(unit=stdout, fmt='(a)') ' v1 is zero '
42 else
43 write(unit=stdout, fmt='(a,1pe10.4)') &
44 ' v1 RMS error/RMS field = ', rms_diff/rms_fild
45 end if
46
47 !-------------------------------------------------------------------------
48 ! [2.0]: Check v2 differences:
49 !-------------------------------------------------------------------------
50
51 diff(its:ite,jts:jte, kts:kte) = vp2 % v2(its:ite,jts:jte,kts:kte) - &
52 vp1 % v2(its:ite,jts:jte,kts:kte)
53
54 rms_fild = sqrt(sum(vp1 % v2(its:ite, jts:jte,kts:kte) &
55 * vp1 % v2(its:ite, jts:jte,kts:kte)) * inv_size)
56 rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
57 * diff(its:ite, jts:jte,kts:kte)) * inv_size)
58
59 if (rms_fild == 0.0) then
60 write(unit=stdout, fmt='(a)') ' v2 is zero '
61 else
62 write(unit=stdout, fmt='(a,1pe10.4)') &
63 ' v2 RMS error/RMS field = ', rms_diff/rms_fild
64 end if
65
66 !-------------------------------------------------------------------------
67 ! [3.0]: Check v3 differences:
68 !-------------------------------------------------------------------------
69
70 diff(its:ite,jts:jte, kts:kte) = vp2 % v3(its:ite,jts:jte,kts:kte) - &
71 vp1 % v3(its:ite,jts:jte,kts:kte)
72
73 rms_fild = sqrt(sum(vp1 % v3(its:ite, jts:jte,kts:kte) &
74 * vp1 % v3(its:ite, jts:jte,kts:kte)) * inv_size)
75 rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
76 * diff(its:ite, jts:jte,kts:kte)) * inv_size)
77
78 if (rms_fild == 0.0) then
79 write(unit=stdout, fmt='(a)') ' v3 is zero '
80 else
81 write(unit=stdout, fmt='(a,1pe10.4)') &
82 ' v3 RMS error/RMS field = ', rms_diff/rms_fild
83 end if
84
85 !-------------------------------------------------------------------------
86 ! [4.0]: Check v4 differences:
87 !-------------------------------------------------------------------------
88
89 diff(its:ite,jts:jte, kts:kte) = vp2 % v4(its:ite,jts:jte,kts:kte) - &
90 vp1 % v4(its:ite,jts:jte,kts:kte)
91
92 rms_fild = sqrt(sum(vp1 % v4(its:ite, jts:jte,kts:kte) &
93 * vp1 % v4(its:ite, jts:jte,kts:kte)) * inv_size)
94 rms_diff = sqrt(sum(diff(its:ite, jts:jte,kts:kte) &
95 * diff(its:ite, jts:jte,kts:kte)) * inv_size)
96
97 if (rms_fild == 0.0) then
98 write(unit=stdout, fmt='(a)') ' v4 is zero '
99 else
100 write(unit=stdout, fmt='(a,1pe10.4)') &
101 ' v4 RMS error/RMS field = ', rms_diff/rms_fild
102 end if
103
104 !-------------------------------------------------------------------------
105 ! [5.0]: Check v5 differences:
106 !-------------------------------------------------------------------------
107
108 inv_size = 1.0 / real((ite-its+1) * (jte-jts+1))
109
110 diff(its:ite, jts:jte,kts:kte) = vp2 % v5(its:ite, jts:jte,kts:kte) - &
111 vp1 % v5(its:ite, jts:jte,kts:kte)
112
113 rms_fild = sqrt(sum(vp1 % v5(its:ite, jts:jte,kts:kte) * &
114 vp1 % v5(its:ite, jts:jte,kts:kte)) * inv_size)
115 rms_diff = sqrt(sum(diff(its:ite, jts:jte,1) * &
116 diff(its:ite, jts:jte,1)) * inv_size)
117
118 if (rms_fild == 0.0) then
119 write(unit=stdout, fmt='(a)') ' v5 is zero '
120 else
121 write(unit=stdout, fmt='(a,1pe10.4)') &
122 ' v5 RMS error/RMS field = ', rms_diff/rms_fild
123 end if
124
125 !-------------------------------------------------------------------------
126 ! [6.0]: Check alpha differences:
127 !-------------------------------------------------------------------------
128
129 if (ne > 0) then
130 inv_size = 1.0 / real((ite-its+1) * (jte-jts+1) * ne)
131 diff_alpha(its:ite,jts:jte,1:ne) = vp2 % alpha(its:ite,jts:jte,1:ne) - &
132 vp1 % alpha(its:ite,jts:jte,1:ne)
133
134 rms_fild = sqrt(sum(vp1 % alpha(its:ite, jts:jte,1:ne) &
135 * vp1 % alpha(its:ite, jts:jte,1:ne)) * inv_size)
136 rms_diff = sqrt(sum(diff_alpha(its:ite, jts:jte,1:ne) &
137 * diff_alpha(its:ite, jts:jte,1:ne)) * inv_size)
138
139 if (rms_fild /= 0.0) then
140 write(unit=stdout, fmt='(a,1pe10.4)') ' alpha RMS error/RMS field = ',&
141 rms_diff/rms_fild
142 end if
143 end if
144
145 if (trace_use) call da_trace_exit("da_check_vp_errors")
146
147 end subroutine da_check_vp_errors
148
149