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