Context navigation

Back to Ticket #646

Ticket #646: ropp_1dvar_levmarq_nr.f90

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

Download in other formats:

Original format