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, is, ie, js, je, ks, ke
40 real :: dtemp1x
41
42
43 if (trace_use) call da_trace_entry("da_calculate_j")
44
45 allocate(cv_xhat_jb(cv_size_jb))
46 allocate(cv_xhat_je(cv_size_je))
47
48
49 je_start = cv_size_jb + 1
50 je_end = cv_size_jb + cv_size_je
51
52 call da_allocate_y(iv, jo_grad_y)
53
54 !-------------------------------------------------------------------------
55 ! [1.0] calculate jo:
56 !-------------------------------------------------------------------------
57
58 ! [1.1] transform from control variable to model grid space:
59
60 if (iter > 0) &
61 call da_transform_vtoy(cv_size, be, grid%ep, xhat, iv, grid%vp, grid%vv,&
62 xbx, y, &
63 grid, config_flags )
64
65 ! [1.2] compute residual (o-a) = (o-b) - h x~
66
67 call da_calculate_residual(iv, y, re)
68
69 ! [1.3] calculate jo:
70
71 call da_jo_and_grady(iv, re, jo_partial, j % jo, jo_grad_y)
72
73 if (test_dm_exact) then
74 ! jo_partial has been already summed at lower level
75 j % jo % total = jo_partial
76 else
77 j % jo % total = wrf_dm_sum_real(jo_partial)
78 end if
79
80 ! [1.4] calculate jc-dfi:
81
82 j % jc = 0.0
83
84 if ( var4d .and. iter > 0 ) then
85
86 ndynopt = grid%dyn_opt
87 grid%dyn_opt = DYN_EM_TL
88 call nl_set_dyn_opt (1 , DYN_EM_TL)
89
90 call da_med_initialdata_input(grid , config_flags, 'tldf')
91
92 grid%dyn_opt = ndynopt
93 call nl_set_dyn_opt (1 , DYN_EM)
94
95 is = grid%xp%ips
96 ie = grid%xp%ipe
97 if ( ie == grid%xp%ide ) ie = ie + 1
98 js = grid%xp%jps
99 je = grid%xp%jpe
100 if ( je == grid%xp%jde ) je = je + 1
101 ks = grid%xp%kps
102 ke = grid%xp%kpe
103
104 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_u_2(is:ie,js:je,ks:ke), grid%g_u_2(is:ie,js:je,ks:ke) )
105 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_v_2(is:ie,js:je,ks:ke), grid%g_v_2(is:ie,js:je,ks:ke) )
106 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_t_2(is:ie,js:je,ks:ke), grid%g_t_2(is:ie,js:je,ks:ke) )
107 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1), grid%g_mu_2(is:ie,js:je), grid%g_mu_2(is:ie,js:je) )
108 j % jc = j % jc + 0.5*grid%jcdfi_gama*da_dot( (ie-is+1)*(je-js+1)*(ke-ks+1), grid%g_moist(is:ie,js:je,ks:ke,P_G_QV), grid%g_moist(is:ie,js:je,ks:ke,P_G_QV) )
109
110 dtemp1x = j % jc
111 ! summation across processors:
112 j % jc = wrf_dm_sum_real(dtemp1x)
113
114 end if
115
116 !-------------------------------------------------------------------------
117 ! [2.0] calculate jb:
118 !-------------------------------------------------------------------------
119
120 cv_jb(1:cv_size_jb) = cv (1:cv_size_jb)
121 xhat_jb(1:cv_size_jb) = xhat (1:cv_size_jb)
122
123 cv_xhat_jb(1:cv_size_jb) = cv_jb(1:cv_size_jb) + xhat_jb(1:cv_size_jb)
124
125 j % jb = 0.5 * da_dot_cv(cv_size_jb, cv_size_domain_jb, &
126 cv_xhat_jb, cv_xhat_jb, grid, &
127 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
128
129 j % jb = jb_factor * j % jb
130
131 !-------------------------------------------------------------------------
132 ! [3.0] calculate je:
133 !-------------------------------------------------------------------------
134
135 j % je = 0.0
136 if (be % ne > 0) then
137 cv_je(1:cv_size_je) = cv(je_start:je_end)
138 xhat_je(1:cv_size_je) = xhat(je_start:je_end)
139 cv_xhat_je(1:cv_size_je) = cv_je(1:cv_size_je) + xhat_je(1:cv_size_je)
140 j % je = 0.5 * da_dot_cv(cv_size_je, cv_size_domain_je, &
141 cv_xhat_je, cv_xhat_je, grid, &
142 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
143 end if
144 j % je = je_factor * j % je
145
146 !-------------------------------------------------------------------------
147 ! [4.0] calculate total cost function j = jo + jb + jc + je:
148 !-------------------------------------------------------------------------
149
150 if (grid%jcdfi_use) then
151 j % total = j % jb + j % jo % total + j % jc + j % je
152 else
153 j % total = j % jb + j % jo % total + j % je
154 end if
155
156 if (it == 1 .and. iter == 0 .and. rootproc) then
157 write(unit=cost_unit,fmt='(a)')'Outer EPS Inner J Jb Jo Jc Je'
158 write(unit=cost_unit,fmt='(a)')'Iter Iter '
159 write(unit=grad_unit,fmt='(a)')'Outer EPS Inner G Gb Go Ge'
160 write(unit=grad_unit,fmt='(a)')'Iter Iter '
161 end if
162
163 if (rootproc) then
164 write(unit=cost_unit,fmt='(2x,i2,1x,e10.3,2x,i4,5(1x,f10.3))') &
165 it, EPS(it), iter, j % total, j % jb, j % jo % total, j % jc, j % je
166 end if
167
168 !-------------------------------------------------------------------------
169 ! [5.0] calculate grad_v (jo):
170 !-------------------------------------------------------------------------
171
172 call da_transform_vtoy_adj(iter, cv_size, be, grid%ep, j_grad, iv, &
173 grid%vp, grid%vv, xbx, jo_grad_y, &
174 grid, config_flags, .false.)
175
176 call da_deallocate_y(jo_grad_y)
177
178 if (print_detail_grad) then
179 write(unit=stdout, fmt='(a,i6)') 'Calculate grad_v(jo) iter=', iter
180 write(unit=stdout, fmt='(a, e24.14)') &
181 ' cv_jb.cv_jb = ',sum(cv_jb(:) * cv_jb(:))
182 write(unit=stdout, fmt='(a, e24.14)') &
183 ' cv_je.cv_je = ',sum(cv_je(:) * cv_je(:))
184 write(unit=stdout, fmt='(a, e24.14)') &
185 ' xhat.xhat = ',sum(xhat(1:cv_size) * xhat(1:cv_size))
186 write(unit=stdout, fmt='(a, e24.14)') &
187 ' j_grad.j_grad = ',sum(j_grad(1:cv_size) * j_grad(1:cv_size))
188 end if
189
190 !-------------------------------------------------------------------------
191 ! [6.0] calculate grad_v (j) = grad_v (jb) + grad_v (jo) + grad_v (je):
192 !-------------------------------------------------------------------------
193
194 gnorm_jo = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, grid, &
195 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
196
197 ! Add Del_v (Jb)
198 j_grad(1:cv_size_jb) = jb_factor * &
199 (cv(1:cv_size_jb) + xhat(1:cv_size_jb)) + &
200 j_grad(1:cv_size_jb)
201
202 ! Add Del_v (Je)
203 if (cv_size_je > 0) then
204 j_grad(je_start:je_end) = je_factor * &
205 (cv(je_start:je_end) + xhat(je_start:je_end)) + &
206 j_grad(je_start:je_end)
207 end if
208 cv_xhat_jb(1:cv_size_jb) = cv_jb(1:cv_size_jb) + xhat_jb(1:cv_size_jb)
209 gnorm_jb = da_dot_cv(cv_size_jb, cv_size_domain_jb, &
210 cv_xhat_jb, cv_xhat_jb, grid, &
211 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
212 cv_xhat_je(1:cv_size_je) = cv_je(1:cv_size_je) + xhat_je(1:cv_size_je)
213 gnorm_je = da_dot_cv(cv_size_je, cv_size_domain_je, &
214 cv_xhat_je, cv_xhat_je, grid, &
215 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
216
217 if (print_detail_grad) then
218 write(unit=stdout, fmt='(a,i6)') 'Calculate grad_v(j) iter=', iter
219 write(unit=stdout, fmt='(a, e24.14)') &
220 ' cv_jb.cv_jb = ', sum(cv_jb(1:cv_size_jb) * cv_jb(1:cv_size_jb))
221 write(unit=stdout, fmt='(a, e24.14)') &
222 ' cv_je.cv_je = ', sum(cv_je(1:cv_size_je) * cv_je(1:cv_size_je))
223 write(unit=stdout, fmt='(a, e24.14)') &
224 ' xhat_jb.xhat_jb = ', sum(xhat_jb(1:cv_size_jb) * xhat_jb(1:cv_size_jb))
225 write(unit=stdout, fmt='(a, e24.14)') &
226 ' xhat_je.xhat_je = ', sum(xhat_je(1:cv_size_je) * xhat_je(1:cv_size_je))
227 write(unit=stdout, fmt='(a, e24.14)') &
228 ' j_grad.j_grad = ',sum(j_grad(1:cv_size) * j_grad(1:cv_size))
229 write (unit=stdout,fmt='(A)') " "
230 end if
231
232 gnorm_j = da_dot_cv(cv_size, cv_size_domain, j_grad, j_grad, grid, &
233 (/ be%v1%mz, be%v2%mz, be%v3%mz, be%v4%mz, be%v5%mz, be%alpha%mz /))
234
235 gnorm_jo= sqrt(gnorm_jo)
236 gnorm_jb= jb_factor * sqrt(gnorm_jb)
237 gnorm_je= je_factor * sqrt(gnorm_je)
238 gnorm_j = sqrt(gnorm_j)
239
240 if (rootproc) then
241 write(grad_unit,fmt='(2x,i2,1x,e10.3,2x,i4,5(1x,f10.3))') &
242 it, eps(it), iter, gnorm_j, gnorm_jb, gnorm_jo, gnorm_je
243 end if
244
245 deallocate(cv_xhat_jb)
246 deallocate(cv_xhat_je)
247
248 if (trace_use) call da_trace_exit("da_calculate_j")
249
250 end subroutine da_calculate_j
251
252