Context navigation

Back to Ticket #646

Ticket #646: ropp_1dvar_levmarq_tr1.f90

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

Download in other formats:

Original format