da_test_vxtransform.inc
References to this file elsewhere.
1 subroutine da_test_vxtransform( xb, xbx, xp, be, &
2 xa, vv, vp)
3
4 !---------------------------------------------------------------------------
5 ! Purpose: Perform inverse/adjoint tests on control variable transform.
6 !
7 ! Method: 1) Test inverse and adjoint of physical variable transform.
8 ! 2) Test inverse and adjoint of vertical transform.
9 ! 3) Perform adjoint test on complete transform: <x,x> = <v_adj,v>.
10 !
11 !---------------------------------------------------------------------------
12
13 implicit none
14
15 type (xb_type), intent(in) :: xb ! first guess (local).
16 type (xbx_type),intent(in) :: xbx ! Header & non-gridded vars.
17 type (xpose_type), intent(inout) :: xp ! Dimensions and xpose buffers.
18 type (be_type), intent(in) :: be ! background error structure.
19 type (x_type), intent(out) :: xa ! analysis increments (local).
20
21 type (vp_type), intent(out) :: vp ! Test CV structure.
22 type (vp_type), intent(out) :: vv ! Test CV structure.
23
24 real, dimension(1:cv_size) :: cv ! Test control variable.
25 integer :: i
26 type (vp_type) :: vp1, vv1 ! Test CV structure.
27 real, dimension(ims:ime, jms:jme) :: field
28
29
30 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_u, xa2_v, xa2_t, &
31 xa2_p, xa2_q, xa2_rh, xa2_rho, &
32 xa2_qt, xa2_qcw, xa2_qrn
33 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_w
34
35 if (trace_use) call da_trace_entry("da_test_vxtransform")
36
37 allocate (vp1 % v1(ims:ime,jms:jme,kms:kme))
38 allocate (vp1 % v2(ims:ime,jms:jme,kms:kme))
39 allocate (vp1 % v3(ims:ime,jms:jme,kms:kme))
40 allocate (vp1 % v4(ims:ime,jms:jme,kms:kme))
41 allocate (vp1 % v5(ims:ime,jms:jme,kms:kme))
42
43 allocate (vv1 % v1(ims:ime,jms:jme,1:kz_vv(1)))
44 allocate (vv1 % v2(ims:ime,jms:jte,1:kz_vv(2)))
45 allocate (vv1 % v3(ims:ime,jms:jme,1:kz_vv(3)))
46 allocate (vv1 % v4(ims:ime,jms:jme,1:kz_vv(4)))
47 allocate (vv1 % v5(ims:ime,jms:jme,1:kz_vv(5)))
48
49 call da_zero_vp_type(vp1)
50 call da_zero_vp_type(vv1)
51
52 write(unit=stdout, fmt='(/a/)') &
53 'da_testvxtransform:'
54
55 write(unit=stdout, fmt='(/a/)') &
56 '---------------------------------------'
57
58 ! Make cv all constant value 1.0
59
60 call random_number(cv(:))
61 cv(:) = cv(:) - 0.5
62 call da_zero_x(xa)
63
64 if ( global ) then
65 write(unit=stdout, fmt='(/a/)') &
66 ' Inverse tests not performed for global application'
67 go to 1111
68 end if
69
70 #ifdef DM_PARALLEL
71 write(unit=stdout, fmt='(/a/)') &
72 ' Inverse tests will not be done as it is MPP run'
73 #else
74 write(unit=stdout, fmt='(/a/)') &
75 ' Inverse tests follows:'
76 call da_transform_vtox( xb, xbx, be, cv, &
77 vv, vp, xp, xa)
78
79 ! Store xa, Vv & Vp for inverse test
80
81 xa2_u(ims:ime,jms:jme,:) = xa % u(ims:ime,jms:jme,:)
82 xa2_v(ims:ime,jms:jme,:) = xa % v(ims:ime,jms:jme,:)
83 xa2_w(ims:ime,jms:jme,:) = xa % w(ims:ime,jms:jme,:)
84 xa2_t(ims:ime,jms:jme,:) = xa % t(ims:ime,jms:jme,:)
85 xa2_p(ims:ime,jms:jme,:) = xa % p(ims:ime,jms:jme,:)
86 xa2_q(ims:ime,jms:jme,:) = xa % q(ims:ime,jms:jme,:)
87 xa2_rho(ims:ime,jms:jme,:) = xa % rho(ims:ime,jms:jme,:)
88
89 if ( cv_options_hum == 2 ) then
90 xa2_rh(ims:ime,jms:jme,:) = xa % rh(ims:ime,jms:jme,:)
91 end if
92 if ( cv_options_hum == 3 ) then
93 xa2_qt (ims:ime,jms:jme,:) = xa % qt (ims:ime,jms:jme,:)
94 xa2_qcw(ims:ime,jms:jme,:) = xa % qcw(ims:ime,jms:jme,:)
95 xa2_qrn(ims:ime,jms:jme,:) = xa % qrn(ims:ime,jms:jme,:)
96 end if
97
98 vv1 % v1(its:ite,jts:jte,1:kz_vv(1)) = vv % v1(its:ite,jts:jte,1:kz_vv(1))
99 vv1 % v2(its:ite,jts:jte,1:kz_vv(2)) = vv % v2(its:ite,jts:jte,1:kz_vv(2))
100 vv1 % v3(its:ite,jts:jte,1:kz_vv(3)) = vv % v3(its:ite,jts:jte,1:kz_vv(3))
101 vv1 % v4(its:ite,jts:jte,1:kz_vv(4)) = vv % v4(its:ite,jts:jte,1:kz_vv(4))
102 vv1 % v5(its:ite,jts:jte,1:kz_vv(5)) = vv % v5(its:ite,jts:jte,1:kz_vv(5))
103
104 vp1 % v1(its:ite,jts:jte,1:kz_vp(1)) = vp % v1(its:ite,jts:jte,1:kz_vp(1))
105 vp1 % v2(its:ite,jts:jte,1:kz_vp(2)) = vp % v2(its:ite,jts:jte,1:kz_vp(2))
106 vp1 % v3(its:ite,jts:jte,1:kz_vp(3)) = vp % v3(its:ite,jts:jte,1:kz_vp(3))
107 vp1 % v4(its:ite,jts:jte,1:kz_vp(4)) = vp % v4(its:ite,jts:jte,1:kz_vp(4))
108 vp1 % v5(its:ite,jts:jte,1:kz_vp(5)) = vp % v5(its:ite,jts:jte,1:kz_vp(5))
109
110 !----------------------------------------------------------------------
111 ! [1.0]: Perform VvToVpToVv test:
112 !----------------------------------------------------------------------
113 ! call da_transform_xtovp( xb, xbx, xa, xp, vp, be)
114
115 if ( vert_corr == 2 ) then
116 ! perform vv = u_v^{-1} vp transform:
117 call da_vertical_transform( 'u_inv', be, &
118 xb % vertical_inner_product, &
119 vv, vp)
120 else
121 vv % v1(its:ite,jts:jte,1:kz_vv(1)) = vp % v1(its:ite,jts:jte,1:kz_vp(1))
122 vv % v2(its:ite,jts:jte,1:kz_vv(2)) = vp % v2(its:ite,jts:jte,1:kz_vp(2))
123 vv % v3(its:ite,jts:jte,1:kz_vv(3)) = vp % v3(its:ite,jts:jte,1:kz_vp(3))
124 vv % v4(its:ite,jts:jte,1:kz_vv(4)) = vp % v4(its:ite,jts:jte,1:kz_vp(4))
125 vv % v5(its:ite,jts:jte,1:kz_vv(5)) = vp % v5(its:ite,jts:jte,1:kz_vp(5))
126 end if
127
128
129 write(unit=stdout, fmt='(/a/)') &
130 'da_check_vvtovptovv_errors'
131
132 call da_check_vp_errors( vv1, vv, kz_vv, its,ite, jts,jte, kts,kte )
133
134 !----------------------------------------------------------------------
135 ! [2.0]: Perform VpToVvToVp test:
136 !----------------------------------------------------------------------
137 if ( vert_corr == 2 ) then
138 ! perform vp = u_v (vv) transform:
139 call da_vertical_transform( 'u', be, &
140 xb % vertical_inner_product, &
141 vv, vp)
142 else
143
144 vp % v1(its:ite,jts:jte,1:kz_vv(1)) = vv % v1(its:ite,jts:jte,1:kz_vv(1))
145 vp % v2(its:ite,jts:jte,1:kz_vv(2)) = vv % v2(its:ite,jts:jte,1:kz_vv(2))
146 vp % v3(its:ite,jts:jte,1:kz_vv(3)) = vv % v3(its:ite,jts:jte,1:kz_vv(3))
147 vp % v4(its:ite,jts:jte,1:kz_vv(4)) = vv % v4(its:ite,jts:jte,1:kz_vv(4))
148 vp % v5(its:ite,jts:jte,1:kz_vv(5)) = vv % v5(its:ite,jts:jte,1:kz_vv(5))
149
150 end if
151
152 ! Check inverse errors:
153
154 write(unit=stdout, fmt='(/a/)') &
155 'da_check_vptovvtovp_errors'
156
157 call da_check_vp_errors( vp1, vp, kz_vv, its,ite, jts,jte, kts,kte )
158
159 !----------------------------------------------------------------------
160 ! [3.0] Check_CvToVv_Adjoint:
161 !----------------------------------------------------------------------
162 call da_check_cvtovv_adjoint( xb, xbx, xp, be, cv, vv)
163
164 !----------------------------------------------------------------------
165 ! [4.0] Test inverse physical variable transform:
166 ! Note: Currently these test are developed for regional only
167 !----------------------------------------------------------------------
168
169 if( .not. global ) then
170 call da_zero_x(xa)
171 call da_transform_vptox( xb, vp, xp, xa, be
172 ! [4.1] Test XToVpToX differences:
173
174 write(unit=stdout, fmt='(/a/)') &
175 'da_check_xtovptox_errors'
176
177 call da_check_xtovptox_errors( xa, xa2_u, xa2_v, xa2_w, xa2_t, &
178 xa2_p, xa2_q, xa2_rho, &
179 xa2_qt, xa2_qcw, xa2_qrn)
180
181 ! [4.2] Perform v_{p} = U_{p}^{-1} x transform (again):
182
183 call da_transform_xtovp( xb, xbx, xa, xp, vv, be)
184
185 ! [2.4] Check inverse errors:
186
187 write(unit=stdout, fmt='(/a/)') &
188 'da_check_vptoxtovp_errors'
189
190 call da_check_vp_errors( vp, vv, kz_vp, its,ite, jts,jte, kts,kte )
191 end if
192
193 !----------------------------------------------------------------------
194 ! [5.0] Perform Vv -> Vp adjoint test:
195 !----------------------------------------------------------------------
196
197 call da_check_vvtovp_adjoint( xb, xp, be, vp, vv)
198
199 !----------------------------------------------------------------------
200 ! [6.0] Perform Vp -> X adjoint tests:
201 !----------------------------------------------------------------------
202
203 call da_check_vptox_adjoint( xb, be, xa, xp, vp)
204
205 #endif
206 1111 continue
207
208 !----------------------------------------------------------------------
209 ! [7.0]: Perform adjoint test on complete transform: <x,x> = <v_adj,v>
210 !----------------------------------------------------------------------
211
212 call da_check_vtox_adjoint( xb, xbx, be, cv, vv, vp, xp, xa)
213
214 !----------------------------------------------------------------------
215 ! [8.0]: Perform Spectral transform tests for Global
216 !----------------------------------------------------------------------
217
218 if (global) then
219 #ifdef DM_PARALLEL
220 write (unit=stdout,fmt=*) &
221 ' Inverse tests for spectral transforms will not be done as it is MPP run'
222 #else
223 call random_number(field(:,:) )
224 field(:,:) = field(:,:) - 0.5
225 call da_test_spectral(xbx, field)
226 #endif
227 end if
228
229 deallocate (vp1 % v1)
230 deallocate (vp1 % v2)
231 deallocate (vp1 % v3)
232 deallocate (vp1 % v4)
233 deallocate (vp1 % v5)
234
235 deallocate (vv1 % v1)
236 deallocate (vv1 % v2)
237 deallocate (vv1 % v3)
238 deallocate (vv1 % v4)
239 deallocate (vv1 % v5)
240
241 if (trace_use) call da_trace_exit("da_test_vxtransform")
242
243 end subroutine da_test_vxtransform
244
245