da_calculate_j.inc
References to this file elsewhere.
1 subroutine da_calculate_j(it, iter, cv_size, cv_size_jb, cv_size_je, &
2 xbx, be, iv, xhat, cv, &
3 re, y, j, j_grad, &
4 grid, config_flags )
5
6 !---------------------------------------------------------------------------
7 ! Purpose: Initialises the Y-array
8 !---------------------------------------------------------------------------
9
10 implicit none
11
12 integer, intent(in) :: it ! external iteration #.
13 integer, intent(in) :: iter ! internal iteration #.
14 integer, intent(in) :: cv_size ! Total cv size.
15 integer, intent(in) :: cv_size_jb ! Jb cv size.
16 integer, intent(in) :: cv_size_je ! Je cv size.
17 type (xbx_type),intent(in) :: xbx ! For header & non-grid arrays.
18 type (be_type), intent(in) :: be ! background error structure.
19 type (iv_type), intent(inout) :: iv ! innovation vector (o-b).
20 real, intent(in) :: xhat(1:cv_size) ! control variables.
21 real, intent(in) :: cv(1:cv_size) ! control variables.
22 type (y_type) , intent(inout) :: re ! residual vector (o-a).
23 type (y_type) , intent(inout) :: y ! y = H(x_inc).
24 real, intent(out) :: j_grad(1:cv_size) ! control variables.
25 type (j_type) , intent(out) :: j ! cost function j
26
27 type(domain), intent(inout) :: grid
28 type(grid_config_rec_type), intent(inout) :: config_flags
29
30 integer :: je_start, je_end ! Start/end indices of Je.
31 real :: jo_partial ! jo for this processor
32 real :: gnorm_jb, gnorm_jo, gnorm_j, gnorm_je
33 type (y_type) :: jo_grad_y ! Grad_y(jo)
34 real :: cv_jb(1:cv_size_jb) ! Jb control variable.
35 real :: xhat_jb(1:cv_size_jb) ! Jb control variable.
36 real :: cv_je(1:cv_size_je) ! Je control variable.
37 real :: xhat_je(1:cv_size_je) ! Je control variable.
38 real, allocatable, dimension(:) :: cv_xhat_jb, cv_xhat_je
39 integer :: ndynopt
40
41
42 if (trace_use) call da_trace_entry("da_calculate_j")
43
44 allocate(cv_xhat_jb(cv_size_jb))
45 allocate(cv_xhat_je(cv_size_je))
46
47
48 je_start = cv_size_jb + 1
49 je_end = cv_size_jb + cv_size_je
50
51 call da_allocate_y(iv, jo_grad_y)
52
53 !-------------------------------------------------------------------------
54 ! [1.0] calculate jo:
55 !-------------------------------------------------------------------------
56
57 ! [1.1] transform from control variable to model grid space:
58
59 if (iter > 0) &
60 call da_transform_vtoy(cv_size, be, grid%ep, xhat, iv, grid%vp, grid%vv,&
61 xbx, y, &
62 grid, config_flags )
63
64 ! [1.2] compute residual (o-a) = (o-b) - h x~
65
66 call da_calculate_residual(iv, y, re)
67
68 ! [1.3] calculate jo:
69
70 call da_jo_and_grady(iv, re, jo_partial, j % jo, jo_grad_y)
71
72 if (test_dm_exact) then
73 ! jo_partial has been already summed at lower level
74 j % jo % total = jo_partial
75 else
76 j % jo % total = wrf_dm_sum_real(jo_partial)
77 end if
78
79 ! [1.4] calculate jc-dfi:
80
81 j % jc = 0.0
82
83 if ( var4d .and. iter > 0 ) then
84
85 ndynopt = grid%dyn_opt
86 grid%dyn_opt = DYN_EM_TL
87 call nl_set_dyn_opt (1 , DYN_EM_TL)
88
89 call da_med_initialdata_input(grid , config_flags, 'tldf')
90
91 grid%dyn_opt = ndynopt
92 call nl_set_dyn_opt (1 , DYN_EM)
93
94 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot(size(grid%g_u_2), grid%g_u_2, grid%g_u_2)
95 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot(size(grid%g_v_2), grid%g_v_2, grid%g_v_2)
96 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot(size(grid%g_t_2), grid%g_t_2, grid%g_t_2)
97 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot(size(grid%g_mu_2), grid%g_mu_2, grid%g_mu_2)
98 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot(size(grid%g_moist), grid%g_moist, grid%g_moist)
99
100 end if
101
102 !-------------------------------------------------------------------------
103 ! [2.0] calculate jb:
104 !-------------------------------------------------------------------------
105
106 cv_jb(1:cv_size_jb) = cv (1:cv_size_jb)
107 xhat_jb(1:cv_size_jb) = xhat (1:cv_size_jb)
108
109 cv_xhat_jb(1:cv_size_jb) = cv_jb(1:cv_size_jb) + xhat_jb(1:cv_size_jb)
110
111 j % jb = 0.5 * da_dot_cv(cv_size_jb, cv_size_domain_jb, &
112 cv_xhat_jb, cv_xhat_jb, grid, &
113 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
114
115 j % jb = jb_factor * j % jb
116
117 !-------------------------------------------------------------------------
118 ! [3.0] calculate je:
119 !-------------------------------------------------------------------------
120
121 j % je = 0.0
122 if (be % ne > 0) then
123 cv_je(1:cv_size_je) = cv(je_start:je_end)
124 xhat_je(1:cv_size_je) = xhat(je_start:je_end)
125 cv_xhat_je(1:cv_size_je) = cv_je(1:cv_size_je) + xhat_je(1:cv_size_je)
126 j % je = 0.5 * da_dot_cv(cv_size_je, cv_size_domain_je, &
127 cv_xhat_je, cv_xhat_je, grid, &
128 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
129 end if
130 j % je = je_factor * j % je
131
132 !-------------------------------------------------------------------------
133 ! [4.0] calculate total cost function j = jo + jb + jc + je:
134 !-------------------------------------------------------------------------
135
136 if (grid%jcdfi_use) then
137 j % total = j % jb + j % jo % total + j % jc + j % je
138 else
139 j % total = j % jb + j % jo % total + j % je
140 end if
141
142 if (it == 1 .and. iter == 0 .and. rootproc) then
143 write(unit=cost_unit,fmt='(a)')'Outer EPS Inner J Jb Jo Jc Je'
144 write(unit=cost_unit,fmt='(a)')'Iter Iter '
145 write(unit=grad_unit,fmt='(a)')'Outer EPS Inner G Gb Go Ge'
146 write(unit=grad_unit,fmt='(a)')'Iter Iter '
147 end if
148
149 if (rootproc) then
150 write(unit=cost_unit,fmt='(2x,i2,1x,e10.3,2x,i4,5(1x,f10.3))') &
151 it, EPS(it), iter, j % total, j % jb, j % jo % total, j % jc, j % je
152 end if
153
154 !-------------------------------------------------------------------------
155 ! [5.0] calculate grad_v (jo):
156 !-------------------------------------------------------------------------
157
158 call da_transform_vtoy_adj(iter, cv_size, be, grid%ep, j_grad, iv, &
159 grid%vp, grid%vv, xbx, jo_grad_y, &
160 grid, config_flags)
161
162 call da_deallocate_y(jo_grad_y)
163
164 if (print_detail_grad) then
165 write(unit=stdout, fmt='(a,i6)') 'Calculate grad_v(jo) iter=', iter
166 write(unit=stdout, fmt='(a, e24.14)') &
167 ' cv_jb.cv_jb = ',sum(cv_jb(:) * cv_jb(:))
168 write(unit=stdout, fmt='(a, e24.14)') &
169 ' cv_je.cv_je = ',sum(cv_je(:) * cv_je(:))
170 write(unit=stdout, fmt='(a, e24.14)') &
171 ' xhat.xhat = ',sum(xhat(1:cv_size) * xhat(1:cv_size))
172 write(unit=stdout, fmt='(a, e24.14)') &
173 ' j_grad.j_grad = ',sum(j_grad(1:cv_size) * j_grad(1:cv_size))
174 end if
175
176 !-------------------------------------------------------------------------
177 ! [6.0] calculate grad_v (j) = grad_v (jb) + grad_v (jo) + grad_v (je):
178 !-------------------------------------------------------------------------
179
180 gnorm_jo = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, grid, &
181 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
182
183 ! Add Del_v (Jb)
184 j_grad(1:cv_size_jb) = jb_factor * &
185 (cv(1:cv_size_jb) + xhat(1:cv_size_jb)) + &
186 j_grad(1:cv_size_jb)
187
188 ! Add Del_v (Je)
189 if (cv_size_je > 0) then
190 j_grad(je_start:je_end) = je_factor * &
191 (cv(je_start:je_end) + xhat(je_start:je_end)) + &
192 j_grad(je_start:je_end)
193 end if
194 cv_xhat_jb(1:cv_size_jb) = cv_jb(1:cv_size_jb) + xhat_jb(1:cv_size_jb)
195 gnorm_jb = da_dot_cv(cv_size_jb, cv_size_domain_jb, &
196 cv_xhat_jb, cv_xhat_jb, grid, &
197 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
198 cv_xhat_je(1:cv_size_je) = cv_je(1:cv_size_je) + xhat_je(1:cv_size_je)
199 gnorm_je = da_dot_cv(cv_size_je, cv_size_domain_je, &
200 cv_xhat_je, cv_xhat_je, grid, &
201 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
202
203 if (print_detail_grad) then
204 write(unit=stdout, fmt='(a,i6)') 'Calculate grad_v(j) iter=', iter
205 write(unit=stdout, fmt='(a, e24.14)') &
206 ' cv_jb.cv_jb = ', sum(cv_jb(1:cv_size_jb) * cv_jb(1:cv_size_jb))
207 write(unit=stdout, fmt='(a, e24.14)') &
208 ' cv_je.cv_je = ', sum(cv_je(1:cv_size_je) * cv_je(1:cv_size_je))
209 write(unit=stdout, fmt='(a, e24.14)') &
210 ' xhat_jb.xhat_jb = ', sum(xhat_jb(1:cv_size_jb) * xhat_jb(1:cv_size_jb))
211 write(unit=stdout, fmt='(a, e24.14)') &
212 ' xhat_je.xhat_je = ', sum(xhat_je(1:cv_size_je) * xhat_je(1:cv_size_je))
213 write(unit=stdout, fmt='(a, e24.14)') &
214 ' j_grad.j_grad = ',sum(j_grad(1:cv_size) * j_grad(1:cv_size))
215 write (unit=stdout,fmt='(A)') " "
216 end if
217
218 gnorm_j = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, grid, &
219 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
220
221 gnorm_jo= sqrt(gnorm_jo)
222 gnorm_jb= jb_factor * sqrt(gnorm_jb)
223 gnorm_je= je_factor * sqrt(gnorm_je)
224 gnorm_j = sqrt(gnorm_j)
225
226 if (rootproc) then
227 write(grad_unit,fmt='(2x,i2,1x,e10.3,2x,i4,5(1x,f10.3))') &
228 it, eps(it), iter, gnorm_j, gnorm_jb, gnorm_jo, gnorm_je
229 end if
230
231 deallocate(cv_xhat_jb)
232 deallocate(cv_xhat_je)
233
234 if (trace_use) call da_trace_exit("da_calculate_j")
235
236 end subroutine da_calculate_j
237
238