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