Context navigation

Back to Ticket #646

Ticket #646: ropp_1dvar_levmarq_tr2.f90

File ropp_1dvar_levmarq_tr2.f90, 25.9 KB (added by Ian Culverwell, 6 years ago)
ropp_1dvar_levmarq_tr2.f90

Line
1	! $Id: ropp_1dvar_levmarq.f90 4452 2015-01-29 14:42:02Z idculv $
2
3	!***s 1DVar/ropp_1dvar_levmarq
4	!
5	! NAME
6	! ropp_1dvar_levmarq - Solve the 1DVar for background data using the
7	! Levenberg-Marquardt minimiser
8	!
9	! SYNOPSIS
10	! CALL ropp_1dvar_levmarq(obs, bg, state, config, diag)
11	!
12	!
13	! DESCRIPTION
14	! This subroutine evaluates a quadratic cost function for a
15	! variational data assimilation procedure.
16	!
17	! More specifically, this routine calculates a cost function
18	!
19	! 1 / \| -1 \| \
20	! J = - < y - H(x) \| O \| y - H(x) > +
21	! 2 \ \| \| /
22	!
23	! 1 / \| -1 \| \
24	! - < x - x \| B \| x - x >
25	! 2 \ b \| \| b/
26	!
27	! where the background state x_b is given by the state vector
28	! state.
29	!
30	! A solution for x is obtained by minimising J using the Levenberg-Marquardt
31	! minimisation method.
32	!
33	! INPUTS
34	! obs Observation data structure.
35	! bg Background data structure.
36	! state State vector structure.
37	! config Configuration structure.
38	! diag Diagnostics structure.
39	!
40	! OUTPUT
41	! state
42	! diag
43	!
44	! REFERENCES
45	! Madsen, K., Nielsen, H.B. & Tingleff, O., Methods for Non-Linear Least
46	! Squares Problems, Lyngby, Denmark, 2004. Available at:
47	! http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3215
48	! Nielsen, H.B., Damping Parameter in Marquardtâs Method, Lyngby, Denmark,
49	! 1999. Available at:
50	! http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=648
51	!
52	! AUTHOR
53	! Met Office, Exeter, UK.
54	! Any comments on this software should be given via the ROM SAF
55	! Helpdesk at http://www.romsaf.org
56	!
57	! COPYRIGHT
58	! (c) EUMETSAT. All rights reserved.
59	! For further details please refer to the file COPYRIGHT
60	! which you should have received as part of this distribution.
61	!
62	!****
63
64	!-------------------------------------------------------------------------------
65	! 1. Bending angle
66	!-------------------------------------------------------------------------------
67
68	SUBROUTINE ropp_1dvar_levmarq_bangle(obs, bg, state, config, diag)
69
70	! 1.1 Declarations
71	! ----------------
72
73	USE typesizes, ONLY: wp => EightByteReal
74	USE ropp_utils
75	USE ropp_fm
76	USE ropp_1dvar, not_this => ropp_1dvar_levmarq_bangle
77	USE matrix
78
79	IMPLICIT NONE
80
81	TYPE(Obs1dBangle), INTENT(inout) :: obs ! Observation data
82	TYPE(State1dFM), INTENT(inout) :: bg ! Background data
83	TYPE(State1dFM), INTENT(inout) :: state ! State vector
84	TYPE(VarConfig), INTENT(in) :: config ! Configuration options
85	TYPE(VarDiag), INTENT(inout) :: diag ! Diagnostic output
86
87	REAL(wp) :: J ! Cost function value
88	TYPE(State1dFM) :: x ! Control vector
89	TYPE(State1dFM) :: x_old ! Control for LM parameter testing
90	TYPE(Obs1dBangle) :: y ! Forward model obs
91	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: K ! K-matrix
92	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_x ! Change of state
93	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_y ! Change of observation
94	REAL(wp), DIMENSION(:), ALLOCATABLE :: x_incr ! Increment of the state vector
95	REAL(wp), DIMENSION(:), ALLOCATABLE :: dJ_dx ! Cost function gradient
96	REAL(wp), DIMENSION(:), ALLOCATABLE :: diag_d2J ! Diagonal d2J/dx2
97	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: d2J_dx2 ! 2nd derivative cost fn
98	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: KO ! K O^-1
99	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Bm1 ! B^-1
100	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Om1 ! O^-1
101	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_J ! Change of cost fn
102	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_last ! Previous state vector
103	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_sigma ! State std deviation
104	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_state ! Change of state vector
105	REAL(wp) :: J_last ! Previous cost function
106	REAL(wp) :: J_min ! Minimum cost function
107	INTEGER :: i_pointer ! Value index
108	INTEGER :: i ! Counter
109	INTEGER :: n_iter ! Number of iterations
110	INTEGER :: n_iter_max ! Maximum number of iterations
111	INTEGER :: n_grad ! Number of gradient updates
112	INTEGER :: nobs ! Number of observations
113	INTEGER :: nstate ! No. of state elements
114	REAL(wp) :: mu ! Levenberg-Marquardt factor
115	REAL(wp) :: rho ! Trust region radius
116	REAL(wp) :: beta ! Increment factor for mu
117	REAL(wp) :: gamma ! Decrement factor for mu
118	REAL(wp) :: nu ! Variation factor for mu
119
120	CHARACTER(len = 4) :: it_str
121	CHARACTER(len = 4) :: gr_str
122	CHARACTER(len = 15) :: ch_str, co_str
123	CHARACTER(len = 256) :: routine
124
125	! 1.2 Message handling
126	! --------------------
127
128	CALL message_get_routine(routine)
129	CALL message_set_routine('ropp_1dvar_levmarq')
130
131	! 1.3 Initialise variables
132	! ------------------------
133
134	i_pointer = 0
135	n_grad = 0
136	n_iter = 0
137	n_iter_max = 50
138
139	mu = 1.0E-2
140	beta = 2.0_wp
141	gamma = 3.0_wp
142	nu = beta
143
144	J = 0.0_wp
145	J_last = 1.0E30_wp
146	J_min = J_last
147
148	nstate = SIZE(bg%state)
149	nobs = SIZE(obs%bangle)
150
151	ALLOCATE(K(nobs, nstate))
152	ALLOCATE(delta_x(nstate))
153	ALLOCATE(delta_y(nobs))
154	ALLOCATE(x_incr(nstate))
155	ALLOCATE(dJ_dx(nstate))
156	ALLOCATE(diag_d2J(nstate))
157	ALLOCATE(d2J_dx2(nstate,nstate))
158	ALLOCATE(KO(nstate,nobs))
159	ALLOCATE(Bm1(nstate,nstate))
160	ALLOCATE(Om1(nobs,nobs))
161
162	IF (ALLOCATED(delta_J)) DEALLOCATE(delta_J)
163	ALLOCATE(delta_J(config%conv_check_n_previous))
164	delta_J(:) = 0.0_wp
165
166	IF (ALLOCATED(state_last)) DEALLOCATE(state_last)
167	ALLOCATE(state_last(SIZE(bg%state)))
168	state_last(:) = 0.0_wp
169
170	IF (ALLOCATED(state_sigma)) DEALLOCATE(state_sigma)
171	ALLOCATE(state_sigma(SIZE(bg%state)))
172	DO i = 1, SIZE(state_sigma)
173	state_sigma(i) = SQRT(bg%cov%d(i + i*(i-1)/2)) ! Direct read from matrix_pp
174	ENDDO
175
176	IF (ALLOCATED(delta_state)) DEALLOCATE(delta_state)
177	ALLOCATE(delta_state(config%conv_check_n_previous))
178	delta_state(:) = 0.0_wp
179
180	! 1.4 Inverse error covariances
181	! -----------------------------
182
183	Bm1 = matrix_invert(bg%cov)
184	Om1 = matrix_invert(obs%cov)
185
186	! 1.5 First guess and pseudo-observations
187	! ---------------------------------------
188
189	x = bg ; y = obs
190
191	! 1.6 Initial cost function
192	! -------------------------
193
194	IF (ASSOCIATED(x%ak)) THEN
195	CALL ropp_fm_state2state_ecmwf(x)
196	ELSE
197	CALL ropp_fm_state2state_meto(x)
198	END IF
199
200	CALL ropp_fm_bangle_1d(x, y)
201
202	delta_x = (x%state - bg%state)
203	delta_y = (y%bangle - obs%bangle) * y%weights
204
205	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) &
206	+ 0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
207
208	diag%J_init = J
209	J_last = J
210	J_min = J
211
212	IF (config%minropp%impres == 0) THEN
213	WRITE(it_str, '(i4)') n_iter
214	WRITE(ch_str, '(g15.5)') J
215	ch_str = ADJUSTL(ch_str)
216	co_str = ' - '
217
218	CALL message(msg_cont, &
219	' n_iter = ' // it_str // ' J = ' // ch_str // &
220	' max(relative change in state) = ' // co_str)
221	END IF
222
223	! 1.7 Main minimisation iteration loop
224	! ------------------------------------
225
226	DO WHILE(n_iter < n_iter_max)
227
228	! 1.7.1 Bookkeeping
229
230	n_grad = n_grad + 1
231
232	! 1.7.2 Compute gradient and Hessian of cost function
233
234	CALL ropp_fm_bangle_1d_grad(x, y, K)
235
236	!TODO: (a) Doe we need this? (b) Does it change the cost function?
237	WHERE(ABS(delta_y) > 50.0_wp)
238	delta_y = 0.0_wp
239	END WHERE
240
241	KO = MATMUL(TRANSPOSE(K), Om1)
242	dJ_dx = MATMUL(KO, delta_y) + MATMUL(Bm1, delta_x)
243	d2J_dx2 = Bm1 + MATMUL(KO, K)
244
245	! 1.7.3 Levenberg-Marquardt adjustment of Hessian diagonal
246
247	DO i = 1, nstate
248	d2J_dx2(i,i) = d2J_dx2(i,i) * (1.0_wp + mu)
249	END DO
250
251	! 1.7.4 Solve matrix equation d2J_dx2 . dx = -dJ_dx and update state
252
253	x_incr = matrix_solve(d2J_dx2, - dJ_dx)
254
255	! 1.7.5 Update test state vector
256
257	x_old = x
258
259	IF (x%use_logq) THEN
260	x%state = x%state + SIGN(MIN(ABS(x_incr), ABS(x%state/2.0_wp)), x_incr)
261	ELSE
262	x%state = x%state + x_incr
263	END IF
264
265	! 1.7.6 Special handling of ionospheric state vector elements
266
267	IF (bg%direct_ion) THEN
268
269	i = nstate - 2
270	IF (x%state(i) < ropp_ZERO) THEN
271	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
272	"Ne_max < 0 ... suggest examining final value. \n")
273	END IF
274
275	i = nstate - 1
276	IF (x%state(i) < 0.01_wp*bg%state(i)) THEN
277	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
278	"H_peak < 1% of background ... resetting to background value. \n")
279	x%state(i) = bg%state(i)
280	END IF
281
282	i = nstate
283	IF (x%state(i) < 0.01_wp*bg%state(i)) THEN
284	CALL message(msg_warn, "Levenberg-Marquardt solver returns " // &
285	"H_width < 1% of background ... resetting to background value. \n")
286	x%state(i) = bg%state(i)
287	END IF
288
289	END IF
290
291	! 1.7.3 Compute cost function
292
293	IF (ASSOCIATED(x%ak)) THEN
294	CALL ropp_fm_state2state_ecmwf(x)
295	ELSE
296	CALL ropp_fm_state2state_meto(x)
297	END IF
298
299	CALL ropp_fm_bangle_1d(x, y)
300
301	delta_x = (x%state - bg%state)
302	delta_y = (y%bangle - obs%bangle) * y%weights
303
304	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) &
305	+ 0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
306
307	! 1.8 Levenberg-Marquardt update and convergence check
308	! ----------------------------------------------------
309
310	! Trust region
311
312	rho = 2.0_wp * (J_last - J) / DOT_PRODUCT(x_incr, mu * x_incr - dJ_dx)
313
314	IF (rho > 0.0_wp) THEN ! Check convergence
315
316	! 1.8.4 Update mu
317
318	mu = mu * max(1.0_wp/gamma, 1.0_wp - (beta - 1.0_wp)(2.0_wprho - 1.0_wp)**5)
319	nu = beta
320
321	! 1.8.5 Keep track of current state
322
323	n_iter = n_iter + 1
324	J_last = J
325
326	IF (config % conv_check_apply) THEN
327	i_pointer = MOD(i_pointer + 1, config % conv_check_n_previous)
328	IF (i_pointer == 0) i_pointer = config % conv_check_n_previous
329
330	delta_J(i_pointer) = J_last - J
331	delta_state(i_pointer) = MAXVAL(ABS(state_last - x%state)/state_sigma)
332
333	state_last = x%state
334	END IF
335
336	IF (J < J_min) THEN
337	J_min = J
338	END IF
339
340	! 1.8.6 Check for convergence
341
342	IF (config % conv_check_apply) THEN
343
344	IF (config%minropp%impres == 0) THEN
345	WRITE(it_str, '(i4)') n_iter
346	WRITE(ch_str, '(g15.5)') J
347	ch_str = ADJUSTL(ch_str)
348	IF (n_iter > 0) THEN
349	WRITE(co_str, '(g15.5)') delta_state(i_pointer)
350	co_str = ADJUSTL(co_str)
351	ELSE
352	co_str = ' - '
353	END IF
354
355	CALL message(msg_cont, &
356	' n_iter = ' // it_str // ' J = ' // ch_str // &
357	' max(relative change in state) = ' // co_str)
358	END IF
359
360	IF (MAXVAL(delta_state) < config%conv_check_max_delta_state &
361	.AND. n_iter > config % conv_check_n_previous) THEN
362	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_state
363	ch_str = ADJUSTL(ch_str)
364	WRITE(it_str, '(i2)') config%conv_check_n_previous
365	it_str = ADJUSTL(it_str)
366	CALL message(msg_cont, '')
367	CALL message(msg_info, &
368	'Convergence assumed to be achieved as the state vector did ' // &
369	'not change by more\n ' // 'than ' // TRIM(ch_str) // ' ' // &
370	'relative to the assumed background errors for the last ' // &
371	TRIM(it_str) // ' iterations.\n')
372	EXIT
373	ELSE IF (MAXVAL(ABS(delta_J)) < config%conv_check_max_delta_J &
374	.AND. n_iter > config % conv_check_n_previous) THEN
375	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_J
376	ch_str = ADJUSTL(ch_str)
377	WRITE(it_str, '(i2)') config%conv_check_n_previous
378	it_str = ADJUSTL(it_str)
379	CALL message(msg_cont, '')
380	CALL message(msg_info, &
381	'Convergence assumed to be achieved as the cost function did ' // &
382	'not change by more\n ' // 'than ' // TRIM(ch_str) // &
383	' for the last ' // TRIM(it_str) // ' iterations.\n')
384	EXIT
385	END IF
386
387	END IF ! Convergence check
388
389	ELSE ! rho < 0; retry with old state and updated mu
390
391	x = x_old
392	mu = mu * nu ; mu = 2.0_wp * nu
393
394	!TODO: Add a configuration option for mu update outputs
395	IF (config%minropp%impres == 0) THEN
396	WRITE(it_str, '(i4)') n_iter
397	WRITE(ch_str, '(g15.5)') J
398	ch_str = ADJUSTL(ch_str)
399	WRITE(co_str, '(g15.5)') mu
400	co_str = ADJUSTL(co_str)
401
402	CALL message(msg_cont, &
403	' n_iter = ' // it_str // ' J = ' // ch_str // &
404	' mu -> ' // co_str)
405	END IF
406
407	END IF ! Levenberg-Marquardt update
408	END DO ! end main iteration loop
409
410	IF (n_iter < 50) THEN
411	WRITE(it_str, '(i4)') n_iter ; it_str = ADJUSTL(it_str)
412	WRITE(gr_str, '(i4)') n_grad ; gr_str = ADJUSTL(gr_str)
413	CALL message(msg_info, 'Finished after ' // TRIM(it_str) // ' iterations (' // &
414	TRIM(gr_str) // ' forward model / gradient evaluations).\n')
415	ELSE
416	WRITE(it_str, '(i4)') n_iter ; it_str = ADJUSTL(it_str)
417	CALL message(msg_warn, &
418	'Iteration ended after ' // TRIM(it_str) // &
419	' iterations without convergence achieved.\n')
420	END IF
421
422	! 1.9 Copy solution back to state
423	! -------------------------------
424
425	state = x
426
427	! 1.10 Diagnostic data
428	! --------------------
429
430	diag%J = J_min
431
432	IF (COUNT(obs%weights > 0.0_wp) > 0) THEN
433	diag%J_scaled = 2.0_wp * J_min / REAL(COUNT(obs%weights > 0.0_wp), wp)
434	ENDIF
435
436	diag%n_iter = n_iter
437
438	ALLOCATE (diag%J_bgr(SIZE(state%state)))
439	delta_x = state%state - bg%state
440	diag%J_bgr = 0.5_wp * delta_x * matrix_solve(bg%cov, delta_x)
441
442	! 1.11 Clean up
443	! -------------
444
445	DEALLOCATE(K)
446	DEALLOCATE(delta_x)
447	DEALLOCATE(delta_y)
448	DEALLOCATE(dJ_dx)
449	DEALLOCATE(x_incr)
450	DEALLOCATE(diag_d2J)
451	DEALLOCATE(d2J_dx2)
452	DEALLOCATE(KO)
453	DEALLOCATE(Bm1)
454	DEALLOCATE(Om1)
455
456	CALL message_set_routine(routine)
457
458	END SUBROUTINE ropp_1dvar_levmarq_bangle
459
460	!***
461
462	!-------------------------------------------------------------------------------
463	! 2. Refractivity
464	!-------------------------------------------------------------------------------
465
466	SUBROUTINE ropp_1dvar_levmarq_refrac(obs, bg, state, config, diag)
467
468	! 2.1 Declarations
469	! ----------------
470
471	USE typesizes, ONLY: wp => EightByteReal
472	USE ropp_utils
473	USE ropp_fm
474	USE ropp_1dvar, not_this => ropp_1dvar_levmarq_refrac
475	USE matrix
476
477	IMPLICIT NONE
478
479	TYPE(Obs1dRefrac), INTENT(inout) :: obs ! Observation data
480	TYPE(State1dFM), INTENT(inout) :: bg ! Background data
481	TYPE(State1dFM), INTENT(inout) :: state ! State vector
482	TYPE(VarConfig), INTENT(in) :: config ! Configuration options
483	TYPE(VarDiag), INTENT(inout) :: diag ! Diagnostic output
484
485	REAL(wp) :: J ! Cost function value
486	TYPE(State1dFM) :: x ! Control vector
487	TYPE(State1dFM) :: xmin ! Control at minimum J
488	TYPE(Obs1dRefrac) :: y ! Forward model obs
489	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: K ! K-matrix
490	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_x ! Change of state
491	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_y ! Change of observation
492	REAL(wp), DIMENSION(:), ALLOCATABLE :: dJ_dx ! Cost function gradient
493	REAL(wp), DIMENSION(:), ALLOCATABLE :: diag_d2J ! Diagonal d2J/dx2
494	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: d2J_dx2 ! 2nd derivative cost
495	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: KO ! K O^-1
496	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Bm1 ! B^-1
497	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Om1 ! O^-1
498	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_J ! Change of cost fn
499	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_last ! Previous state vector
500	REAL(wp), DIMENSION(:), ALLOCATABLE :: state_sigma ! State std deviation
501	REAL(wp), DIMENSION(:), ALLOCATABLE :: delta_state ! Change of state vector
502	REAL(wp) :: J_last ! Previous cost function
503	REAL(wp) :: J_min ! Minimum cost function
504	INTEGER :: i_pointer ! Value index
505	INTEGER :: i ! Counter
506	INTEGER :: n_iter ! Number of iterations
507	INTEGER :: nobs ! Number of obserations
508	INTEGER :: nstate ! No. of state elements
509	LOGICAL :: marq ! Flag to L-M minimise
510	REAL(wp) :: lambda ! Iteration factor
511
512	CHARACTER(len = 4) :: it_str
513	CHARACTER(len = 15) :: ch_str, co_str
514	CHARACTER(len = 256) :: routine
515
516	! 2.2 Message handling
517	! --------------------
518
519	CALL message_get_routine(routine)
520	CALL message_set_routine('ropp_1dvar_levmarq')
521
522	! 2.3 Initialise rolling buffers and pointer for convergence checks
523	! -----------------------------------------------------------------
524
525	i_pointer = 0
526	n_iter = 0
527	lambda = 1.0E-4_wp
528	marq = .FALSE.
529	J = 0.0_wp
530	J_last = 1.0E30_wp
531	J_min = J_last
532
533	nstate = SIZE(bg%state)
534	nobs = SIZE(obs%refrac)
535
536	ALLOCATE(K(nobs, nstate))
537	ALLOCATE(delta_x(nstate))
538	ALLOCATE(delta_y(nobs))
539	ALLOCATE(dJ_dx(nstate))
540	ALLOCATE(diag_d2J(nstate))
541	ALLOCATE(d2J_dx2(nstate,nstate))
542	ALLOCATE(KO(nstate,nobs))
543	ALLOCATE(Bm1(nstate,nstate))
544	ALLOCATE(Om1(nobs,nobs))
545
546	Bm1 = matrix_invert(bg%cov)
547	Om1 = matrix_invert(obs%cov)
548
549	IF (ALLOCATED(delta_J)) DEALLOCATE(delta_J)
550	ALLOCATE(delta_J(config%conv_check_n_previous))
551	delta_J(:) = 0.0_wp
552
553	IF (ALLOCATED(state_last)) DEALLOCATE(state_last)
554	ALLOCATE(state_last(SIZE(bg%state)))
555	state_last(:) = 0.0_wp
556
557	IF (ALLOCATED(state_sigma)) DEALLOCATE(state_sigma)
558	ALLOCATE(state_sigma(SIZE(bg%state)))
559	DO i = 1, SIZE(state_sigma)
560	state_sigma(i) = SQRT(bg%cov%d(i + i*(i-1)/2)) ! Direct read from matrix_pp
561	ENDDO
562
563	IF (ALLOCATED(delta_state)) DEALLOCATE(delta_state)
564	ALLOCATE(delta_state(config%conv_check_n_previous))
565	delta_state(:) = 0.0_wp
566
567	! 2.4 First guess state
568	! ---------------------
569
570	x = bg
571	xmin = bg
572
573	! 2.5 Main minimisation iteration loop
574	! ------------------------------------
575
576	DO WHILE (J <= J_min)
577
578	! 2.5.1 Update iteration counter
579
580	n_iter = n_iter + 1
581
582	! 2.6 Compute cost function
583	! -------------------------
584
585	! 2.6.1 Calculate pseudo observations
586	! -----------------------------------
587
588	ALLOCATE(y%refrac(SIZE(obs%refrac)))
589
590	CALL copy_alloc(obs % geop, y % geop)
591	CALL copy_alloc(obs % weights, y % weights)
592
593	! 2.6.2 Forward model
594	! -------------------
595
596	IF(ASSOCIATED(x%ak))THEN
597	CALL ropp_fm_state2state_ecmwf(x)
598	ELSE
599	CALL ropp_fm_state2state_meto(x)
600	ENDIF
601
602	IF (x%new_ref_op) THEN
603	CALL ropp_fm_refrac_1d_new(x, y)
604	ELSE
605	CALL ropp_fm_refrac_1d(x, y)
606	END IF
607
608	! 2.6.3 Compute cost function
609	! ---------------------------
610
611	delta_x = (x%state - bg%state)
612	delta_y = (y%refrac - obs%refrac) * y%weights
613
614	J = 0.5_wp * DOT_PRODUCT(delta_y, MATMUL(Om1, delta_y)) &
615	+ 0.5_wp * DOT_PRODUCT(delta_x, MATMUL(Bm1, delta_x))
616
617	IF (n_iter == 1) diag % J_init = J
618
619	IF (J > J_last) marq = .TRUE.
620
621	! 2.7 Levenberg-Marquardt minimisation
622	! ------------------------------------
623
624	IF(.NOT. marq)THEN
625
626	! 2.7.1 Normal Newtonian iteration
627	! --------------------------------
628
629	lambda = 0.1_wp * lambda
630	marq = .FALSE.
631
632	! 2.7.2 Evaluate K gradient matrix for current x
633	! ----------------------------------------------
634
635	CALL ropp_fm_refrac_1d_grad(x, y, K)
636
637	! 2.7.3 Calculate -dJ_dx vector and d2J_dx2 matrix at x
638	! -----------------------------------------------------
639
640	delta_x = x%state - bg%state
641	delta_y = (obs%refrac - y%refrac) * y%weights
642
643	WHERE(ABS(delta_y) > 50.0_wp)
644	delta_y = 0.0_wp
645	END WHERE
646
647	KO = MATMUL(TRANSPOSE(K), Om1)
648	dJ_dx = MATMUL(KO, delta_y) - MATMUL(Bm1, delta_x)
649	d2J_dx2 = Bm1 + MATMUL(KO, K)
650
651	! 2.7.4 Store inverse of solution covariance matrix
652	! -------------------------------------------------
653
654	DO i=1,SIZE(bg%state)
655	diag_d2J(i) = d2J_dx2(i,i)
656	ENDDO
657
658	ELSE
659
660	! 2.7.5 Levenberg-Marquardt iteration
661	! -----------------------------------
662	! The previous increment increased the value of the penalty function.
663	! Use previous values of -dJ_dx, d2J_dx2 and adjust value of lambda
664
665	marq = .TRUE.
666	lambda = 100.0_wp * lambda
667
668	ENDIF
669
670	! 2.7.6 Levenberg-Marquardt adjustment to diagonal terms
671	! ------------------------------------------------------
672
673	DO i=1, SIZE(bg%state)
674	d2J_dx2(i,i) = diag_d2J(i) * (lambda+1.0_wp)
675	ENDDO
676
677	! 2.7.7 Solve matrix equation d2J_dx2 . dx = -dJ_dx and update state
678	! ------------------------------------------------------------------
679
680	if (x%use_logq) then
681	delta_x = matrix_solve(d2J_dx2, dJ_dx)
682	x%state = x%state + SIGN( MIN(ABS(delta_x),ABS(x%state/2.0_wp)), &
683	delta_x)
684	else
685	x%state = x%state + matrix_solve(d2J_dx2, dJ_dx)
686	endif
687
688	! 2.8 Update state vector variables - model dependent
689	! ---------------------------------
690
691	IF(ASSOCIATED(x%ak))THEN
692	CALL ropp_fm_state2state_ecmwf(x)
693	ELSE
694	CALL ropp_fm_state2state_meto(x)
695	ENDIF
696
697	! 2.9 Keep track of current state
698	! -------------------------------
699
700	IF( J_min > J ) THEN
701	J_min = J
702	xmin = x
703	ENDIF
704
705	IF (config % conv_check_apply) THEN
706
707	i_pointer = MOD(i_pointer + 1, config % conv_check_n_previous)
708	IF (i_pointer == 0) i_pointer = config % conv_check_n_previous
709
710	delta_J(i_pointer) = J_last - J
711	J_last = J
712	delta_state(i_pointer) = MAXVAL(ABS(state_last - x%state)/state_sigma)
713	state_last = x%state
714
715	ENDIF
716
717	! 2.10 Check for convergence
718	! --------------------------
719
720	IF (config % conv_check_apply) THEN
721
722	IF (config%minropp%impres == 0)THEN
723	WRITE(it_str, '(i4)') n_iter
724	WRITE(ch_str, '(g15.5)') J
725	ch_str = ADJUSTL(ch_str)
726	IF (n_iter > 1) THEN
727	WRITE(co_str, '(g15.5)') delta_state(i_pointer)
728	co_str = ADJUSTL(co_str)
729	ELSE
730	co_str = ' - '
731	ENDIF
732
733	CALL message(msg_cont, &
734	' n_iter = ' // it_str // ' J = ' // ch_str // &
735	' max(relative change in state) = ' // co_str)
736	ENDIF
737
738
739	IF (MAXVAL(delta_state) < config%conv_check_max_delta_state &
740	.AND. n_iter > config % conv_check_n_previous) THEN
741	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_state
742	ch_str = ADJUSTL(ch_str)
743	WRITE(it_str, '(i2)') config%conv_check_n_previous
744	it_str = ADJUSTL(it_str)
745	CALL message(msg_cont, '')
746	CALL message(msg_info, &
747	'Convergence assumed to be achieved as the state vector did ' // &
748	'not change by more\n ' // 'than ' // TRIM(ch_str) // ' ' // &
749	'relative to the assumed background errors for the last ' // &
750	TRIM(it_str) // ' iterations.\n')
751	EXIT
752	ELSE IF (MAXVAL(ABS(delta_J)) < config%conv_check_max_delta_J &
753	.AND. n_iter > config % conv_check_n_previous) THEN
754	WRITE(ch_str, '(g15.5)') config%conv_check_max_delta_J
755	ch_str = ADJUSTL(ch_str)
756	WRITE(it_str, '(i2)') config%conv_check_n_previous
757	it_str = ADJUSTL(it_str)
758	CALL message(msg_cont, '')
759	CALL message(msg_info, &
760	'Convergence assumed to be achieved as the cost function did ' // &
761	'not change by more\n ' // 'than ' // TRIM(ch_str) // &
762	' for the last ' // TRIM(it_str) // ' iterations.\n')
763	EXIT
764	ENDIF
765
766	ENDIF
767
768	ENDDO ! end main iteration loop
769
770	! 2.11 Copy solution back to state
771	! --------------------------------
772
773	state = xmin
774
775	! 2.12 Diagnostic data
776	! --------------------
777
778	diag % J = J_min
779
780	IF (COUNT(obs % weights > 0.0_wp) > 0) THEN
781	diag % J_scaled = 2.0_wp * J_min / REAL(COUNT(obs % weights > 0.0_wp), wp)
782	ENDIF
783
784	diag % n_iter = n_iter
785
786	ALLOCATE (diag%J_bgr(SIZE(state%state)))
787	delta_x = state%state - bg%state
788	diag % J_bgr = 0.5_wp * delta_x * matrix_solve(bg%cov, delta_x)
789
790	! 2.13 Clean up
791	! -------------
792
793	DEALLOCATE(K)
794	DEALLOCATE(delta_x)
795	DEALLOCATE(delta_y)
796	DEALLOCATE(dJ_dx)
797	DEALLOCATE(diag_d2J)
798	DEALLOCATE(d2J_dx2)
799	DEALLOCATE(KO)
800	DEALLOCATE(Bm1)
801	DEALLOCATE(Om1)
802
803	CALL message_set_routine(routine)
804
805	END SUBROUTINE ropp_1dvar_levmarq_refrac

Download in other formats:

Original format