ROMS
Loading...
Searching...
No Matches
posterior.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE posterior_mod
3
4#if defined WEAK_CONSTRAINT && \
5 (defined posterior_eofs || defined posterior_error_i || \
6 defined posterior_error_f)
7!
8!git $Id$
9!================================================== Hernan G. Arango ===
10! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
11! Licensed under a MIT/X style license !
12! See License_ROMS.md !
13!=======================================================================
14! !
15! This module computes the Lanczos vectors and eigenvectors of the !
16! posterior analysis error covariance matrix. !
17! !
18!=======================================================================
19!
20 USE mod_param
21 USE mod_parallel
22# ifdef ADJUST_BOUNDARY
23 USE mod_boundary
24# endif
25# ifdef SOLVE3D
26 USE mod_coupling
27# endif
28# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
29 USE mod_forces
30# endif
31 USE mod_fourdvar
32 USE mod_grid
33 USE mod_iounits
34 USE mod_ncparam
35 USE mod_netcdf
36 USE mod_ocean
37# if defined PIO_LIB && defined DISTRIBUTE
38 USE mod_pio_netcdf
39# endif
40 USE mod_stepping
41 USE mod_scalars
42!
43# ifdef DISTRIBUTE
44 USE distribute_mod, ONLY : mp_bcastf, mp_bcasti
45# endif
46 USE lapack_mod, ONLY : dsteqr
47 USE state_addition_mod, ONLY : state_addition
48 USE state_copy_mod, ONLY : state_copy
49 USE state_dotprod_mod, ONLY : state_dotprod
50 USE state_initialize_mod, ONLY : state_initialize
51 USE state_read_mod, ONLY : state_read
52 USE state_scale_mod, ONLY : state_scale
53 USE strings_mod, ONLY : founderror
54 USE wrt_hessian_mod, ONLY : wrt_hessian
55!
56 implicit none
57!
58 PRIVATE
59 PUBLIC :: posterior
60!
61 CONTAINS
62!
63!***********************************************************************
64 SUBROUTINE posterior (ng, tile, model, innLoop, outLoop, Ltrace)
65!***********************************************************************
66!
67! Imported variable declarations.
68!
69 integer, intent(in) :: ng, tile, model, innloop, outloop
70 logical, intent(in) :: ltrace
71!
72! Local variable declarations.
73!
74 character (len=*), parameter :: myfile = &
75 & __FILE__
76!
77# include "tile.h"
78!
79# ifdef PROFILE
80 CALL wclock_on (ng, model, 83, __line__, myfile)
81# endif
82!
83 CALL posterior_tile (ng, tile, model, &
84 & lbi, ubi, lbj, ubj, lbij, ubij, &
85 & imins, imaxs, jmins, jmaxs, &
86 & lold(ng), lnew(ng), &
87 & innloop, outloop, ltrace, &
88# ifdef MASKING
89 & grid(ng) % rmask, &
90 & grid(ng) % umask, &
91 & grid(ng) % vmask, &
92# endif
93# ifdef ADJUST_BOUNDARY
94# ifdef SOLVE3D
95 & boundary(ng) % t_obc, &
96 & boundary(ng) % u_obc, &
97 & boundary(ng) % v_obc, &
98# endif
99 & boundary(ng) % ubar_obc, &
100 & boundary(ng) % vbar_obc, &
101 & boundary(ng) % zeta_obc, &
102# endif
103# ifdef ADJUST_WSTRESS
104 & forces(ng) % ustr, &
105 & forces(ng) % vstr, &
106# endif
107# ifdef SOLVE3D
108# ifdef ADJUST_STFLUX
109 & forces(ng) % tflux, &
110# endif
111 & ocean(ng) % t, &
112 & ocean(ng) % u, &
113 & ocean(ng) % v, &
114# if defined WEAK_CONSTRAINT && defined TIME_CONV
115 & ocean(ng) % ubar, &
116 & ocean(ng) % vbar, &
117# endif
118# else
119 & ocean(ng) % ubar, &
120 & ocean(ng) % vbar, &
121# endif
122 & ocean(ng) % zeta, &
123# ifdef ADJUST_BOUNDARY
124# ifdef SOLVE3D
125 & boundary(ng) % tl_t_obc, &
126 & boundary(ng) % tl_u_obc, &
127 & boundary(ng) % tl_v_obc, &
128# endif
129 & boundary(ng) % tl_ubar_obc, &
130 & boundary(ng) % tl_vbar_obc, &
131 & boundary(ng) % tl_zeta_obc, &
132# endif
133# ifdef ADJUST_WSTRESS
134 & forces(ng) % tl_ustr, &
135 & forces(ng) % tl_vstr, &
136# endif
137# ifdef SOLVE3D
138# ifdef ADJUST_STFLUX
139 & forces(ng) % tl_tflux, &
140# endif
141 & ocean(ng) % tl_t, &
142 & ocean(ng) % tl_u, &
143 & ocean(ng) % tl_v, &
144# if defined WEAK_CONSTRAINT && defined TIME_CONV
145 & ocean(ng) % tl_ubar, &
146 & ocean(ng) % tl_vbar, &
147# endif
148# else
149 & ocean(ng) % tl_ubar, &
150 & ocean(ng) % tl_vbar, &
151# endif
152 & ocean(ng) % tl_zeta, &
153# ifdef ADJUST_BOUNDARY
154# ifdef SOLVE3D
155 & boundary(ng) % d_t_obc, &
156 & boundary(ng) % d_u_obc, &
157 & boundary(ng) % d_v_obc, &
158# endif
159 & boundary(ng) % d_ubar_obc, &
160 & boundary(ng) % d_vbar_obc, &
161 & boundary(ng) % d_zeta_obc, &
162# endif
163# ifdef ADJUST_WSTRESS
164 & forces(ng) % d_sustr, &
165 & forces(ng) % d_svstr, &
166# endif
167# ifdef SOLVE3D
168# ifdef ADJUST_STFLUX
169 & forces(ng) % d_stflx, &
170# endif
171 & ocean(ng) % d_t, &
172 & ocean(ng) % d_u, &
173 & ocean(ng) % d_v, &
174# if defined WEAK_CONSTRAINT && defined TIME_CONV
175 & ocean(ng) % d_ubar, &
176 & ocean(ng) % d_vbar, &
177# endif
178# else
179 & ocean(ng) % d_ubar, &
180 & ocean(ng) % d_vbar, &
181# endif
182 & ocean(ng) % d_zeta, &
183# ifdef ADJUST_BOUNDARY
184# ifdef SOLVE3D
185 & boundary(ng) % ad_t_obc, &
186 & boundary(ng) % ad_u_obc, &
187 & boundary(ng) % ad_v_obc, &
188# endif
189 & boundary(ng) % ad_ubar_obc, &
190 & boundary(ng) % ad_vbar_obc, &
191 & boundary(ng) % ad_zeta_obc, &
192# endif
193# ifdef ADJUST_WSTRESS
194 & forces(ng) % ad_ustr, &
195 & forces(ng) % ad_vstr, &
196# endif
197# ifdef SOLVE3D
198# ifdef ADJUST_STFLUX
199 & forces(ng) % ad_tflux, &
200# endif
201 & ocean(ng) % ad_t, &
202 & ocean(ng) % ad_u, &
203 & ocean(ng) % ad_v, &
204# if defined WEAK_CONSTRAINT && defined TIME_CONV
205 & ocean(ng) % ad_ubar, &
206 & ocean(ng) % ad_vbar, &
207# endif
208# else
209 & ocean(ng) % ad_ubar, &
210 & ocean(ng) % ad_vbar, &
211# endif
212 & ocean(ng) % ad_zeta)
213# ifdef PROFILE
214 CALL wclock_off (ng, model, 83, __line__, myfile)
215# endif
216!
217 RETURN
218 END SUBROUTINE posterior
219!
220!***********************************************************************
221 SUBROUTINE posterior_tile (ng, tile, model, &
222 & LBi, UBi, LBj, UBj, LBij, UBij, &
223 & IminS, ImaxS, JminS, JmaxS, &
224 & Lold, Lnew, &
225 & innLoop, outLoop, Ltrace, &
226# ifdef MASKING
227 & rmask, umask, vmask, &
228# endif
229# ifdef ADJUST_BOUNDARY
230# ifdef SOLVE3D
231 & nl_t_obc, nl_u_obc, nl_v_obc, &
232# endif
233 & nl_ubar_obc, nl_vbar_obc, &
234 & nl_zeta_obc, &
235# endif
236# ifdef ADJUST_WSTRESS
237 & nl_ustr, nl_vstr, &
238# endif
239# ifdef SOLVE3D
240# ifdef ADJUST_STFLUX
241 & nl_tflux, &
242# endif
243 & nl_t, nl_u, nl_v, &
244# if defined WEAK_CONSTRAINT && defined TIME_CONV
245 & nl_ubar, nl_vbar, &
246# endif
247# else
248 & nl_ubar, nl_vbar, &
249# endif
250 & nl_zeta, &
251# ifdef ADJUST_BOUNDARY
252# ifdef SOLVE3D
253 & tl_t_obc, tl_u_obc, tl_v_obc, &
254# endif
255 & tl_ubar_obc, tl_vbar_obc, &
256 & tl_zeta_obc, &
257# endif
258# ifdef ADJUST_WSTRESS
259 & tl_ustr, tl_vstr, &
260# endif
261# ifdef SOLVE3D
262# ifdef ADJUST_STFLUX
263 & tl_tflux, &
264# endif
265 & tl_t, tl_u, tl_v, &
266# if defined WEAK_CONSTRAINT && defined TIME_CONV
267 & tl_ubar, tl_vbar, &
268# endif
269# else
270 & tl_ubar, tl_vbar, &
271# endif
272 & tl_zeta, &
273# ifdef ADJUST_BOUNDARY
274# ifdef SOLVE3D
275 & d_t_obc, d_u_obc, d_v_obc, &
276# endif
277 & d_ubar_obc, d_vbar_obc, &
278 & d_zeta_obc, &
279# endif
280# ifdef ADJUST_WSTRESS
281 & d_sustr, d_svstr, &
282# endif
283# ifdef SOLVE3D
284# ifdef ADJUST_STFLUX
285 & d_stflx, &
286# endif
287 & d_t, d_u, d_v, &
288# if defined WEAK_CONSTRAINT && defined TIME_CONV
289 & d_ubar, d_vbar, &
290# endif
291# else
292 & d_ubar, d_vbar, &
293# endif
294 & d_zeta, &
295# ifdef ADJUST_BOUNDARY
296# ifdef SOLVE3D
297 & ad_t_obc, ad_u_obc, ad_v_obc, &
298# endif
299 & ad_ubar_obc, ad_vbar_obc, &
300 & ad_zeta_obc, &
301# endif
302# ifdef ADJUST_WSTRESS
303 & ad_ustr, ad_vstr, &
304# endif
305# ifdef SOLVE3D
306# ifdef ADJUST_STFLUX
307 & ad_tflux, &
308# endif
309 & ad_t, ad_u, ad_v, &
310# if defined WEAK_CONSTRAINT && defined TIME_CONV
311 & ad_ubar, ad_vbar, &
312# endif
313# else
314 & ad_ubar, ad_vbar, &
315# endif
316 & ad_zeta)
317!***********************************************************************
318!
319! Imported variable declarations.
320!
321 integer, intent(in) :: ng, tile, model
322 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
323 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
324 integer, intent(in) :: Lold, Lnew
325 integer, intent(in) :: innLoop, outLoop
326 logical, intent(in) :: Ltrace
327!
328# ifdef ASSUMED_SHAPE
329# ifdef MASKING
330 real(r8), intent(in) :: rmask(LBi:,LBj:)
331 real(r8), intent(in) :: umask(LBi:,LBj:)
332 real(r8), intent(in) :: vmask(LBi:,LBj:)
333# endif
334# ifdef ADJUST_BOUNDARY
335# ifdef SOLVE3D
336 real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
337 real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
338 real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
339# endif
340 real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
341 real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
342 real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
343# endif
344# ifdef ADJUST_WSTRESS
345 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
346 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
347# endif
348# ifdef SOLVE3D
349# ifdef ADJUST_STFLUX
350 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
351# endif
352 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
353 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
354 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
355# if defined WEAK_CONSTRAINT && defined TIME_CONV
356 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
357 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
358# endif
359# else
360 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
361 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
362# endif
363 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
364# ifdef ADJUST_BOUNDARY
365# ifdef SOLVE3D
366 real(r8), intent(inout) :: d_t_obc(LBij:,:,:,:,:)
367 real(r8), intent(inout) :: d_u_obc(LBij:,:,:,:)
368 real(r8), intent(inout) :: d_v_obc(LBij:,:,:,:)
369# endif
370 real(r8), intent(inout) :: d_ubar_obc(LBij:,:,:)
371 real(r8), intent(inout) :: d_vbar_obc(LBij:,:,:)
372 real(r8), intent(inout) :: d_zeta_obc(LBij:,:,:)
373# endif
374# ifdef ADJUST_WSTRESS
375 real(r8), intent(inout) :: d_sustr(LBi:,LBj:,:)
376 real(r8), intent(inout) :: d_svstr(LBi:,LBj:,:)
377# endif
378# ifdef SOLVE3D
379# ifdef ADJUST_STFLUX
380 real(r8), intent(inout) :: d_stflx(LBi:,LBj:,:,:)
381# endif
382 real(r8), intent(inout) :: d_t(LBi:,LBj:,:,:)
383 real(r8), intent(inout) :: d_u(LBi:,LBj:,:)
384 real(r8), intent(inout) :: d_v(LBi:,LBj:,:)
385# if defined WEAK_CONSTRAINT && defined TIME_CONV
386 real(r8), intent(inout) :: d_ubar(LBi:,LBj:)
387 real(r8), intent(inout) :: d_vbar(LBi:,LBj:)
388# endif
389# else
390 real(r8), intent(inout) :: d_ubar(LBi:,LBj:)
391 real(r8), intent(inout) :: d_vbar(LBi:,LBj:)
392# endif
393 real(r8), intent(inout) :: d_zeta(LBi:,LBj:)
394# ifdef ADJUST_BOUNDARY
395# ifdef SOLVE3D
396 real(r8), intent(inout) :: nl_t_obc(LBij:,:,:,:,:,:)
397 real(r8), intent(inout) :: nl_u_obc(LBij:,:,:,:,:)
398 real(r8), intent(inout) :: nl_v_obc(LBij:,:,:,:,:)
399# endif
400 real(r8), intent(inout) :: nl_ubar_obc(LBij:,:,:,:)
401 real(r8), intent(inout) :: nl_vbar_obc(LBij:,:,:,:)
402 real(r8), intent(inout) :: nl_zeta_obc(LBij:,:,:,:)
403# endif
404# ifdef ADJUST_WSTRESS
405 real(r8), intent(inout) :: nl_ustr(LBi:,LBj:,:,:)
406 real(r8), intent(inout) :: nl_vstr(LBi:,LBj:,:,:)
407# endif
408# ifdef SOLVE3D
409# ifdef ADJUST_STFLUX
410 real(r8), intent(inout) :: nl_tflux(LBi:,LBj:,:,:,:)
411# endif
412 real(r8), intent(inout) :: nl_t(LBi:,LBj:,:,:,:)
413 real(r8), intent(inout) :: nl_u(LBi:,LBj:,:,:)
414 real(r8), intent(inout) :: nl_v(LBi:,LBj:,:,:)
415# if defined WEAK_CONSTRAINT && defined TIME_CONV
416 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
417 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
418# endif
419# else
420 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
421 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
422# endif
423 real(r8), intent(inout) :: nl_zeta(LBi:,LBj:,:)
424# ifdef ADJUST_BOUNDARY
425# ifdef SOLVE3D
426 real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
427 real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
428 real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
429# endif
430 real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
431 real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
432 real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
433# endif
434# ifdef ADJUST_WSTRESS
435 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
436 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
437# endif
438# ifdef SOLVE3D
439# ifdef ADJUST_STFLUX
440 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
441# endif
442 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
443 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
444 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
445# if defined WEAK_CONSTRAINT && defined TIME_CONV
446 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
447 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
448# endif
449# else
450 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
451 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
452# endif
453 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
454
455# else
456
457# ifdef MASKING
458 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
459 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
460 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
461# endif
462# ifdef ADJUST_BOUNDARY
463# ifdef SOLVE3D
464 real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, &
465 & Nbrec(ng),2,NT(ng))
466 real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
467 real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
468# endif
469 real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
470 real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
471 real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
472# endif
473# ifdef ADJUST_WSTRESS
474 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
475 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
476# endif
477# ifdef SOLVE3D
478# ifdef ADJUST_STFLUX
479 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, &
480 & Nfrec(ng),2,NT(ng))
481# endif
482 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
483 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
484 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
485# if defined WEAK_CONSTRAINT && defined TIME_CONV
486 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
487 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
488# endif
489# else
490 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
491 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
492# endif
493 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
494# ifdef ADJUST_BOUNDARY
495# ifdef SOLVE3D
496 real(r8), intent(inout) :: d_t_obc(LBij:UBij,N(ng),4, &
497 & Nbrec(ng),NT(ng))
498 real(r8), intent(inout) :: d_u_obc(LBij:UBij,N(ng),4,Nbrec(ng))
499 real(r8), intent(inout) :: d_v_obc(LBij:UBij,N(ng),4,Nbrec(ng))
500# endif
501 real(r8), intent(inout) :: d_ubar_obc(LBij:UBij,4,Nbrec(ng))
502 real(r8), intent(inout) :: d_vbar_obc(LBij:UBij,4,Nbrec(ng))
503 real(r8), intent(inout) :: d_zeta_obc(LBij:UBij,4,Nbrec(ng))
504# endif
505# ifdef ADJUST_WSTRESS
506 real(r8), intent(inout) :: d_sustr(LBi:UBi,LBj:UBj,Nfrec(ng))
507 real(r8), intent(inout) :: d_svstr(LBi:UBi,LBj:UBj,Nfrec(ng))
508# endif
509# ifdef SOLVE3D
510# ifdef ADJUST_STFLUX
511 real(r8), intent(inout) :: d_stflx(LBi:UBi,LBj:UBj, &
512 & Nfrec(ng),NT(ng))
513# endif
514 real(r8), intent(inout) :: d_t(LBi:UBi,LBj:UBj,N(ng),NT(ng))
515 real(r8), intent(inout) :: d_u(LBi:UBi,LBj:UBj,N(ng))
516 real(r8), intent(inout) :: d_v(LBi:UBi,LBj:UBj,N(ng))
517# if defined WEAK_CONSTRAINT && defined TIME_CONV
518 real(r8), intent(inout) :: d_ubar(LBi:UBi,LBj:UBj)
519 real(r8), intent(inout) :: d_vbar(LBi:UBi,LBj:UBj)
520# endif
521# else
522 real(r8), intent(inout) :: d_ubar(LBi:UBi,LBj:UBj)
523 real(r8), intent(inout) :: d_vbar(LBi:UBi,LBj:UBj)
524# endif
525 real(r8), intent(inout) :: d_zeta(LBi:UBi,LBj:UBj)
526# ifdef ADJUST_BOUNDARY
527# ifdef SOLVE3D
528 real(r8), intent(inout) :: nl_t_obc(LBij:UBij,N(ng),4, &
529 & Nbrec(ng),2,NT(ng))
530 real(r8), intent(inout) :: nl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
531 real(r8), intent(inout) :: nl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
532# endif
533 real(r8), intent(inout) :: nl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
534 real(r8), intent(inout) :: nl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
535 real(r8), intent(inout) :: nl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
536# endif
537# ifdef ADJUST_WSTRESS
538 real(r8), intent(inout) :: nl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
539 real(r8), intent(inout) :: nl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
540# endif
541# ifdef SOLVE3D
542# ifdef ADJUST_STFLUX
543 real(r8), intent(inout) :: nl_tflux(LBi:UBi,LBj:UBj, &
544 & Nfrec(ng),2,NT(ng))
545# endif
546 real(r8), intent(inout) :: nl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
547 real(r8), intent(inout) :: nl_u(LBi:UBi,LBj:UBj,N(ng),2)
548 real(r8), intent(inout) :: nl_v(LBi:UBi,LBj:UBj,N(ng),2)
549# if defined WEAK_CONSTRAINT && defined TIME_CONV
550 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,:)
551 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,:)
552# endif
553# else
554 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,:)
555 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,:)
556# endif
557 real(r8), intent(inout) :: nl_zeta(LBi:UBi,LBj:UBj,:)
558# ifdef ADJUST_BOUNDARY
559# ifdef SOLVE3D
560 real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, &
561 & Nbrec(ng),2,NT(ng))
562 real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
563 real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
564# endif
565 real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
566 real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
567 real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
568# endif
569# ifdef ADJUST_WSTRESS
570 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
571 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
572# endif
573# ifdef SOLVE3D
574# ifdef ADJUST_STFLUX
575 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, &
576 & Nfrec(ng),2,NT(ng))
577# endif
578 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
579 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
580 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
581# if defined WEAK_CONSTRAINT && defined TIME_CONV
582 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
583 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
584# endif
585# else
586 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
587 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
588# endif
589 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
590# endif
591!
592! Local variable declarations.
593!
594 logical :: Ltrans
595!
596 integer :: L1 = 1
597 integer :: L2 = 2
598
599 integer :: Linp, Lout, Lscale, Lwrk, Lwrk1, i, j, ic
600 integer :: info, itheta1
601!
602 real(r8) :: norm, zbeta
603
604 real(r8), dimension(2*NpostI-2) :: work
605# ifdef BEOFS_ONLY
606 integer :: LAMM
607 real(r8), dimension(0:NstateVar(ng)) :: dot
608# endif
609!
610 character (len=13) :: string
611
612 character (len=*), parameter :: MyFile = &
613 & __FILE__//", posterior_tile"
614
615# include "set_bounds.h"
616# ifdef BEOFS_ONLY
617!
618!-----------------------------------------------------------------------
619! Compute the action of the background error covariance matrix.
620!-----------------------------------------------------------------------
621!
622! NOTE: In the case of WEAK_CONSTRAINT and TIME_CONV, tl_ubar, tl_vbar
623! ad_ubar and ad_vbar are only passed as required but are not
624! used in subsequent calculations.
625!
626! Copy tl_var(Lnew) into ad_var(Lnew).
627!
628 IF (innloop.eq.0) THEN
629 lamm=lold
630 ELSE
631 lamm=lnew
632 END IF
633
634 CALL state_copy (ng, tile, &
635 & lbi, ubi, lbj, ubj, lbij, ubij, &
636 & lamm, lnew, &
637# ifdef ADJUST_BOUNDARY
638# ifdef SOLVE3D
639 & ad_t_obc, tl_t_obc, &
640 & ad_u_obc, tl_u_obc, &
641 & ad_v_obc, tl_v_obc, &
642# endif
643 & ad_ubar_obc, tl_ubar_obc, &
644 & ad_vbar_obc, tl_vbar_obc, &
645 & ad_zeta_obc, tl_zeta_obc, &
646# endif
647# ifdef ADJUST_WSTRESS
648 & ad_ustr, tl_ustr, &
649 & ad_vstr, tl_vstr, &
650# endif
651# ifdef SOLVE3D
652# ifdef ADJUST_STFLUX
653 & ad_tflux, tl_tflux, &
654# endif
655 & ad_t, tl_t, &
656 & ad_u, tl_u, &
657 & ad_v, tl_v, &
658# if defined WEAK_CONSTRAINT && defined TIME_CONV
659 & ad_ubar, tl_ubar, &
660 & ad_vbar, tl_vbar, &
661# endif
662# else
663 & ad_ubar, tl_ubar, &
664 & ad_vbar, tl_vbar, &
665# endif
666 & ad_zeta, tl_zeta)
667!
668 lwrk=1
669 IF ((innloop.gt.0).or.ltrace) THEN
670 linp=1
671 lout=2
672!
673!-----------------------------------------------------------------------
674! Compute norm Delta(k) as the dot-product between the new vector
675! and previous Lanczos vector.
676!-----------------------------------------------------------------------
677!
678! Compute current iteration norm Delta(k) used to compute tri-diagonal
679! matrix T(k) in the Lanczos recurrence.
680!
681 CALL state_dotprod (ng, tile, model, &
682 & lbi, ubi, lbj, ubj, lbij, ubij, &
683 & nstatevar(ng), dot(0:), &
684# ifdef MASKING
685 & rmask, umask, vmask, &
686# endif
687# ifdef ADJUST_BOUNDARY
688# ifdef SOLVE3D
689 & tl_t_obc(:,:,:,:,lold,:), &
690 & tl_t_obc(:,:,:,:,lnew,:), &
691 & tl_u_obc(:,:,:,:,lold), &
692 & tl_u_obc(:,:,:,:,lnew), &
693 & tl_v_obc(:,:,:,:,lold), &
694 & tl_v_obc(:,:,:,:,lnew), &
695# endif
696 & tl_ubar_obc(:,:,:,lold), &
697 & tl_ubar_obc(:,:,:,lnew), &
698 & tl_vbar_obc(:,:,:,lold), &
699 & tl_vbar_obc(:,:,:,lnew), &
700 & tl_zeta_obc(:,:,:,lold), &
701 & tl_zeta_obc(:,:,:,lnew), &
702# endif
703# ifdef ADJUST_WSTRESS
704 & tl_ustr(:,:,:,lold), tl_ustr(:,:,:,lnew), &
705 & tl_vstr(:,:,:,lold), tl_vstr(:,:,:,lnew), &
706# endif
707# ifdef SOLVE3D
708# ifdef ADJUST_STFLUX
709 & tl_tflux(:,:,:,lold,:), &
710 & tl_tflux(:,:,:,lnew,:), &
711# endif
712 & tl_t(:,:,:,lold,:), tl_t(:,:,:,lnew,:), &
713 & tl_u(:,:,:,lold), tl_u(:,:,:,lnew), &
714 & tl_v(:,:,:,lold), tl_v(:,:,:,lnew), &
715# else
716 & tl_ubar(:,:,lold), tl_ubar(:,:,lnew), &
717 & tl_vbar(:,:,lold), tl_vbar(:,:,lnew), &
718# endif
719 & tl_zeta(:,:,lold), tl_zeta(:,:,lnew))
720
721 ae_delta(innloop,outloop)=dot(0)
722# else
723!
724!-----------------------------------------------------------------------
725! Compute the action of the act analysis error covariance matrix.
726!-----------------------------------------------------------------------
727!
728! NOTE: In the case of weak constraint ("WEAK_CONSTRAINT") and
729! and time convolutions ("TIME_CONV"), the state arrays
730! tl_ubar, tl_vbar, ad_ubar, and ad_vbar are only passed
731! as required by the "state" operators routines but they
732! are not used in subsequent calculations.
733!
734! Copy tl_var(Lold) into ad_var(Lnew). See NOTE above.
735!
736 CALL state_copy (ng, tile, &
737 & lbi, ubi, lbj, ubj, lbij, ubij, &
738 & lold, lnew, &
739# ifdef ADJUST_BOUNDARY
740# ifdef SOLVE3D
741 & ad_t_obc, tl_t_obc, &
742 & ad_u_obc, tl_u_obc, &
743 & ad_v_obc, tl_v_obc, &
744# endif
745 & ad_ubar_obc, tl_ubar_obc, &
746 & ad_vbar_obc, tl_vbar_obc, &
747 & ad_zeta_obc, tl_zeta_obc, &
748# endif
749# ifdef ADJUST_WSTRESS
750 & ad_ustr, tl_ustr, &
751 & ad_vstr, tl_vstr, &
752# endif
753# ifdef SOLVE3D
754# ifdef ADJUST_STFLUX
755 & ad_tflux, tl_tflux, &
756# endif
757 & ad_t, tl_t, &
758 & ad_u, tl_u, &
759 & ad_v, tl_v, &
760# if defined WEAK_CONSTRAINT && defined TIME_CONV
761 & ad_ubar, tl_ubar, &
762 & ad_vbar, tl_vbar, &
763# endif
764# else
765 & ad_ubar, tl_ubar, &
766 & ad_vbar, tl_vbar, &
767# endif
768 & ad_zeta, tl_zeta)
769!
770 lwrk=1
771 IF ((innloop.gt.0).or.ltrace) THEN
772 linp=1
773 lout=2
774 CALL analysis_error (ng, tile, model, &
775 & lbi, ubi, lbj, ubj, lbij, ubij, &
776 & imins, imaxs, jmins, jmaxs, &
777 & linp, lout, lwrk, &
778 & innloop, outloop, &
779# ifdef MASKING
780 & rmask, umask, vmask, &
781# endif
782# ifdef ADJUST_BOUNDARY
783# ifdef SOLVE3D
784 & ad_t_obc, ad_u_obc, ad_v_obc, &
785# endif
786 & ad_ubar_obc, ad_vbar_obc, &
787 & ad_zeta_obc, &
788# endif
789# ifdef ADJUST_WSTRESS
790 & ad_ustr, ad_vstr, &
791# endif
792# ifdef SOLVE3D
793# ifdef ADJUST_STFLUX
794 & ad_tflux, &
795# endif
796 & ad_t, ad_u, ad_v, &
797# else
798 & ad_ubar, ad_vbar, &
799# endif
800 & ad_zeta, &
801# ifdef ADJUST_BOUNDARY
802# ifdef SOLVE3D
803 & tl_t_obc, tl_u_obc, tl_v_obc, &
804# endif
805 & tl_ubar_obc, tl_vbar_obc, &
806 & tl_zeta_obc, &
807# endif
808# ifdef ADJUST_WSTRESS
809 & tl_ustr, tl_vstr, &
810# endif
811# ifdef SOLVE3D
812# ifdef ADJUST_STFLUX
813 & tl_tflux, &
814# endif
815 & tl_t, tl_u, tl_v, &
816# else
817 & tl_ubar, tl_vbar, &
818# endif
819 & tl_zeta, &
820# ifdef ADJUST_BOUNDARY
821# ifdef SOLVE3D
822 & nl_t_obc, nl_u_obc, nl_v_obc, &
823# endif
824 & nl_ubar_obc, nl_vbar_obc, &
825 & nl_zeta_obc, &
826# endif
827# ifdef ADJUST_WSTRESS
828 & nl_ustr, nl_vstr, &
829# endif
830# ifdef SOLVE3D
831# ifdef ADJUST_STFLUX
832 & nl_tflux, &
833# endif
834 & nl_t, nl_u, nl_v, &
835# else
836 & nl_ubar, nl_vbar, &
837# endif
838 & nl_zeta)
839# endif
840 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
841!
842! Check that the analysis error covariance matrix is positive
843! definite, or report the trace estimate.
844!
845 IF (ltrace) THEN
846 ae_trace(innloop+1)=ae_delta(innloop,outloop)
847!
848 SELECT CASE (hss(ng)%IOtype)
849 CASE (io_nf90)
850 CALL netcdf_put_fvar (ng, model, hss(ng)%name, &
851 & 'ae_trace', ae_trace(innloop+1:), &
852 & (/innloop+1/), (/1/), &
853 & ncid = hss(ng)%ncid)
854
855# if defined PIO_LIB && defined DISTRIBUTE
856 CASE (io_pio)
857 CALL pio_netcdf_put_fvar (ng, model, hss(ng)%name, &
858 & 'ae_trace', &
859 & ae_trace(innloop+1:), &
860 & (/innloop+1/), (/1/), &
861 & piofile = hss(ng)%pioFile)
862# endif
863 END SELECT
864 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
865!
866 IF (master) THEN
867 WRITE (stdout,10) outloop, innloop, &
868 & ae_trace(innloop+1)
869 10 FORMAT (1x,'(',i3.3,',',i3.3,'): ', &
870 & 'Analysis Error Trace Estimate, ae_trace = ', &
871 & 1p,e14.7)
872 END IF
873 RETURN
874 END IF
875 IF (ae_delta(innloop,outloop).le.0.0_r8) THEN
876 WRITE (stdout,*) ' AE_DELTA not positive.'
877 WRITE (stdout,*) ' AE_DELTA = ', ae_delta(innloop,outloop), &
878 & ', outer = ', outloop, ', inner = ', innloop
879 exit_flag=8
880 RETURN
881 END IF
882 END IF
883!
884! Apply the Lanczos recurrence and orthonormalize.
885!
886 linp=1
887 lout=2
888 CALL lanczos (ng, tile, model, &
889 & lbi, ubi, lbj, ubj, lbij, ubij, &
890 & imins, imaxs, jmins, jmaxs, &
891 & linp, lout, lwrk, &
892 & innloop, outloop, &
893# ifdef MASKING
894 & rmask, umask, vmask, &
895# endif
896# ifdef ADJUST_BOUNDARY
897# ifdef SOLVE3D
898 & tl_t_obc, tl_u_obc, tl_v_obc, &
899# endif
900 & tl_ubar_obc, tl_vbar_obc, &
901 & tl_zeta_obc, &
902# endif
903# ifdef ADJUST_WSTRESS
904 & tl_ustr, tl_vstr, &
905# endif
906# ifdef SOLVE3D
907# ifdef ADJUST_STFLUX
908 & tl_tflux, &
909# endif
910 & tl_t, tl_u, tl_v, &
911# else
912 & tl_ubar, tl_vbar, &
913# endif
914 & tl_zeta, &
915# ifdef ADJUST_BOUNDARY
916# ifdef SOLVE3D
917 & ad_t_obc, ad_u_obc, ad_v_obc, &
918# endif
919 & ad_ubar_obc, ad_vbar_obc, &
920 & ad_zeta_obc, &
921# endif
922# ifdef ADJUST_WSTRESS
923 & ad_ustr, ad_vstr, &
924# endif
925# ifdef SOLVE3D
926# ifdef ADJUST_STFLUX
927 & ad_tflux, &
928# endif
929 & ad_t, ad_u, ad_v, &
930# else
931 & ad_ubar, ad_vbar, &
932# endif
933 & ad_zeta)
934 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
935!
936! Compute new direction, d(k+1).
937!
938 linp=1
939 lout=2
940 CALL new_direction (ng, tile, model, &
941 & lbi, ubi, lbj, ubj, lbij, ubij, &
942 & imins, imaxs, jmins, jmaxs, &
943 & linp, lout, &
944# ifdef MASKING
945 & rmask, umask, vmask, &
946# endif
947# ifdef ADJUST_BOUNDARY
948# ifdef SOLVE3D
949 & ad_t_obc, ad_u_obc, ad_v_obc, &
950# endif
951 & ad_ubar_obc, ad_vbar_obc, &
952 & ad_zeta_obc, &
953# endif
954# ifdef ADJUST_WSTRESS
955 & ad_ustr, ad_vstr, &
956# endif
957# ifdef SOLVE3D
958# ifdef ADJUST_STFLUX
959 & ad_tflux, &
960# endif
961 & ad_t, ad_u, ad_v, &
962# else
963 & ad_ubar, ad_vbar, &
964# endif
965 & ad_zeta, &
966# ifdef ADJUST_BOUNDARY
967# ifdef SOLVE3D
968 & d_t_obc, d_u_obc, d_v_obc, &
969# endif
970 & d_ubar_obc, d_vbar_obc, &
971 & d_zeta_obc, &
972# endif
973# ifdef ADJUST_WSTRESS
974 & d_sustr, d_svstr, &
975# endif
976# ifdef SOLVE3D
977# ifdef ADJUST_STFLUX
978 & d_stflx, &
979# endif
980 & d_t, d_u, d_v, &
981# else
982 & d_ubar, d_vbar, &
983# endif
984 & d_zeta)
985 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
986!
987!-----------------------------------------------------------------------
988! Determine the eigenvalues and eigenvectors of the tridiagonal matrix.
989! These will be used on the last inner-loop to compute the eigenvectors
990! of the Hessian.
991!-----------------------------------------------------------------------
992!
993 IF (innloop.gt.0) THEN
994 DO i=1,innloop
995 ae_ritz(i,outloop)=ae_delta(i,outloop)
996 END DO
997 DO i=1,innloop-1
998 ae_tmatrix(i,1)=ae_beta(i+1,outloop)
999 END DO
1000!
1001! Use the LAPACK routine DSTEQR to compute the eigenvectors and
1002! eigenvalues of the tridiagonal matrix. If applicable, the
1003! eigenpairs is computed by master thread only. Notice that on
1004! exit, the matrix "ae_Tmatrix" is destroyed.
1005!
1006 IF (master) THEN
1007 CALL dsteqr ('I', innloop, ae_ritz(1,outloop), ae_tmatrix, &
1008 & ae_zv, nposti, work, info)
1009 END IF
1010# ifdef DISTRIBUTE
1011 CALL mp_bcasti (ng, model, info)
1012# endif
1013 IF (info.ne.0) THEN
1014 WRITE (stdout,*) ' Error in DSTEQR: info = ', info
1015 exit_flag=8
1016 RETURN
1017 END IF
1018# ifdef DISTRIBUTE
1019 CALL mp_bcastf (ng, model, ae_ritz(:,outloop))
1020 CALL mp_bcastf (ng, model, ae_zv)
1021# endif
1022!
1023! Estimate the Ritz value error bounds.
1024!
1025 DO i=1,innloop
1026 ae_ritzerr(i,outloop)=abs(ae_beta(innloop+1,outloop)* &
1027 & ae_zv(innloop,i))
1028 END DO
1029!
1030! Check for exploding or negative Ritz values.
1031!
1032 DO i=1,innloop
1033 IF (ae_ritz(i,outloop).lt.0.0_r8) THEN
1034 WRITE (stdout,*) ' Negative Ritz value found.'
1035 exit_flag=8
1036 RETURN
1037 END IF
1038 END DO
1039!
1040! Calculate the converged eigenvectors of the Hessian.
1041!
1042 IF (innloop.eq.nposti) THEN
1043 ritzmaxerr=hevecerr
1044 DO i=1,innloop
1045 ae_ritzerr(i,outloop)=ae_ritzerr(i,outloop)/ &
1046 & ae_ritz(nposti,outloop)
1047 END DO
1048 lwrk=2
1049 linp=1
1050 lout=2
1051 CALL posterior_eofs (ng, tile, model, &
1052 & lbi, ubi, lbj, ubj, lbij, ubij, &
1053 & imins, imaxs, jmins, jmaxs, &
1054 & linp, lout, lwrk, &
1055 & innloop, outloop, &
1056# ifdef MASKING
1057 & rmask, umask, vmask, &
1058# endif
1059# ifdef ADJUST_BOUNDARY
1060# ifdef SOLVE3D
1061 & nl_t_obc, nl_u_obc, nl_v_obc, &
1062# endif
1063 & nl_ubar_obc, nl_vbar_obc, &
1064 & nl_zeta_obc, &
1065# endif
1066# ifdef ADJUST_WSTRESS
1067 & nl_ustr, nl_vstr, &
1068# endif
1069# ifdef SOLVE3D
1070# ifdef ADJUST_STFLUX
1071 & nl_tflux, &
1072# endif
1073 & nl_t, nl_u, nl_v, &
1074# else
1075 & nl_ubar, nl_vbar, &
1076# endif
1077 & nl_zeta, &
1078# ifdef ADJUST_BOUNDARY
1079# ifdef SOLVE3D
1080 & tl_t_obc, tl_u_obc, tl_v_obc, &
1081# endif
1082 & tl_ubar_obc, tl_vbar_obc, &
1083 & tl_zeta_obc, &
1084# endif
1085# ifdef ADJUST_WSTRESS
1086 & tl_ustr, tl_vstr, &
1087# endif
1088# ifdef SOLVE3D
1089# ifdef ADJUST_STFLUX
1090 & tl_tflux, &
1091# endif
1092 & tl_t, tl_u, tl_v, &
1093# else
1094 & tl_ubar, tl_vbar, &
1095# endif
1096 & tl_zeta, &
1097# ifdef ADJUST_BOUNDARY
1098# ifdef SOLVE3D
1099 & ad_t_obc, ad_u_obc, ad_v_obc, &
1100# endif
1101 & ad_ubar_obc, ad_vbar_obc, &
1102 & ad_zeta_obc, &
1103# endif
1104# ifdef ADJUST_WSTRESS
1105 & ad_ustr, ad_vstr, &
1106# endif
1107# ifdef SOLVE3D
1108# ifdef ADJUST_STFLUX
1109 & ad_tflux, &
1110# endif
1111 & ad_t, ad_u, ad_v, &
1112# else
1113 & ad_ubar, ad_vbar, &
1114# endif
1115 & ad_zeta)
1116 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1117 IF (master.and.(nconvritz.eq.0)) THEN
1118 WRITE(stdout,*) ' No converged Hesssian eigenvectors found.'
1119 END IF
1120 END IF
1121 END IF
1122!
1123!-----------------------------------------------------------------------
1124! Set TLM initial conditions for next inner loop, X(k+1).
1125!-----------------------------------------------------------------------
1126!
1127! X(k+1) = tau(k+1) * d(k+1)
1128!
1129 linp=1
1130 lout=2
1131 CALL tl_new_vector (ng, tile, model, &
1132 & lbi, ubi, lbj, ubj, lbij, ubij, &
1133 & imins, imaxs, jmins, jmaxs, &
1134 & linp, lout, &
1135 & innloop, outloop, &
1136# ifdef MASKING
1137 & rmask, umask, vmask, &
1138# endif
1139# ifdef ADJUST_BOUNDARY
1140# ifdef SOLVE3D
1141 & d_t_obc, d_u_obc, d_v_obc, &
1142# endif
1143 & d_ubar_obc, d_vbar_obc, &
1144 & d_zeta_obc, &
1145# endif
1146# ifdef ADJUST_WSTRESS
1147 & d_sustr, d_svstr, &
1148# endif
1149# ifdef SOLVE3D
1150# ifdef ADJUST_STFLUX
1151 & d_stflx, &
1152# endif
1153 & d_t, d_u, d_v, &
1154# else
1155 & d_ubar, d_vbar, &
1156# endif
1157 & d_zeta, &
1158# ifdef ADJUST_BOUNDARY
1159# ifdef SOLVE3D
1160 & tl_t_obc, tl_u_obc, tl_v_obc, &
1161# endif
1162 & tl_ubar_obc, tl_vbar_obc, &
1163 & tl_zeta_obc, &
1164# endif
1165# ifdef ADJUST_WSTRESS
1166 & tl_ustr, tl_vstr, &
1167# endif
1168# ifdef SOLVE3D
1169# ifdef ADJUST_STFLUX
1170 & tl_tflux, &
1171# endif
1172 & tl_t, tl_u, tl_v, &
1173# else
1174 & tl_ubar, tl_vbar, &
1175# endif
1176 & tl_zeta, &
1177# ifdef ADJUST_BOUNDARY
1178# ifdef SOLVE3D
1179 & ad_t_obc, ad_u_obc, ad_v_obc, &
1180# endif
1181 & ad_ubar_obc, ad_vbar_obc, &
1182 & ad_zeta_obc, &
1183# endif
1184# ifdef ADJUST_WSTRESS
1185 & ad_ustr, ad_vstr, &
1186# endif
1187# ifdef SOLVE3D
1188# ifdef ADJUST_STFLUX
1189 & ad_tflux, &
1190# endif
1191 & ad_t, ad_u, ad_v, &
1192# else
1193 & ad_ubar, ad_vbar, &
1194# endif
1195 & ad_zeta)
1196!
1197!-----------------------------------------------------------------------
1198! Report posterior analysis error covariance estimation parameters.
1199!-----------------------------------------------------------------------
1200!
1201 IF (master) THEN
1202 IF (inner.eq.0) THEN
1203 WRITE (stdout,20) outloop, innloop, &
1204 & ae_gnorm(outloop)
1205 20 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1206 & 'Analysis Error gradient norm, ae_Gnorm = ', &
1207 & 1p,e14.7,/)
1208 END IF
1209 IF (innloop.gt.0) THEN
1210 WRITE (stdout,30) ritzmaxerr
1211 30 FORMAT (/,' Ritz Eigenvalues and relative accuracy: ', &
1212 & 'RitzMaxErr = ',1p,e14.7,/)
1213 ic=0
1214 DO i=1,innloop
1215 IF (ae_ritzerr(i,outloop).le.ritzmaxerr) THEN
1216 string='converged'
1217 ic=ic+1
1218 WRITE (stdout,40) i, ae_ritz(i,outloop), &
1219 & ae_ritzerr(i,outloop), &
1220 & trim(adjustl(string)), ic
1221 40 FORMAT(5x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a,2x, &
1222 & '(Good='i3.3,')')
1223 ELSE
1224 string='not converged'
1225 WRITE (stdout,50) i, ae_ritz(i,outloop), &
1226 & ae_ritzerr(i,outloop), &
1227 & trim(adjustl(string))
1228 50 FORMAT(5x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a)
1229 END IF
1230 END DO
1231 WRITE (stdout,'(/)')
1232 END IF
1233 END IF
1234
1235 RETURN
1236 END SUBROUTINE posterior_tile
1237!
1238!***********************************************************************
1239 SUBROUTINE tl_new_vector (ng, tile, model, &
1240 & LBi, UBi, LBj, UBj, LBij, UBij, &
1241 & IminS, ImaxS, JminS, JmaxS, &
1242 & Linp, Lout, &
1243 & innLoop, outLoop, &
1244# ifdef MASKING
1245 & rmask, umask, vmask, &
1246# endif
1247# ifdef ADJUST_BOUNDARY
1248# ifdef SOLVE3D
1249 & d_t_obc, d_u_obc, d_v_obc, &
1250# endif
1251 & d_ubar_obc, d_vbar_obc, &
1252 & d_zeta_obc, &
1253# endif
1254# ifdef ADJUST_WSTRESS
1255 & d_sustr, d_svstr, &
1256# endif
1257# ifdef SOLVE3D
1258# ifdef ADJUST_STFLUX
1259 & d_stflx, &
1260# endif
1261 & d_t, d_u, d_v, &
1262# else
1263 & d_ubar, d_vbar, &
1264# endif
1265 & d_zeta, &
1266# ifdef ADJUST_BOUNDARY
1267# ifdef SOLVE3D
1268 & tl_t_obc, tl_u_obc, tl_v_obc, &
1269# endif
1270 & tl_ubar_obc, tl_vbar_obc, &
1271 & tl_zeta_obc, &
1272# endif
1273# ifdef ADJUST_WSTRESS
1274 & tl_ustr, tl_vstr, &
1275# endif
1276# ifdef SOLVE3D
1277# ifdef ADJUST_STFLUX
1278 & tl_tflux, &
1279# endif
1280 & tl_t, tl_u, tl_v, &
1281# else
1282 & tl_ubar, tl_vbar, &
1283# endif
1284 & tl_zeta, &
1285# ifdef ADJUST_BOUNDARY
1286# ifdef SOLVE3D
1287 & ad_t_obc, ad_u_obc, ad_v_obc, &
1288# endif
1289 & ad_ubar_obc, ad_vbar_obc, &
1290 & ad_zeta_obc, &
1291# endif
1292# ifdef ADJUST_WSTRESS
1293 & ad_ustr, ad_vstr, &
1294# endif
1295# ifdef SOLVE3D
1296# ifdef ADJUST_STFLUX
1297 & ad_tflux, &
1298# endif
1299 & ad_t, ad_u, ad_v, &
1300# else
1301 & ad_ubar, ad_vbar, &
1302# endif
1303 & ad_zeta)
1304!***********************************************************************
1305!
1306! Imported variable declarations.
1307!
1308 integer, intent(in) :: ng, tile, model
1309 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
1310 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
1311 integer, intent(in) :: Linp, Lout
1312 integer, intent(in) :: innLoop, outLoop
1313!
1314# ifdef ASSUMED_SHAPE
1315# ifdef MASKING
1316 real(r8), intent(in) :: rmask(LBi:,LBj:)
1317 real(r8), intent(in) :: umask(LBi:,LBj:)
1318 real(r8), intent(in) :: vmask(LBi:,LBj:)
1319# endif
1320# ifdef ADJUST_BOUNDARY
1321# ifdef SOLVE3D
1322 real(r8), intent(inout) :: d_t_obc(LBij:,:,:,:,:)
1323 real(r8), intent(inout) :: d_u_obc(LBij:,:,:,:)
1324 real(r8), intent(inout) :: d_v_obc(LBij:,:,:,:)
1325# endif
1326 real(r8), intent(inout) :: d_ubar_obc(LBij:,:,:)
1327 real(r8), intent(inout) :: d_vbar_obc(LBij:,:,:)
1328 real(r8), intent(inout) :: d_zeta_obc(LBij:,:,:)
1329# endif
1330# ifdef ADJUST_WSTRESS
1331 real(r8), intent(in) :: d_sustr(LBi:,LBj:,:)
1332 real(r8), intent(in) :: d_svstr(LBi:,LBj:,:)
1333# endif
1334# ifdef SOLVE3D
1335# ifdef ADJUST_STFLUX
1336 real(r8), intent(in) :: d_stflx(LBi:,LBj:,:,:)
1337# endif
1338 real(r8), intent(in) :: d_t(LBi:,LBj:,:,:)
1339 real(r8), intent(in) :: d_u(LBi:,LBj:,:)
1340 real(r8), intent(in) :: d_v(LBi:,LBj:,:)
1341# else
1342 real(r8), intent(in) :: d_ubar(LBi:,LBj:)
1343 real(r8), intent(in) :: d_vbar(LBi:,LBj:)
1344# endif
1345 real(r8), intent(in) :: d_zeta(LBi:,LBj:)
1346# ifdef ADJUST_BOUNDARY
1347# ifdef SOLVE3D
1348 real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
1349 real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
1350 real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
1351# endif
1352 real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
1353 real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
1354 real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
1355# endif
1356# ifdef ADJUST_WSTRESS
1357 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
1358 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
1359# endif
1360# ifdef SOLVE3D
1361# ifdef ADJUST_STFLUX
1362 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
1363# endif
1364 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
1365 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
1366 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
1367# else
1368 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
1369 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
1370# endif
1371 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
1372# ifdef ADJUST_BOUNDARY
1373# ifdef SOLVE3D
1374 real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
1375 real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
1376 real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
1377# endif
1378 real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
1379 real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
1380 real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
1381# endif
1382# ifdef ADJUST_WSTRESS
1383 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
1384 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
1385# endif
1386# ifdef SOLVE3D
1387# ifdef ADJUST_STFLUX
1388 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
1389# endif
1390 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
1391 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
1392 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
1393# else
1394 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
1395 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
1396# endif
1397 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
1398
1399# else
1400
1401# ifdef MASKING
1402 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1403 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1404 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1405# endif
1406# ifdef ADJUST_BOUNDARY
1407# ifdef SOLVE3D
1408 real(r8), intent(in) :: d_t_obc(LBij:UBij,N(ng),4, &
1409 & Nbrec(ng),NT(ng))
1410 real(r8), intent(in) :: d_u_obc(LBij:UBij,N(ng),4,Nbrec(ng))
1411 real(r8), intent(in) :: d_v_obc(LBij:UBij,N(ng),4,Nbrec(ng))
1412# endif
1413 real(r8), intent(in) :: d_ubar_obc(LBij:UBij,4,Nbrec(ng))
1414 real(r8), intent(in) :: d_vbar_obc(LBij:UBij,4,Nbrec(ng))
1415 real(r8), intent(in) :: d_zeta_obc(LBij:UBij,4,Nbrec(ng))
1416# endif
1417# ifdef ADJUST_WSTRESS
1418 real(r8), intent(in) :: d_sustr(LBi:UBi,LBj:UBj,Nfrec(ng))
1419 real(r8), intent(in) :: d_svstr(LBi:UBi,LBj:UBj,Nfrec(ng))
1420# endif
1421# ifdef SOLVE3D
1422# ifdef ADJUST_STFLUX
1423 real(r8), intent(in) :: d_stflx(LBi:UBi,LBj:UBj, &
1424 & Nfrec(ng),NT(ng))
1425# endif
1426 real(r8), intent(in) :: d_t(LBi:UBi,LBj:UBj,N(ng),NT(ng))
1427 real(r8), intent(in) :: d_u(LBi:UBi,LBj:UBj,N(ng))
1428 real(r8), intent(in) :: d_v(LBi:UBi,LBj:UBj,N(ng))
1429# else
1430 real(r8), intent(in) :: d_ubar(LBi:UBi,LBj:UBj)
1431 real(r8), intent(in) :: d_vbar(LBi:UBi,LBj:UBj)
1432# endif
1433 real(r8), intent(in) :: d_zeta(LBi:UBi,LBj:UBj)
1434# ifdef ADJUST_BOUNDARY
1435# ifdef SOLVE3D
1436 real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, &
1437 & Nbrec(ng),2,NT(ng))
1438 real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
1439 real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
1440# endif
1441 real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
1442 real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
1443 real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
1444# endif
1445# ifdef ADJUST_WSTRESS
1446 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
1447 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
1448# endif
1449# ifdef SOLVE3D
1450# ifdef ADJUST_STFLUX
1451 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, &
1452 & Nfrec(ng),2,NT(ng))
1453# endif
1454 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1455 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
1456 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
1457# else
1458 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
1459 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
1460# endif
1461 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
1462# ifdef ADJUST_BOUNDARY
1463# ifdef SOLVE3D
1464 real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, &
1465 & Nbrec(ng),2,NT(ng))
1466 real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
1467 real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
1468# endif
1469 real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
1470 real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
1471 real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
1472# endif
1473# ifdef ADJUST_WSTRESS
1474 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
1475 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
1476# endif
1477# ifdef SOLVE3D
1478# ifdef ADJUST_STFLUX
1479 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, &
1480 & Nfrec(ng),2,NT(ng))
1481# endif
1482 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1483 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
1484 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
1485# else
1486 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
1487 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
1488# endif
1489 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
1490# endif
1491!
1492! Local variable declarations.
1493!
1494 integer :: i, j, k, lstr, rec
1495 integer :: ib, ir, it
1496
1497 real(r8) :: fac, fac1, fac2
1498
1499# include "set_bounds.h"
1500!
1501!-----------------------------------------------------------------------
1502! Compute new starting tangent linear state vector, X(k+1).
1503!-----------------------------------------------------------------------
1504!
1505! Free-surface.
1506!
1507 DO j=jstrt,jendt
1508 DO i=istrt,iendt
1509 tl_zeta(i,j,lout)=d_zeta(i,j)
1510# ifdef MASKING
1511 tl_zeta(i,j,lout)=tl_zeta(i,j,lout)*rmask(i,j)
1512# endif
1513 END DO
1514 END DO
1515
1516# ifdef ADJUST_BOUNDARY
1517!
1518! Free-surface open boundaries.
1519!
1520 IF (any(lobc(:,isfsur,ng))) THEN
1521 DO ir=1,nbrec(ng)
1522 IF ((lobc(iwest,isfsur,ng)).and. &
1523 & domain(ng)%Western_Edge(tile)) THEN
1524 ib=iwest
1525 DO j=jstr,jend
1526 tl_zeta_obc(j,ib,ir,lout)=d_zeta_obc(j,ib,ir)
1527# ifdef MASKING
1528 tl_zeta_obc(j,ib,ir,lout)=tl_zeta_obc(j,ib,ir,lout)* &
1529 & rmask(istr-1,j)
1530# endif
1531 END DO
1532 END IF
1533 IF ((lobc(ieast,isfsur,ng)).and. &
1534 & domain(ng)%Eastern_Edge(tile)) THEN
1535 ib=ieast
1536 DO j=jstr,jend
1537 tl_zeta_obc(j,ib,ir,lout)=d_zeta_obc(j,ib,ir)
1538# ifdef MASKING
1539 tl_zeta_obc(j,ib,ir,lout)=tl_zeta_obc(j,ib,ir,lout)* &
1540 & rmask(iend+1,j)
1541# endif
1542 END DO
1543 END IF
1544 IF ((lobc(isouth,isfsur,ng)).and. &
1545 & domain(ng)%Southern_Edge(tile)) THEN
1546 ib=isouth
1547 DO i=istr,iend
1548 tl_zeta_obc(i,ib,ir,lout)=d_zeta_obc(i,ib,ir)
1549# ifdef MASKING
1550 tl_zeta_obc(i,ib,ir,lout)=tl_zeta_obc(i,ib,ir,lout)* &
1551 & rmask(i,jstr-1)
1552# endif
1553 END DO
1554 END IF
1555 IF ((lobc(inorth,isfsur,ng)).and. &
1556 & domain(ng)%Northern_Edge(tile)) THEN
1557 ib=inorth
1558 DO i=istr,iend
1559 tl_zeta_obc(i,ib,ir,lout)=d_zeta_obc(i,ib,ir)
1560# ifdef MASKING
1561 tl_zeta_obc(i,ib,ir,lout)=tl_zeta_obc(i,ib,ir,lout)* &
1562 & rmask(i,jend+1)
1563# endif
1564 END DO
1565 END IF
1566 END DO
1567 END IF
1568# endif
1569
1570# ifndef SOLVE3D
1571!
1572! 2D U-momentum.
1573!
1574 DO j=jstrt,jendt
1575 DO i=istrp,iendt
1576 tl_ubar(i,j,lout)=d_ubar(i,j)
1577# ifdef MASKING
1578 tl_ubar(i,j,lout)=tl_ubar(i,j,lout)*umask(i,j)
1579# endif
1580 END DO
1581 END DO
1582# endif
1583
1584# ifdef ADJUST_BOUNDARY
1585!
1586! 2D U-momentum open boundaries.
1587!
1588 IF (any(lobc(:,isubar,ng))) THEN
1589 DO ir=1,nbrec(ng)
1590 IF ((lobc(iwest,isubar,ng)).and. &
1591 & domain(ng)%Western_Edge(tile)) THEN
1592 ib=iwest
1593 DO j=jstr,jend
1594 tl_ubar_obc(j,ib,ir,lout)=d_ubar_obc(j,ib,ir)
1595# ifdef MASKING
1596 tl_ubar_obc(j,ib,ir,lout)=tl_ubar_obc(j,ib,ir,lout)* &
1597 & umask(istr,j)
1598# endif
1599 END DO
1600 END IF
1601 IF ((lobc(ieast,isubar,ng)).and. &
1602 & domain(ng)%Eastern_Edge(tile)) THEN
1603 ib=ieast
1604 DO j=jstr,jend
1605 tl_ubar_obc(j,ib,ir,lout)=d_ubar_obc(j,ib,ir)
1606# ifdef MASKING
1607 tl_ubar_obc(j,ib,ir,lout)=tl_ubar_obc(j,ib,ir,lout)* &
1608 & umask(iend+1,j)
1609# endif
1610 END DO
1611 END IF
1612 IF ((lobc(isouth,isubar,ng)).and. &
1613 & domain(ng)%Southern_Edge(tile)) THEN
1614 ib=isouth
1615 DO i=istru,iend
1616 tl_ubar_obc(i,ib,ir,lout)=d_ubar_obc(i,ib,ir)
1617# ifdef MASKING
1618 tl_ubar_obc(i,ib,ir,lout)=tl_ubar_obc(i,ib,ir,lout)* &
1619 & umask(i,jstr-1)
1620# endif
1621 END DO
1622 END IF
1623 IF ((lobc(inorth,isubar,ng)).and. &
1624 & domain(ng)%Northern_Edge(tile)) THEN
1625 ib=inorth
1626 DO i=istru,iend
1627 tl_ubar_obc(i,ib,ir,lout)=d_ubar_obc(i,ib,ir)
1628# ifdef MASKING
1629 tl_ubar_obc(i,ib,ir,lout)=tl_ubar_obc(i,ib,ir,lout)* &
1630 & umask(i,jend+1)
1631# endif
1632 END DO
1633 END IF
1634 END DO
1635 END IF
1636# endif
1637
1638# ifndef SOLVE3D
1639!
1640! 2D V-momentum.
1641!
1642 DO j=jstrp,jendt
1643 DO i=istrt,iendt
1644 tl_vbar(i,j,lout)=d_vbar(i,j)
1645# ifdef MASKING
1646 tl_vbar(i,j,lout)=tl_vbar(i,j,lout)*vmask(i,j)
1647# endif
1648 END DO
1649 END DO
1650# endif
1651
1652# ifdef ADJUST_BOUNDARY
1653!
1654! 2D V-momentum open boundaries.
1655!
1656 IF (any(lobc(:,isvbar,ng))) THEN
1657 DO ir=1,nbrec(ng)
1658 IF ((lobc(iwest,isvbar,ng)).and. &
1659 & domain(ng)%Western_Edge(tile)) THEN
1660 ib=iwest
1661 DO j=jstrv,jend
1662 tl_vbar_obc(j,ib,ir,lout)=d_vbar_obc(j,ib,ir)
1663# ifdef MASKING
1664 tl_vbar_obc(j,ib,ir,lout)=tl_vbar_obc(j,ib,ir,lout)* &
1665 & vmask(istr-1,j)
1666# endif
1667 END DO
1668 END IF
1669 IF ((lobc(ieast,isvbar,ng)).and. &
1670 & domain(ng)%Eastern_Edge(tile)) THEN
1671 ib=ieast
1672 DO j=jstrv,jend
1673 tl_vbar_obc(j,ib,ir,lout)=d_vbar_obc(j,ib,ir)
1674# ifdef MASKING
1675 tl_vbar_obc(j,ib,ir,lout)=tl_vbar_obc(j,ib,ir,lout)* &
1676 & vmask(iend+1,j)
1677# endif
1678 END DO
1679 END IF
1680 IF ((lobc(isouth,isvbar,ng)).and. &
1681 & domain(ng)%Southern_Edge(tile)) THEN
1682 ib=isouth
1683 DO i=istr,iend
1684 tl_vbar_obc(i,ib,ir,lout)=d_vbar_obc(i,ib,ir)
1685# ifdef MASKING
1686 tl_vbar_obc(i,ib,ir,lout)=tl_vbar_obc(i,ib,ir,lout)* &
1687 & vmask(i,jstr)
1688# endif
1689 END DO
1690 END IF
1691 IF ((lobc(inorth,isvbar,ng)).and. &
1692 & domain(ng)%Northern_Edge(tile)) THEN
1693 ib=inorth
1694 DO i=istr,iend
1695 tl_vbar_obc(i,ib,ir,lout)=d_vbar_obc(i,ib,ir)
1696# ifdef MASKING
1697 tl_vbar_obc(i,ib,ir,lout)=tl_vbar_obc(i,ib,ir,lout)* &
1698 & vmask(i,jend+1)
1699# endif
1700 END DO
1701 END IF
1702 END DO
1703 END IF
1704# endif
1705
1706# ifdef ADJUST_WSTRESS
1707!
1708! Surface momentum stress.
1709!
1710 DO ir=1,nfrec(ng)
1711 DO j=jstrt,jendt
1712 DO i=istrp,iendt
1713 tl_ustr(i,j,ir,lout)=d_sustr(i,j,ir)
1714# ifdef MASKING
1715 tl_ustr(i,j,ir,lout)=tl_ustr(i,j,ir,lout)*umask(i,j)
1716# endif
1717 END DO
1718 END DO
1719 DO j=jstrp,jendt
1720 DO i=istrt,iendt
1721 tl_vstr(i,j,ir,lout)=d_svstr(i,j,ir)
1722# ifdef MASKING
1723 tl_vstr(i,j,ir,lout)=tl_vstr(i,j,ir,lout)*vmask(i,j)
1724# endif
1725 END DO
1726 END DO
1727 END DO
1728# endif
1729
1730# ifdef SOLVE3D
1731!
1732! 3D U-momentum.
1733!
1734 DO k=1,n(ng)
1735 DO j=jstrt,jendt
1736 DO i=istrp,iendt
1737 tl_u(i,j,k,lout)=d_u(i,j,k)
1738# ifdef MASKING
1739 tl_u(i,j,k,lout)=tl_u(i,j,k,lout)*umask(i,j)
1740# endif
1741 END DO
1742 END DO
1743 END DO
1744
1745# ifdef ADJUST_BOUNDARY
1746!
1747! 3D U-momentum open boundaries.
1748!
1749 IF (any(lobc(:,isuvel,ng))) THEN
1750 DO ir=1,nbrec(ng)
1751 IF ((lobc(iwest,isuvel,ng)).and. &
1752 & domain(ng)%Western_Edge(tile)) THEN
1753 ib=iwest
1754 DO k=1,n(ng)
1755 DO j=jstr,jend
1756 tl_u_obc(j,k,ib,ir,lout)=d_u_obc(j,k,ib,ir)
1757# ifdef MASKING
1758 tl_u_obc(j,k,ib,ir,lout)=tl_u_obc(j,k,ib,ir,lout)* &
1759 & umask(istr,j)
1760# endif
1761 END DO
1762 END DO
1763 END IF
1764 IF ((lobc(ieast,isuvel,ng)).and. &
1765 & domain(ng)%Eastern_Edge(tile)) THEN
1766 ib=ieast
1767 DO k=1,n(ng)
1768 DO j=jstr,jend
1769 tl_u_obc(j,k,ib,ir,lout)=d_u_obc(j,k,ib,ir)
1770# ifdef MASKING
1771 tl_u_obc(j,k,ib,ir,lout)=tl_u_obc(j,k,ib,ir,lout)* &
1772 & umask(iend+1,j)
1773# endif
1774 END DO
1775 END DO
1776 END IF
1777 IF ((lobc(isouth,isuvel,ng)).and. &
1778 & domain(ng)%Southern_Edge(tile)) THEN
1779 ib=isouth
1780 DO k=1,n(ng)
1781 DO i=istru,iend
1782 tl_u_obc(i,k,ib,ir,lout)=d_u_obc(i,k,ib,ir)
1783# ifdef MASKING
1784 tl_u_obc(i,k,ib,ir,lout)=tl_u_obc(i,k,ib,ir,lout)* &
1785 & umask(i,jstr-1)
1786# endif
1787 END DO
1788 END DO
1789 END IF
1790 IF ((lobc(inorth,isuvel,ng)).and. &
1791 & domain(ng)%Northern_Edge(tile)) THEN
1792 ib=inorth
1793 DO k=1,n(ng)
1794 DO i=istru,iend
1795 tl_u_obc(i,k,ib,ir,lout)=d_u_obc(i,k,ib,ir)
1796# ifdef MASKING
1797 tl_u_obc(i,k,ib,ir,lout)=tl_u_obc(i,k,ib,ir,lout)* &
1798 & umask(i,jend+1)
1799# endif
1800 END DO
1801 END DO
1802 END IF
1803 END DO
1804 END IF
1805# endif
1806!
1807! 3D V-momentum.
1808!
1809 DO k=1,n(ng)
1810 DO j=jstrp,jendt
1811 DO i=istrt,iendt
1812 tl_v(i,j,k,lout)=d_v(i,j,k)
1813# ifdef MASKING
1814 tl_v(i,j,k,lout)=tl_v(i,j,k,lout)*vmask(i,j)
1815# endif
1816 END DO
1817 END DO
1818 END DO
1819
1820# ifdef ADJUST_BOUNDARY
1821!
1822! 3D V-momentum open boundaries.
1823!
1824 IF (any(lobc(:,isvvel,ng))) THEN
1825 DO ir=1,nbrec(ng)
1826 IF ((lobc(iwest,isvvel,ng)).and. &
1827 & domain(ng)%Western_Edge(tile)) THEN
1828 ib=iwest
1829 DO k=1,n(ng)
1830 DO j=jstrv,jend
1831 tl_v_obc(j,k,ib,ir,lout)=d_v_obc(j,k,ib,ir)
1832# ifdef MASKING
1833 tl_v_obc(j,k,ib,ir,lout)=tl_v_obc(j,k,ib,ir,lout)* &
1834 & vmask(istr-1,j)
1835# endif
1836 END DO
1837 END DO
1838 END IF
1839 IF ((lobc(ieast,isvvel,ng)).and. &
1840 & domain(ng)%Eastern_Edge(tile)) THEN
1841 ib=ieast
1842 DO k=1,n(ng)
1843 DO j=jstrv,jend
1844 tl_v_obc(j,k,ib,ir,lout)=d_v_obc(j,k,ib,ir)
1845# ifdef MASKING
1846 tl_v_obc(j,k,ib,ir,lout)=tl_v_obc(j,k,ib,ir,lout)* &
1847 & vmask(iend+1,j)
1848# endif
1849 END DO
1850 END DO
1851 END IF
1852 IF ((lobc(isouth,isvvel,ng)).and. &
1853 & domain(ng)%Southern_Edge(tile)) THEN
1854 ib=isouth
1855 DO k=1,n(ng)
1856 DO i=istr,iend
1857 tl_v_obc(i,k,ib,ir,lout)=d_v_obc(i,k,ib,ir)
1858# ifdef MASKING
1859 tl_v_obc(i,k,ib,ir,lout)=tl_v_obc(i,k,ib,ir,lout)* &
1860 & vmask(i,jstr)
1861# endif
1862 END DO
1863 END DO
1864 END IF
1865 IF ((lobc(inorth,isvvel,ng)).and. &
1866 & domain(ng)%Northern_Edge(tile)) THEN
1867 ib=inorth
1868 DO k=1,n(ng)
1869 DO i=istr,iend
1870 tl_v_obc(i,k,ib,ir,lout)=d_v_obc(i,k,ib,ir)
1871# ifdef MASKING
1872 tl_v_obc(i,k,ib,ir,lout)=tl_v_obc(i,k,ib,ir,lout)* &
1873 & vmask(i,jend+1)
1874# endif
1875 END DO
1876 END DO
1877 END IF
1878 END DO
1879 END IF
1880# endif
1881!
1882! Tracers.
1883!
1884 DO it=1,nt(ng)
1885 DO k=1,n(ng)
1886 DO j=jstrt,jendt
1887 DO i=istrt,iendt
1888 tl_t(i,j,k,lout,it)=d_t(i,j,k,it)
1889# ifdef MASKING
1890 tl_t(i,j,k,lout,it)=tl_t(i,j,k,lout,it)*rmask(i,j)
1891# endif
1892 END DO
1893 END DO
1894 END DO
1895 END DO
1896
1897# ifdef ADJUST_BOUNDARY
1898!
1899! Tracers open boundaries.
1900!
1901 DO it=1,nt(ng)
1902 IF (any(lobc(:,istvar(it),ng))) THEN
1903 DO ir=1,nbrec(ng)
1904 IF ((lobc(iwest,istvar(it),ng)).and. &
1905 & domain(ng)%Western_Edge(tile)) THEN
1906 ib=iwest
1907 DO k=1,n(ng)
1908 DO j=jstr,jend
1909 tl_t_obc(j,k,ib,ir,lout,it)=d_t_obc(j,k,ib,ir,it)
1910# ifdef MASKING
1911 tl_t_obc(j,k,ib,ir,lout,it)= &
1912 & tl_t_obc(j,k,ib,ir,lout,it)*rmask(istr-1,j)
1913# endif
1914 END DO
1915 END DO
1916 END IF
1917 IF ((lobc(ieast,istvar(it),ng)).and. &
1918 & domain(ng)%Eastern_Edge(tile)) THEN
1919 ib=ieast
1920 DO k=1,n(ng)
1921 DO j=jstr,jend
1922 tl_t_obc(j,k,ib,ir,lout,it)=d_t_obc(j,k,ib,ir,it)
1923# ifdef MASKING
1924 tl_t_obc(j,k,ib,ir,lout,it)= &
1925 & tl_t_obc(j,k,ib,ir,lout,it)*rmask(iend+1,j)
1926# endif
1927 END DO
1928 END DO
1929 END IF
1930 IF ((lobc(isouth,istvar(it),ng)).and. &
1931 & domain(ng)%Southern_Edge(tile)) THEN
1932 ib=isouth
1933 DO k=1,n(ng)
1934 DO i=istr,iend
1935 tl_t_obc(i,k,ib,ir,lout,it)=d_t_obc(i,k,ib,ir,it)
1936# ifdef MASKING
1937 tl_t_obc(i,k,ib,ir,lout,it)= &
1938 & tl_t_obc(i,k,ib,ir,lout,it)*rmask(i,jstr-1)
1939# endif
1940 END DO
1941 END DO
1942 END IF
1943 IF ((lobc(inorth,istvar(it),ng)).and. &
1944 & domain(ng)%Northern_Edge(tile)) THEN
1945 ib=inorth
1946 DO k=1,n(ng)
1947 DO i=istr,iend
1948 tl_t_obc(i,k,ib,ir,lout,it)=d_t_obc(i,k,ib,ir,it)
1949# ifdef MASKING
1950 tl_t_obc(i,k,ib,ir,lout,it)= &
1951 & tl_t_obc(i,k,ib,ir,lout,it)*rmask(i,jend+1)
1952# endif
1953 END DO
1954 END DO
1955 END IF
1956 END DO
1957 END IF
1958 END DO
1959# endif
1960
1961# ifdef ADJUST_STFLUX
1962!
1963! Surface tracers flux.
1964!
1965 DO it=1,nt(ng)
1966 IF (lstflux(it,ng)) THEN
1967 DO ir=1,nfrec(ng)
1968 DO j=jstrt,jendt
1969 DO i=istrt,iendt
1970 tl_tflux(i,j,ir,lout,it)=d_stflx(i,j,ir,it)
1971# ifdef MASKING
1972 tl_tflux(i,j,ir,lout,it)=tl_tflux(i,j,ir,lout,it)* &
1973 & rmask(i,j)
1974# endif
1975 END DO
1976 END DO
1977 END DO
1978 END IF
1979 END DO
1980# endif
1981
1982# endif
1983!
1984 RETURN
1985 END SUBROUTINE tl_new_vector
1986!
1987!***********************************************************************
1988 SUBROUTINE new_direction (ng, tile, model, &
1989 & LBi, UBi, LBj, UBj, LBij, UBij, &
1990 & IminS, ImaxS, JminS, JmaxS, &
1991 & Lold, Lnew, &
1992# ifdef MASKING
1993 & rmask, umask, vmask, &
1994# endif
1995# ifdef ADJUST_BOUNDARY
1996# ifdef SOLVE3D
1997 & ad_t_obc, ad_u_obc, ad_v_obc, &
1998# endif
1999 & ad_ubar_obc, ad_vbar_obc, &
2000 & ad_zeta_obc, &
2001# endif
2002# ifdef ADJUST_WSTRESS
2003 & ad_ustr, ad_vstr, &
2004# endif
2005# ifdef SOLVE3D
2006# ifdef ADJUST_STFLUX
2007 & ad_tflux, &
2008# endif
2009 & ad_t, ad_u, ad_v, &
2010# else
2011 & ad_ubar, ad_vbar, &
2012# endif
2013 & ad_zeta, &
2014# ifdef ADJUST_BOUNDARY
2015# ifdef SOLVE3D
2016 & d_t_obc, d_u_obc, d_v_obc, &
2017# endif
2018 & d_ubar_obc, d_vbar_obc, &
2019 & d_zeta_obc, &
2020# endif
2021# ifdef ADJUST_WSTRESS
2022 & d_sustr, d_svstr, &
2023# endif
2024# ifdef SOLVE3D
2025# ifdef ADJUST_STFLUX
2026 & d_stflx, &
2027# endif
2028 & d_t, d_u, d_v, &
2029# else
2030 & d_ubar, d_vbar, &
2031# endif
2032 & d_zeta)
2033!***********************************************************************
2034!
2035! Imported variable declarations.
2036!
2037 integer, intent(in) :: ng, tile, model
2038 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
2039 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
2040 integer, intent(in) :: Lold, Lnew
2041!
2042# ifdef ASSUMED_SHAPE
2043# ifdef MASKING
2044 real(r8), intent(in) :: rmask(LBi:,LBj:)
2045 real(r8), intent(in) :: umask(LBi:,LBj:)
2046 real(r8), intent(in) :: vmask(LBi:,LBj:)
2047# endif
2048# ifdef ADJUST_BOUNDARY
2049# ifdef SOLVE3D
2050 real(r8), intent(in) :: ad_t_obc(LBij:,:,:,:,:,:)
2051 real(r8), intent(in) :: ad_u_obc(LBij:,:,:,:,:)
2052 real(r8), intent(in) :: ad_v_obc(LBij:,:,:,:,:)
2053# endif
2054 real(r8), intent(in) :: ad_ubar_obc(LBij:,:,:,:)
2055 real(r8), intent(in) :: ad_vbar_obc(LBij:,:,:,:)
2056 real(r8), intent(in) :: ad_zeta_obc(LBij:,:,:,:)
2057# endif
2058# ifdef ADJUST_WSTRESS
2059 real(r8), intent(in) :: ad_ustr(LBi:,LBj:,:,:)
2060 real(r8), intent(in) :: ad_vstr(LBi:,LBj:,:,:)
2061# endif
2062# ifdef SOLVE3D
2063# ifdef ADJUST_STFLUX
2064 real(r8), intent(in) :: ad_tflux(LBi:,LBj:,:,:,:)
2065# endif
2066 real(r8), intent(in) :: ad_t(LBi:,LBj:,:,:,:)
2067 real(r8), intent(in) :: ad_u(LBi:,LBj:,:,:)
2068 real(r8), intent(in) :: ad_v(LBi:,LBj:,:,:)
2069# else
2070 real(r8), intent(in) :: ad_ubar(LBi:,LBj:,:)
2071 real(r8), intent(in) :: ad_vbar(LBi:,LBj:,:)
2072# endif
2073 real(r8), intent(in) :: ad_zeta(LBi:,LBj:,:)
2074# ifdef ADJUST_BOUNDARY
2075# ifdef SOLVE3D
2076 real(r8), intent(inout) :: d_t_obc(LBij:,:,:,:,:)
2077 real(r8), intent(inout) :: d_u_obc(LBij:,:,:,:)
2078 real(r8), intent(inout) :: d_v_obc(LBij:,:,:,:)
2079# endif
2080 real(r8), intent(inout) :: d_ubar_obc(LBij:,:,:)
2081 real(r8), intent(inout) :: d_vbar_obc(LBij:,:,:)
2082 real(r8), intent(inout) :: d_zeta_obc(LBij:,:,:)
2083# endif
2084# ifdef ADJUST_WSTRESS
2085 real(r8), intent(inout) :: d_sustr(LBi:,LBj:,:)
2086 real(r8), intent(inout) :: d_svstr(LBi:,LBj:,:)
2087# endif
2088# ifdef SOLVE3D
2089# ifdef ADJUST_STFLUX
2090 real(r8), intent(inout) :: d_stflx(LBi:,LBj:,:,:)
2091# endif
2092 real(r8), intent(inout) :: d_t(LBi:,LBj:,:,:)
2093 real(r8), intent(inout) :: d_u(LBi:,LBj:,:)
2094 real(r8), intent(inout) :: d_v(LBi:,LBj:,:)
2095# else
2096 real(r8), intent(inout) :: d_ubar(LBi:,LBj:)
2097 real(r8), intent(inout) :: d_vbar(LBi:,LBj:)
2098# endif
2099 real(r8), intent(inout) :: d_zeta(LBi:,LBj:)
2100
2101# else
2102
2103# ifdef MASKING
2104 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
2105 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
2106 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
2107# endif
2108# ifdef ADJUST_BOUNDARY
2109# ifdef SOLVE3D
2110 real(r8), intent(in) :: ad_t_obc(LBij:UBij,N(ng),4, &
2111 & Nbrec(ng),2,NT(ng))
2112 real(r8), intent(in) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2113 real(r8), intent(in) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2114# endif
2115 real(r8), intent(in) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
2116 real(r8), intent(in) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
2117 real(r8), intent(in) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
2118# endif
2119# ifdef ADJUST_WSTRESS
2120 real(r8), intent(in) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
2121 real(r8), intent(in) :: ad_vstr(LBi:UBI,LBj:UBj,Nfrec(ng),2)
2122# endif
2123# ifdef SOLVE3D
2124# ifdef ADJUST_STFLUX
2125 real(r8), intent(in) :: ad_tflux(LBi:UBi,LBj:UBj, &
2126 & Nfrec(ng),2,NT(ng))
2127# endif
2128 real(r8), intent(in) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2129 real(r8), intent(in) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
2130 real(r8), intent(in) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
2131# else
2132 real(r8), intent(in) :: ad_ubar(LBi:UBi,LBj:UBj,:)
2133 real(r8), intent(in) :: ad_vbar(LBi:UBi,LBj:UBj,:)
2134# endif
2135 real(r8), intent(in) :: ad_zeta(LBi:UBi,LBj:UBj,:)
2136# ifdef ADJUST_BOUNDARY
2137# ifdef SOLVE3D
2138 real(r8), intent(inout) :: d_t_obc(LBij:UBij,N(ng),4, &
2139 & Nbrec(ng),NT(ng))
2140 real(r8), intent(inout) :: d_u_obc(LBij:UBij,N(ng),4,Nbrec(ng))
2141 real(r8), intent(inout) :: d_v_obc(LBij:UBij,N(ng),4,Nbrec(ng))
2142# endif
2143 real(r8), intent(inout) :: d_ubar_obc(LBij:UBij,4,Nbrec(ng))
2144 real(r8), intent(inout) :: d_vbar_obc(LBij:UBij,4,Nbrec(ng))
2145 real(r8), intent(inout) :: d_zeta_obc(LBij:UBij,4,Nbrec(ng))
2146# endif
2147# ifdef ADJUST_WSTRESS
2148 real(r8), intent(inout) :: d_sustr(LBi:UBi,LBj:UBj,Nfrec(ng))
2149 real(r8), intent(inout) :: d_svstr(LBi:UBI,LBj:UBj,Nfrec(ng))
2150# endif
2151# ifdef SOLVE3D
2152# ifdef ADJUST_STFLUX
2153 real(r8), intent(inout) :: d_stflx(LBi:UBi,LBj:UBj, &
2154 & Nfrec(ng),NT(ng))
2155# endif
2156 real(r8), intent(inout) :: d_t(LBi:UBi,LBj:UBj,N(ng),NT(ng))
2157 real(r8), intent(inout) :: d_u(LBi:UBi,LBj:UBj,N(ng))
2158 real(r8), intent(inout) :: d_v(LBi:UBi,LBj:UBj,N(ng))
2159# else
2160 real(r8), intent(inout) :: d_ubar(LBi:UBi,LBj:UBj)
2161 real(r8), intent(inout) :: d_vbar(LBi:UBi,LBj:UBj)
2162# endif
2163 real(r8), intent(inout) :: d_zeta(LBi:UBi,LBj:UBj)
2164# endif
2165!
2166! Local variable declarations.
2167!
2168 integer :: i, j, k
2169 integer :: ib, ir, it
2170
2171# include "set_bounds.h"
2172!
2173!-----------------------------------------------------------------------
2174! Compute new conjugate descent direction, d(k+1). Notice that the old
2175! descent direction is overwritten.
2176!-----------------------------------------------------------------------
2177!
2178! Free-sruface.
2179!
2180 DO j=jstrt,jendt
2181 DO i=istrt,iendt
2182 d_zeta(i,j)=ad_zeta(i,j,lnew)
2183# ifdef MASKING
2184 d_zeta(i,j)=d_zeta(i,j)*rmask(i,j)
2185# endif
2186 END DO
2187 END DO
2188
2189# ifdef ADJUST_BOUNDARY
2190!
2191! Free-surface open boundaries.
2192!
2193 IF (any(lobc(:,isfsur,ng))) THEN
2194 DO ir=1,nbrec(ng)
2195 IF ((lobc(iwest,isfsur,ng)).and. &
2196 & domain(ng)%Western_Edge(tile)) THEN
2197 ib=iwest
2198 DO j=jstr,jend
2199 d_zeta_obc(j,ib,ir)=ad_zeta_obc(j,ib,ir,lnew)
2200# ifdef MASKING
2201 d_zeta_obc(j,ib,ir)=d_zeta_obc(j,ib,ir)* &
2202 & rmask(istr-1,j)
2203# endif
2204 END DO
2205 END IF
2206 IF ((lobc(ieast,isfsur,ng)).and. &
2207 & domain(ng)%Eastern_Edge(tile)) THEN
2208 ib=ieast
2209 DO j=jstr,jend
2210 d_zeta_obc(j,ib,ir)=ad_zeta_obc(j,ib,ir,lnew)
2211# ifdef MASKING
2212 d_zeta_obc(j,ib,ir)=d_zeta_obc(j,ib,ir)* &
2213 & rmask(iend+1,j)
2214# endif
2215 END DO
2216 END IF
2217 IF ((lobc(isouth,isfsur,ng)).and. &
2218 & domain(ng)%Southern_Edge(tile)) THEN
2219 ib=isouth
2220 DO i=istr,iend
2221 d_zeta_obc(i,ib,ir)=ad_zeta_obc(i,ib,ir,lnew)
2222# ifdef MASKING
2223 d_zeta_obc(i,ib,ir)=d_zeta_obc(i,ib,ir)* &
2224 & rmask(i,jstr-1)
2225# endif
2226 END DO
2227 END IF
2228 IF ((lobc(inorth,isfsur,ng)).and. &
2229 & domain(ng)%Northern_Edge(tile)) THEN
2230 ib=inorth
2231 DO i=istr,iend
2232 d_zeta_obc(i,ib,ir)=ad_zeta_obc(i,ib,ir,lnew)
2233# ifdef MASKING
2234 d_zeta_obc(i,ib,ir)=d_zeta_obc(i,ib,ir)* &
2235 & rmask(i,jend+1)
2236# endif
2237 END DO
2238 END IF
2239 END DO
2240 END IF
2241# endif
2242
2243# ifndef SOLVE3D
2244!
2245! 2D U-momentum.
2246!
2247 DO j=jstrt,jendt
2248 DO i=istrp,iendt
2249 d_ubar(i,j)=ad_ubar(i,j,lnew)
2250# ifdef MASKING
2251 d_ubar(i,j)=d_ubar(i,j)*umask(i,j)
2252# endif
2253 END DO
2254 END DO
2255# endif
2256
2257# ifdef ADJUST_BOUNDARY
2258!
2259! 2D U-momentum open boundaries.
2260!
2261 IF (any(lobc(:,isubar,ng))) THEN
2262 DO ir=1,nbrec(ng)
2263 IF ((lobc(iwest,isubar,ng)).and. &
2264 & domain(ng)%Western_Edge(tile)) THEN
2265 ib=iwest
2266 DO j=jstr,jend
2267 d_ubar_obc(j,ib,ir)=ad_ubar_obc(j,ib,ir,lnew)
2268# ifdef MASKING
2269 d_ubar_obc(j,ib,ir)=d_ubar_obc(j,ib,ir)* &
2270 & umask(istr,j)
2271# endif
2272 END DO
2273 END IF
2274 IF ((lobc(ieast,isubar,ng)).and. &
2275 & domain(ng)%Eastern_Edge(tile)) THEN
2276 ib=ieast
2277 DO j=jstr,jend
2278 d_ubar_obc(j,ib,ir)=ad_ubar_obc(j,ib,ir,lnew)
2279# ifdef MASKING
2280 d_ubar_obc(j,ib,ir)=d_ubar_obc(j,ib,ir)* &
2281 & umask(iend+1,j)
2282# endif
2283 END DO
2284 END IF
2285 IF ((lobc(isouth,isubar,ng)).and. &
2286 & domain(ng)%Southern_Edge(tile)) THEN
2287 ib=isouth
2288 DO i=istru,iend
2289 d_ubar_obc(i,ib,ir)=ad_ubar_obc(i,ib,ir,lnew)
2290# ifdef MASKING
2291 d_ubar_obc(i,ib,ir)=d_ubar_obc(i,ib,ir)* &
2292 & umask(i,jstr-1)
2293# endif
2294 END DO
2295 END IF
2296 IF ((lobc(inorth,isubar,ng)).and. &
2297 & domain(ng)%Northern_Edge(tile)) THEN
2298 ib=inorth
2299 DO i=istru,iend
2300 d_ubar_obc(i,ib,ir)=ad_ubar_obc(i,ib,ir,lnew)
2301# ifdef MASKING
2302 d_ubar_obc(i,ib,ir)=d_ubar_obc(i,ib,ir)* &
2303 & umask(i,jend+1)
2304# endif
2305 END DO
2306 END IF
2307 END DO
2308 END IF
2309# endif
2310
2311# ifndef SOLVE3D
2312!
2313! 2D V-momentum.
2314!
2315 DO j=jstrp,jendt
2316 DO i=istrt,iendt
2317 d_vbar(i,j)=ad_vbar(i,j,lnew)
2318# ifdef MASKING
2319 d_vbar(i,j)=d_vbar(i,j)*vmask(i,j)
2320# endif
2321 END DO
2322 END DO
2323# endif
2324
2325# ifdef ADJUST_BOUNDARY
2326!
2327! 2D V-momentum open boundaries.
2328!
2329 IF (any(lobc(:,isvbar,ng))) THEN
2330 DO ir=1,nbrec(ng)
2331 IF ((lobc(iwest,isvbar,ng)).and. &
2332 & domain(ng)%Western_Edge(tile)) THEN
2333 ib=iwest
2334 DO j=jstrv,jend
2335 d_vbar_obc(j,ib,ir)=ad_vbar_obc(j,ib,ir,lnew)
2336# ifdef MASKING
2337 d_vbar_obc(j,ib,ir)=d_vbar_obc(j,ib,ir)* &
2338 & vmask(istr-1,j)
2339# endif
2340 END DO
2341 END IF
2342 IF ((lobc(ieast,isvbar,ng)).and. &
2343 & domain(ng)%Eastern_Edge(tile)) THEN
2344 ib=ieast
2345 DO j=jstrv,jend
2346 d_vbar_obc(j,ib,ir)=ad_vbar_obc(j,ib,ir,lnew)
2347# ifdef MASKING
2348 d_vbar_obc(j,ib,ir)=d_vbar_obc(j,ib,ir)* &
2349 & vmask(iend+1,j)
2350# endif
2351 END DO
2352 END IF
2353 IF ((lobc(isouth,isvbar,ng)).and. &
2354 & domain(ng)%Southern_Edge(tile)) THEN
2355 ib=isouth
2356 DO i=istr,iend
2357 d_vbar_obc(i,ib,ir)=ad_vbar_obc(i,ib,ir,lnew)
2358# ifdef MASKING
2359 d_vbar_obc(i,ib,ir)=d_vbar_obc(i,ib,ir)* &
2360 & vmask(i,jstr)
2361# endif
2362 END DO
2363 END IF
2364 IF ((lobc(inorth,isvbar,ng)).and. &
2365 & domain(ng)%Northern_Edge(tile)) THEN
2366 ib=inorth
2367 DO i=istr,iend
2368 d_vbar_obc(i,ib,ir)=ad_vbar_obc(i,ib,ir,lnew)
2369# ifdef MASKING
2370 d_vbar_obc(i,ib,ir)=d_vbar_obc(i,ib,ir)* &
2371 & vmask(i,jend+1)
2372# endif
2373 END DO
2374 END IF
2375 END DO
2376 END IF
2377# endif
2378
2379# ifdef ADJUST_WSTRESS
2380!
2381! Surface momentum stress.
2382!
2383 DO ir=1,nfrec(ng)
2384 DO j=jstrt,jendt
2385 DO i=istrp,iendt
2386 d_sustr(i,j,ir)=ad_ustr(i,j,ir,lnew)
2387# ifdef MASKING
2388 d_sustr(i,j,ir)=d_sustr(i,j,ir)*umask(i,j)
2389# endif
2390 END DO
2391 END DO
2392 DO j=jstrp,jendt
2393 DO i=istrt,iendt
2394 d_svstr(i,j,ir)=ad_vstr(i,j,ir,lnew)
2395# ifdef MASKING
2396 d_svstr(i,j,ir)=d_svstr(i,j,ir)*vmask(i,j)
2397# endif
2398 END DO
2399 END DO
2400 END DO
2401# endif
2402
2403# ifdef SOLVE3D
2404!
2405! 3D U-momentum.
2406!
2407 DO k=1,n(ng)
2408 DO j=jstrt,jendt
2409 DO i=istrp,iendt
2410 d_u(i,j,k)=ad_u(i,j,k,lnew)
2411# ifdef MASKING
2412 d_u(i,j,k)=d_u(i,j,k)*umask(i,j)
2413# endif
2414 END DO
2415 END DO
2416 END DO
2417
2418# ifdef ADJUST_BOUNDARY
2419!
2420! 3D U-momentum open boundaries.
2421!
2422 IF (any(lobc(:,isuvel,ng))) THEN
2423 DO ir=1,nbrec(ng)
2424 IF ((lobc(iwest,isuvel,ng)).and. &
2425 & domain(ng)%Western_Edge(tile)) THEN
2426 ib=iwest
2427 DO k=1,n(ng)
2428 DO j=jstr,jend
2429 d_u_obc(j,k,ib,ir)=ad_u_obc(j,k,ib,ir,lnew)
2430# ifdef MASKING
2431 d_u_obc(j,k,ib,ir)=d_u_obc(j,k,ib,ir)* &
2432 & umask(istr,j)
2433# endif
2434 END DO
2435 END DO
2436 END IF
2437 IF ((lobc(ieast,isuvel,ng)).and. &
2438 & domain(ng)%Eastern_Edge(tile)) THEN
2439 ib=ieast
2440 DO k=1,n(ng)
2441 DO j=jstr,jend
2442 d_u_obc(j,k,ib,ir)=ad_u_obc(j,k,ib,ir,lnew)
2443# ifdef MASKING
2444 d_u_obc(j,k,ib,ir)=d_u_obc(j,k,ib,ir)* &
2445 & umask(iend+1,j)
2446# endif
2447 END DO
2448 END DO
2449 END IF
2450 IF ((lobc(isouth,isuvel,ng)).and. &
2451 & domain(ng)%Southern_Edge(tile)) THEN
2452 ib=isouth
2453 DO k=1,n(ng)
2454 DO i=istru,iend
2455 d_u_obc(i,k,ib,ir)=ad_u_obc(i,k,ib,ir,lnew)
2456# ifdef MASKING
2457 d_u_obc(i,k,ib,ir)=d_u_obc(i,k,ib,ir)* &
2458 & umask(i,jstr-1)
2459# endif
2460 END DO
2461 END DO
2462 END IF
2463 IF ((lobc(inorth,isuvel,ng)).and. &
2464 & domain(ng)%Northern_Edge(tile)) THEN
2465 ib=inorth
2466 DO k=1,n(ng)
2467 DO i=istru,iend
2468 d_u_obc(i,k,ib,ir)=ad_u_obc(i,k,ib,ir,lnew)
2469# ifdef MASKING
2470 d_u_obc(i,k,ib,ir)=d_u_obc(i,k,ib,ir)* &
2471 & umask(i,jend+1)
2472# endif
2473 END DO
2474 END DO
2475 END IF
2476 END DO
2477 END IF
2478# endif
2479!
2480! 3D V-momentum.
2481!
2482 DO k=1,n(ng)
2483 DO j=jstrp,jendt
2484 DO i=istrt,iendt
2485 d_v(i,j,k)=ad_v(i,j,k,lnew)
2486# ifdef MASKING
2487 d_v(i,j,k)=d_v(i,j,k)*vmask(i,j)
2488# endif
2489 END DO
2490 END DO
2491 END DO
2492
2493# ifdef ADJUST_BOUNDARY
2494!
2495! 3D V-momentum open boundaries.
2496!
2497 IF (any(lobc(:,isvvel,ng))) THEN
2498 DO ir=1,nbrec(ng)
2499 IF ((lobc(iwest,isvvel,ng)).and. &
2500 & domain(ng)%Western_Edge(tile)) THEN
2501 ib=iwest
2502 DO k=1,n(ng)
2503 DO j=jstrv,jend
2504 d_v_obc(j,k,ib,ir)=ad_v_obc(j,k,ib,ir,lnew)
2505# ifdef MASKING
2506 d_v_obc(j,k,ib,ir)=d_v_obc(j,k,ib,ir)* &
2507 & vmask(istr-1,j)
2508# endif
2509 END DO
2510 END DO
2511 END IF
2512 IF ((lobc(ieast,isvvel,ng)).and. &
2513 & domain(ng)%Eastern_Edge(tile)) THEN
2514 ib=ieast
2515 DO k=1,n(ng)
2516 DO j=jstrv,jend
2517 d_v_obc(j,k,ib,ir)=ad_v_obc(j,k,ib,ir,lnew)
2518# ifdef MASKING
2519 d_v_obc(j,k,ib,ir)=d_v_obc(j,k,ib,ir)* &
2520 & vmask(iend+1,j)
2521# endif
2522 END DO
2523 END DO
2524 END IF
2525 IF ((lobc(isouth,isvvel,ng)).and. &
2526 & domain(ng)%Southern_Edge(tile)) THEN
2527 ib=isouth
2528 DO k=1,n(ng)
2529 DO i=istr,iend
2530 d_v_obc(i,k,ib,ir)=ad_v_obc(i,k,ib,ir,lnew)
2531# ifdef MASKING
2532 d_v_obc(i,k,ib,ir)=d_v_obc(i,k,ib,ir)* &
2533 & vmask(i,jstr)
2534# endif
2535 END DO
2536 END DO
2537 END IF
2538 IF ((lobc(inorth,isvvel,ng)).and. &
2539 & domain(ng)%Northern_Edge(tile)) THEN
2540 ib=inorth
2541 DO k=1,n(ng)
2542 DO i=istr,iend
2543 d_v_obc(i,k,ib,ir)=ad_v_obc(i,k,ib,ir,lnew)
2544# ifdef MASKING
2545 d_v_obc(i,k,ib,ir)=d_v_obc(i,k,ib,ir)* &
2546 & vmask(i,jend+1)
2547# endif
2548 END DO
2549 END DO
2550 END IF
2551 END DO
2552 END IF
2553# endif
2554!
2555! Tracers.
2556!
2557 DO it=1,nt(ng)
2558 DO k=1,n(ng)
2559 DO j=jstrt,jendt
2560 DO i=istrt,iendt
2561 d_t(i,j,k,it)=ad_t(i,j,k,lnew,it)
2562# ifdef MASKING
2563 d_t(i,j,k,it)=d_t(i,j,k,it)*rmask(i,j)
2564# endif
2565 END DO
2566 END DO
2567 END DO
2568 END DO
2569
2570# ifdef ADJUST_BOUNDARY
2571!
2572! Tracers open boundaries.
2573!
2574 DO it=1,nt(ng)
2575 IF (any(lobc(:,istvar(it),ng))) THEN
2576 DO ir=1,nbrec(ng)
2577 IF ((lobc(iwest,istvar(it),ng)).and. &
2578 & domain(ng)%Western_Edge(tile)) THEN
2579 ib=iwest
2580 DO k=1,n(ng)
2581 DO j=jstr,jend
2582 d_t_obc(j,k,ib,ir,it)=ad_t_obc(j,k,ib,ir,lnew,it)
2583# ifdef MASKING
2584 d_t_obc(j,k,ib,ir,it)=d_t_obc(j,k,ib,ir,it)* &
2585 & rmask(istr-1,j)
2586# endif
2587 END DO
2588 END DO
2589 END IF
2590 IF ((lobc(ieast,istvar(it),ng)).and. &
2591 & domain(ng)%Eastern_Edge(tile)) THEN
2592 ib=ieast
2593 DO k=1,n(ng)
2594 DO j=jstr,jend
2595 d_t_obc(j,k,ib,ir,it)=ad_t_obc(j,k,ib,ir,lnew,it)
2596# ifdef MASKING
2597 d_t_obc(j,k,ib,ir,it)=d_t_obc(j,k,ib,ir,it)* &
2598 & rmask(iend+1,j)
2599# endif
2600 END DO
2601 END DO
2602 END IF
2603 IF ((lobc(isouth,istvar(it),ng)).and. &
2604 & domain(ng)%Southern_Edge(tile)) THEN
2605 ib=isouth
2606 DO k=1,n(ng)
2607 DO i=istr,iend
2608 d_t_obc(i,k,ib,ir,it)=ad_t_obc(i,k,ib,ir,lnew,it)
2609# ifdef MASKING
2610 d_t_obc(i,k,ib,ir,it)=d_t_obc(i,k,ib,ir,it)* &
2611 & rmask(i,jstr-1)
2612# endif
2613 END DO
2614 END DO
2615 END IF
2616 IF ((lobc(inorth,istvar(it),ng)).and. &
2617 & domain(ng)%Northern_Edge(tile)) THEN
2618 ib=inorth
2619 DO k=1,n(ng)
2620 DO i=istr,iend
2621 d_t_obc(i,k,ib,ir,it)=ad_t_obc(i,k,ib,ir,lnew,it)
2622# ifdef MASKING
2623 d_t_obc(i,k,ib,ir,it)=d_t_obc(i,k,ib,ir,it)* &
2624 & rmask(i,jend+1)
2625# endif
2626 END DO
2627 END DO
2628 END IF
2629 END DO
2630 END IF
2631 END DO
2632# endif
2633
2634# ifdef ADJUST_STFLUX
2635!
2636! Surface tracers flux.
2637!
2638 DO it=1,nt(ng)
2639 IF (lstflux(it,ng)) THEN
2640 DO ir=1,nfrec(ng)
2641 DO j=jstrt,jendt
2642 DO i=istrt,iendt
2643 d_stflx(i,j,ir,it)=ad_tflux(i,j,ir,lnew,it)
2644# ifdef MASKING
2645 d_stflx(i,j,ir,it)=d_stflx(i,j,ir,it)*rmask(i,j)
2646# endif
2647 END DO
2648 END DO
2649 END DO
2650 END IF
2651 END DO
2652# endif
2653# endif
2654!
2655 RETURN
2656 END SUBROUTINE new_direction
2657!
2658!***********************************************************************
2659 SUBROUTINE analysis_error (ng, tile, model, &
2660 & LBi, UBi, LBj, UBj, LBij, UBij, &
2661 & IminS, ImaxS, JminS, JmaxS, &
2662 & Lold, Lnew, Lwrk, &
2663 & innLoop, outLoop, &
2664# ifdef MASKING
2665 & rmask, umask, vmask, &
2666# endif
2667# ifdef ADJUST_BOUNDARY
2668# ifdef SOLVE3D
2669 & ad_t_obc, ad_u_obc, ad_v_obc, &
2670# endif
2671 & ad_ubar_obc, ad_vbar_obc, &
2672 & ad_zeta_obc, &
2673# endif
2674# ifdef ADJUST_WSTRESS
2675 & ad_ustr, ad_vstr, &
2676# endif
2677# ifdef SOLVE3D
2678# ifdef ADJUST_STFLUX
2679 & ad_tflux, &
2680# endif
2681 & ad_t, ad_u, ad_v, &
2682# if defined WEAK_CONSTRAINT && defined TIME_CONV
2683 & ad_ubar, ad_vbar, &
2684# endif
2685# else
2686 & ad_ubar, ad_vbar, &
2687# endif
2688 & ad_zeta, &
2689# ifdef ADJUST_BOUNDARY
2690# ifdef SOLVE3D
2691 & tl_t_obc, tl_u_obc, tl_v_obc, &
2692# endif
2693 & tl_ubar_obc, tl_vbar_obc, &
2694 & tl_zeta_obc, &
2695# endif
2696# ifdef ADJUST_WSTRESS
2697 & tl_ustr, tl_vstr, &
2698# endif
2699# ifdef SOLVE3D
2700# ifdef ADJUST_STFLUX
2701 & tl_tflux, &
2702# endif
2703 & tl_t, tl_u, tl_v, &
2704# if defined WEAK_CONSTRAINT && defined TIME_CONV
2705 & tl_ubar, tl_vbar, &
2706# endif
2707# else
2708 & tl_ubar, tl_vbar, &
2709# endif
2710 & tl_zeta, &
2711# ifdef ADJUST_BOUNDARY
2712# ifdef SOLVE3D
2713 & nl_t_obc, nl_u_obc, nl_v_obc, &
2714# endif
2715 & nl_ubar_obc, nl_vbar_obc, &
2716 & nl_zeta_obc, &
2717# endif
2718# ifdef ADJUST_WSTRESS
2719 & nl_ustr, nl_vstr, &
2720# endif
2721# ifdef SOLVE3D
2722# ifdef ADJUST_STFLUX
2723 & nl_tflux, &
2724# endif
2725 & nl_t, nl_u, nl_v, &
2726# if defined WEAK_CONSTRAINT && defined TIME_CONV
2727 & nl_ubar, nl_vbar, &
2728# endif
2729# else
2730 & nl_ubar, nl_vbar, &
2731# endif
2732 & nl_zeta)
2733!***********************************************************************
2734!
2735! Imported variable declarations.
2736!
2737 integer, intent(in) :: ng, tile, model
2738 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
2739 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
2740 integer, intent(in) :: Lold, Lnew, Lwrk
2741 integer, intent(in) :: innLoop, outLoop
2742!
2743# ifdef ASSUMED_SHAPE
2744# ifdef MASKING
2745 real(r8), intent(in) :: rmask(LBi:,LBj:)
2746 real(r8), intent(in) :: umask(LBi:,LBj:)
2747 real(r8), intent(in) :: vmask(LBi:,LBj:)
2748# endif
2749# ifdef ADJUST_BOUNDARY
2750# ifdef SOLVE3D
2751 real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
2752 real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
2753 real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
2754# endif
2755 real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
2756 real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
2757 real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
2758# endif
2759# ifdef ADJUST_WSTRESS
2760 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
2761 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
2762# endif
2763# ifdef SOLVE3D
2764# ifdef ADJUST_STFLUX
2765 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
2766# endif
2767 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
2768 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
2769 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
2770# if defined WEAK_CONSTRAINT && defined TIME_CONV
2771 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
2772 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
2773# endif
2774# else
2775 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
2776 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
2777# endif
2778 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
2779# ifdef ADJUST_BOUNDARY
2780# ifdef SOLVE3D
2781 real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
2782 real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
2783 real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
2784# endif
2785 real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
2786 real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
2787 real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
2788# endif
2789# ifdef ADJUST_WSTRESS
2790 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
2791 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
2792# endif
2793# ifdef SOLVE3D
2794# ifdef ADJUST_STFLUX
2795 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
2796# endif
2797 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
2798 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
2799 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
2800# if defined WEAK_CONSTRAINT && defined TIME_CONV
2801 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
2802 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
2803# endif
2804# else
2805 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
2806 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
2807# endif
2808 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
2809# ifdef ADJUST_BOUNDARY
2810# ifdef SOLVE3D
2811 real(r8), intent(inout) :: nl_t_obc(LBij:,:,:,:,:,:)
2812 real(r8), intent(inout) :: nl_u_obc(LBij:,:,:,:,:)
2813 real(r8), intent(inout) :: nl_v_obc(LBij:,:,:,:,:)
2814# endif
2815 real(r8), intent(inout) :: nl_ubar_obc(LBij:,:,:,:)
2816 real(r8), intent(inout) :: nl_vbar_obc(LBij:,:,:,:)
2817 real(r8), intent(inout) :: nl_zeta_obc(LBij:,:,:,:)
2818# endif
2819# ifdef ADJUST_WSTRESS
2820 real(r8), intent(inout) :: nl_ustr(LBi:,LBj:,:,:)
2821 real(r8), intent(inout) :: nl_vstr(LBi:,LBj:,:,:)
2822# endif
2823# ifdef SOLVE3D
2824# ifdef ADJUST_STFLUX
2825 real(r8), intent(inout) :: nl_tflux(LBi:,LBj:,:,:,:)
2826# endif
2827 real(r8), intent(inout) :: nl_t(LBi:,LBj:,:,:,:)
2828 real(r8), intent(inout) :: nl_u(LBi:,LBj:,:,:)
2829 real(r8), intent(inout) :: nl_v(LBi:,LBj:,:,:)
2830# if defined WEAK_CONSTRAINT && defined TIME_CONV
2831 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
2832 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
2833# endif
2834# else
2835 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
2836 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
2837# endif
2838 real(r8), intent(inout) :: nl_zeta(LBi:,LBj:,:)
2839
2840# else
2841
2842# ifdef MASKING
2843 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
2844 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
2845 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
2846# endif
2847# ifdef ADJUST_BOUNDARY
2848# ifdef SOLVE3D
2849 real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, &
2850 & Nbrec(ng),2,NT(ng))
2851 real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2852 real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2853# endif
2854 real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
2855 real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
2856 real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
2857# endif
2858# ifdef ADJUST_WSTRESS
2859 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
2860 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
2861# endif
2862# ifdef SOLVE3D
2863# ifdef ADJUST_STFLUX
2864 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, &
2865 & Nfrec(ng),2,NT(ng))
2866# endif
2867 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2868 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
2869 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
2870# if defined WEAK_CONSTRAINT && defined TIME_CONV
2871 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
2872 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
2873# endif
2874# else
2875 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
2876 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
2877# endif
2878 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
2879# ifdef ADJUST_BOUNDARY
2880# ifdef SOLVE3D
2881 real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, &
2882 & Nbrec(ng),2,NT(ng))
2883 real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2884 real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2885# endif
2886 real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
2887 real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
2888 real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
2889# endif
2890# ifdef ADJUST_WSTRESS
2891 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
2892 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
2893# endif
2894# ifdef SOLVE3D
2895# ifdef ADJUST_STFLUX
2896 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, &
2897 & Nfrec(ng),2,NT(ng))
2898# endif
2899 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2900 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
2901 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
2902# if defined WEAK_CONSTRAINT && defined TIME_CONV
2903 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
2904 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
2905# endif
2906# else
2907 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
2908 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
2909# endif
2910 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
2911# ifdef ADJUST_BOUNDARY
2912# ifdef SOLVE3D
2913 real(r8), intent(inout) :: nl_t_obc(LBij:UBij,N(ng),4, &
2914 & Nbrec(ng),2,NT(ng))
2915 real(r8), intent(inout) :: nl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2916 real(r8), intent(inout) :: nl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
2917# endif
2918 real(r8), intent(inout) :: nl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
2919 real(r8), intent(inout) :: nl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
2920 real(r8), intent(inout) :: nl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
2921# endif
2922# ifdef ADJUST_WSTRESS
2923 real(r8), intent(inout) :: nl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
2924 real(r8), intent(inout) :: nl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
2925# endif
2926# ifdef SOLVE3D
2927# ifdef ADJUST_STFLUX
2928 real(r8), intent(inout) :: nl_tflux(LBi:UBi,LBj:UBj, &
2929 & Nfrec(ng),2,NT(ng))
2930# endif
2931 real(r8), intent(inout) :: nl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2932 real(r8), intent(inout) :: nl_u(LBi:UBi,LBj:UBj,N(ng),2)
2933 real(r8), intent(inout) :: nl_v(LBi:UBi,LBj:UBj,N(ng),2)
2934# if defined WEAK_CONSTRAINT && defined TIME_CONV
2935 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,:)
2936 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,:)
2937# endif
2938# else
2939 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,:)
2940 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,:)
2941# endif
2942 real(r8), intent(inout) :: nl_zeta(LBi:UBi,LBj:UBj,:)
2943# endif
2944!
2945! Local variable declarations.
2946!
2947 integer :: i, j, k, lstr
2948 integer :: ib, ir, it, ivec, Ltmp1, Ltmp2, rec
2949!
2950 real(r8) :: fac, fac1, fac2, zbet
2951
2952 real(r8), dimension(0:NstateVar(ng)) :: dot
2953 real(r8), dimension(1:Ninner) :: DotProd
2954
2955 real(r8), dimension(Ninner) :: zu, zgam
2956!
2957 character (len=256) :: ncname
2958
2959 character (len=*), parameter :: MyFile = &
2960 & __FILE__//", analysis_error"
2961
2962# include "set_bounds.h"
2963!
2964!-----------------------------------------------------------------------
2965! Compute analysis error covariance matrix.
2966!-----------------------------------------------------------------------
2967!
2968! NOTE: In the case of weak constraint ("WEAK_CONSTRAINT") and
2969! and time convolutions ("TIME_CONV"), the state arrays
2970! tl_ubar, tl_vbar, ad_ubar, and ad_vbar are only passed
2971! as required by the "state" operators routines but they
2972! are not used in subsequent calculations.
2973!
2974 ltmp1=1
2975 ltmp2=2
2976!
2977! Compute the dot-product of ad_var(Lnew) with each evolved and
2978! convolved Lanczos vector of the stabilized representer matrix
2979! which are temporarily stored in the hessian NetCDF file.
2980!
2981 ncname=hss(ng)%name
2982 DO ivec=1,ninner
2983 rec=ivec
2984!
2985! Read vectors stored in hessian NetCDF file using nl_var(Ltmp1) as
2986! temporary storage.
2987!
2988 CALL state_read (ng, tile, model, hss(ng)%IOtype, &
2989 & lbi, ubi, lbj, ubj, lbij, ubij, &
2990 & ltmp1, rec, &
2991 & 0, hss(ng)%ncid, ncname, &
2992# ifdef MASKING
2993 & rmask, umask, vmask, &
2994# endif
2995# ifdef ADJUST_BOUNDARY
2996# ifdef SOLVE3D
2997 & nl_t_obc, nl_u_obc, nl_v_obc, &
2998# endif
2999 & nl_ubar_obc, nl_vbar_obc, &
3000 & nl_zeta_obc, &
3001# endif
3002# ifdef ADJUST_WSTRESS
3003 & nl_ustr, nl_vstr, &
3004# endif
3005# ifdef SOLVE3D
3006# ifdef ADJUST_STFLUX
3007 & nl_tflux, &
3008# endif
3009 & nl_t, nl_u, nl_v, &
3010# else
3011 & nl_ubar, nl_vbar, &
3012# endif
3013 & nl_zeta)
3014!
3015! Compute dot product.
3016!
3017 CALL state_dotprod (ng, tile, model, &
3018 & lbi, ubi, lbj, ubj, lbij, ubij, &
3019 & nstatevar(ng), dot(0:), &
3020# ifdef MASKING
3021 & rmask, umask, vmask, &
3022# endif
3023# ifdef ADJUST_BOUNDARY
3024# ifdef SOLVE3D
3025 & ad_t_obc(:,:,:,:,lnew,:), &
3026 & nl_t_obc(:,:,:,:,ltmp1,:), &
3027 & ad_u_obc(:,:,:,:,lnew), &
3028 & nl_u_obc(:,:,:,:,ltmp1), &
3029 & ad_v_obc(:,:,:,:,lnew), &
3030 & nl_v_obc(:,:,:,:,ltmp1), &
3031# endif
3032 & ad_ubar_obc(:,:,:,lnew), &
3033 & nl_ubar_obc(:,:,:,ltmp1), &
3034 & ad_vbar_obc(:,:,:,lnew), &
3035 & nl_vbar_obc(:,:,:,ltmp1), &
3036 & ad_zeta_obc(:,:,:,lnew), &
3037 & nl_zeta_obc(:,:,:,ltmp1), &
3038# endif
3039# ifdef ADJUST_WSTRESS
3040 & ad_ustr(:,:,:,lnew), nl_ustr(:,:,:,ltmp1), &
3041 & ad_vstr(:,:,:,lnew), nl_vstr(:,:,:,ltmp1), &
3042# endif
3043# ifdef SOLVE3D
3044# ifdef ADJUST_STFLUX
3045 & ad_tflux(:,:,:,lnew,:), &
3046 & nl_tflux(:,:,:,ltmp1,:), &
3047# endif
3048 & ad_t(:,:,:,lnew,:), nl_t(:,:,:,ltmp1,:), &
3049 & ad_u(:,:,:,lnew), nl_u(:,:,:,ltmp1), &
3050 & ad_v(:,:,:,lnew), nl_v(:,:,:,ltmp1), &
3051# else
3052 & ad_ubar(:,:,lnew), nl_ubar(:,:,ltmp1), &
3053 & ad_vbar(:,:,lnew), nl_vbar(:,:,ltmp1), &
3054# endif
3055 & ad_zeta(:,:,lnew), nl_zeta(:,:,ltmp1))
3056
3057 dotprod(ivec)=dot(0)
3058 END DO
3059!
3060! Now multiply the result by the inverse tridiagonal matrix of
3061! stabilized representer matrix Lanczos vector coefficients.
3062!
3063 zbet=cg_delta(1,outloop)
3064 zu(1)=dotprod(1)/zbet
3065 DO ivec=2,ninner
3066 zgam(ivec)=cg_beta(ivec,outloop)/zbet
3067 zbet=cg_delta(ivec,outloop)-cg_beta(ivec,outloop)*zgam(ivec)
3068 zu(ivec)=(dotprod(ivec)-cg_beta(ivec,outloop)* &
3069 & zu(ivec-1))/zbet
3070 END DO
3071
3072 DO ivec=ninner-1,1,-1
3073 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
3074 END DO
3075!
3076! Now form the weighted sum of the evolved and convolved Lanczos vector
3077! of the stabilized representer matrix. Use nl_var(2) as temporary
3078! storage.
3079!
3080! Initialize nonl-linear state arrays: nl_var(Ltmp2) = fac
3081!
3082 fac=0.0_r8
3083
3084 CALL state_initialize (ng, tile, &
3085 & lbi, ubi, lbj, ubj, lbij, ubij, &
3086 & ltmp2, fac, &
3087# ifdef MASKING
3088 & rmask, umask, vmask, &
3089# endif
3090# ifdef ADJUST_BOUNDARY
3091# ifdef SOLVE3D
3092 & nl_t_obc, nl_u_obc, nl_v_obc, &
3093# endif
3094 & nl_ubar_obc, nl_vbar_obc, &
3095 & nl_zeta_obc, &
3096# endif
3097# ifdef ADJUST_WSTRESS
3098 & nl_ustr, nl_vstr, &
3099# endif
3100# ifdef SOLVE3D
3101# ifdef ADJUST_STFLUX
3102 & nl_tflux, &
3103# endif
3104 & nl_t, nl_u, nl_v, &
3105# else
3106 & nl_ubar, nl_vbar, &
3107# endif
3108 & nl_zeta)
3109!
3110 ncname=hss(ng)%name
3111 DO ivec=1,ninner
3112 rec=ivec
3113!
3114! Read vectors stored in hessian NetCDF file using nl_var(Ltmp1) as
3115! temporary storage.
3116!
3117 CALL state_read (ng, tile, model, hss(ng)%IOtype, &
3118 & lbi, ubi, lbj, ubj, lbij, ubij, &
3119 & ltmp1, rec, &
3120 & 0, hss(ng)%ncid, ncname, &
3121# ifdef MASKING
3122 & rmask, umask, vmask, &
3123# endif
3124# ifdef ADJUST_BOUNDARY
3125# ifdef SOLVE3D
3126 & nl_t_obc, nl_u_obc, nl_v_obc, &
3127# endif
3128 & nl_ubar_obc, nl_vbar_obc, &
3129 & nl_zeta_obc, &
3130# endif
3131# ifdef ADJUST_WSTRESS
3132 & nl_ustr, nl_vstr, &
3133# endif
3134# ifdef SOLVE3D
3135# ifdef ADJUST_STFLUX
3136 & nl_tflux, &
3137# endif
3138 & nl_t, nl_u, nl_v, &
3139# else
3140 & nl_ubar, nl_vbar, &
3141# endif
3142 & nl_zeta)
3143!
3144! Compute the sum: (See NOTE above)
3145!
3146! nl_var(Ltmp2) = fac1 * nl_var(Ltmp2) + fac2 * nl_var(Ltmp1)
3147!
3148 fac1=1.0_r8
3149 fac2=zu(ivec)
3150!
3151 CALL state_addition (ng, tile, &
3152 & lbi, ubi, lbj, ubj, lbij, ubij, &
3153 & ltmp2, ltmp1, ltmp2, fac1, fac2, &
3154# ifdef MASKING
3155 & rmask, umask, vmask, &
3156# endif
3157# ifdef ADJUST_BOUNDARY
3158# ifdef SOLVE3D
3159 & nl_t_obc, nl_t_obc, &
3160 & nl_u_obc, nl_u_obc, &
3161 & nl_v_obc, nl_v_obc, &
3162# endif
3163 & nl_ubar_obc, nl_ubar_obc, &
3164 & nl_vbar_obc, nl_vbar_obc, &
3165 & nl_zeta_obc, nl_zeta_obc, &
3166# endif
3167# ifdef ADJUST_WSTRESS
3168 & nl_ustr, nl_ustr, &
3169 & nl_vstr, nl_vstr, &
3170# endif
3171# ifdef SOLVE3D
3172# ifdef ADJUST_STFLUX
3173 & nl_tflux, nl_tflux, &
3174# endif
3175 & nl_t, nl_t, &
3176 & nl_u, nl_u, &
3177 & nl_v, nl_v, &
3178# if defined WEAK_CONSTRAINT && defined TIME_CONV
3179 & nl_ubar, nl_ubar, &
3180 & nl_vbar, nl_vbar, &
3181# endif
3182# else
3183 & nl_ubar, nl_ubar, &
3184 & nl_vbar, nl_vbar, &
3185# endif
3186 & nl_zeta, nl_zeta)
3187 END DO
3188!
3189! Copy nl_var(Ltmp2) into ad_var(Lnew). See NOTE above.
3190!
3191 CALL state_copy (ng, tile, &
3192 & lbi, ubi, lbj, ubj, lbij, ubij, &
3193 & ltmp2, lnew, &
3194# ifdef ADJUST_BOUNDARY
3195# ifdef SOLVE3D
3196 & ad_t_obc, nl_t_obc, &
3197 & ad_u_obc, nl_u_obc, &
3198 & ad_v_obc, nl_v_obc, &
3199# endif
3200 & ad_ubar_obc, nl_ubar_obc, &
3201 & ad_vbar_obc, nl_vbar_obc, &
3202 & ad_zeta_obc, nl_zeta_obc, &
3203# endif
3204# ifdef ADJUST_WSTRESS
3205 & ad_ustr, nl_ustr, &
3206 & ad_vstr, nl_vstr, &
3207# endif
3208# ifdef SOLVE3D
3209# ifdef ADJUST_STFLUX
3210 & ad_tflux, nl_tflux, &
3211# endif
3212 & ad_t, nl_t, &
3213 & ad_u, nl_u, &
3214 & ad_v, nl_v, &
3215# if defined WEAK_CONSTRAINT && defined TIME_CONV
3216 & ad_ubar, nl_ubar, &
3217 & ad_vbar, nl_vbar, &
3218# endif
3219# else
3220 & ad_ubar, nl_ubar, &
3221 & ad_vbar, nl_vbar, &
3222# endif
3223 & ad_zeta, nl_zeta)
3224!
3225! Now form the final analysis error covariance matrix-vector product:
3226!
3227! y=(D-DG'VT^-1V'GD)x
3228!
3229! where tl_var(Lnew)=Dx and ad_var(Lnew)=DG'VT^-1V'GDx
3230!
3231! Free-surface.
3232!
3233 DO j=jstrt,jendt
3234 DO i=istrt,iendt
3235 ad_zeta(i,j,lnew)=tl_zeta(i,j,lnew)-ad_zeta(i,j,lnew)
3236# ifdef MASKING
3237 ad_zeta(i,j,lnew)=ad_zeta(i,j,lnew)*rmask(i,j)
3238# endif
3239 END DO
3240 END DO
3241
3242# ifdef ADJUST_BOUNDARY
3243!
3244! Free-surface open boundaries.
3245!
3246 IF (any(lobc(:,isfsur,ng))) THEN
3247 DO ir=1,nbrec(ng)
3248 IF ((lobc(iwest,isfsur,ng)).and. &
3249 & domain(ng)%Western_Edge(tile)) THEN
3250 ib=iwest
3251 DO j=jstr,jend
3252 ad_zeta_obc(j,ib,ir,lnew)=tl_zeta_obc(j,ib,ir,lnew)- &
3253 & ad_zeta_obc(j,ib,ir,lnew)
3254# ifdef MASKING
3255 ad_zeta_obc(j,ib,ir,lnew)=ad_zeta_obc(j,ib,ir,lnew)* &
3256 & rmask(istr-1,j)
3257# endif
3258 END DO
3259 END IF
3260 IF ((lobc(ieast,isfsur,ng)).and. &
3261 & domain(ng)%Eastern_Edge(tile)) THEN
3262 ib=ieast
3263 DO j=jstr,jend
3264 ad_zeta_obc(j,ib,ir,lnew)=tl_zeta_obc(j,ib,ir,lnew)- &
3265 & ad_zeta_obc(j,ib,ir,lnew)
3266# ifdef MASKING
3267 ad_zeta_obc(j,ib,ir,lnew)=ad_zeta_obc(j,ib,ir,lnew)* &
3268 & rmask(iend+1,j)
3269# endif
3270 END DO
3271 END IF
3272 IF ((lobc(isouth,isfsur,ng)).and. &
3273 & domain(ng)%Southern_Edge(tile)) THEN
3274 ib=isouth
3275 DO i=istr,iend
3276 ad_zeta_obc(i,ib,ir,lnew)=tl_zeta_obc(i,ib,ir,lnew)- &
3277 & ad_zeta_obc(i,ib,ir,lnew)
3278# ifdef MASKING
3279 ad_zeta_obc(i,ib,ir,lnew)=ad_zeta_obc(i,ib,ir,lnew)* &
3280 & rmask(i,jstr-1)
3281# endif
3282 END DO
3283 END IF
3284 IF ((lobc(inorth,isfsur,ng)).and. &
3285 & domain(ng)%Northern_Edge(tile)) THEN
3286 ib=inorth
3287 DO i=istr,iend
3288 ad_zeta_obc(i,ib,ir,lnew)=tl_zeta_obc(i,ib,ir,lnew)- &
3289 & ad_zeta_obc(i,ib,ir,lnew)
3290# ifdef MASKING
3291 ad_zeta_obc(i,ib,ir,lnew)=ad_zeta_obc(i,ib,ir,lnew)* &
3292 & rmask(i,jend+1)
3293# endif
3294 END DO
3295 END IF
3296 END DO
3297 END IF
3298# endif
3299
3300# ifndef SOLVE3D
3301!
3302! 2D U-momentum.
3303!
3304 DO j=jstrt,jendt
3305 DO i=istrp,iendt
3306 ad_ubar(i,j,lnew)=tl_ubar(i,j,lnew)-ad_ubar(i,j,lnew)
3307# ifdef MASKING
3308 ad_ubar(i,j,lnew)=ad_ubar(i,j,lnew)*umask(i,j)
3309# endif
3310 END DO
3311 END DO
3312# endif
3313
3314# ifdef ADJUST_BOUNDARY
3315!
3316! 2D U-momentum open boundaries.
3317!
3318 IF (any(lobc(:,isubar,ng))) THEN
3319 DO ir=1,nbrec(ng)
3320 IF ((lobc(iwest,isubar,ng)).and. &
3321 & domain(ng)%Western_Edge(tile)) THEN
3322 ib=iwest
3323 DO j=jstr,jend
3324 ad_ubar_obc(j,ib,ir,lnew)=tl_ubar_obc(j,ib,ir,lnew)- &
3325 & ad_ubar_obc(j,ib,ir,lnew)
3326# ifdef MASKING
3327 ad_ubar_obc(j,ib,ir,lnew)=ad_ubar_obc(j,ib,ir,lnew)* &
3328 & umask(istr,j)
3329# endif
3330 END DO
3331 END IF
3332 IF ((lobc(ieast,isubar,ng)).and. &
3333 & domain(ng)%Eastern_Edge(tile)) THEN
3334 ib=ieast
3335 DO j=jstr,jend
3336 ad_ubar_obc(j,ib,ir,lnew)=tl_ubar_obc(j,ib,ir,lnew)- &
3337 & ad_ubar_obc(j,ib,ir,lnew)
3338# ifdef MASKING
3339 ad_ubar_obc(j,ib,ir,lnew)=ad_ubar_obc(j,ib,ir,lnew)* &
3340 & umask(iend+1,j)
3341# endif
3342 END DO
3343 END IF
3344 IF ((lobc(isouth,isubar,ng)).and. &
3345 & domain(ng)%Southern_Edge(tile)) THEN
3346 ib=isouth
3347 DO i=istru,iend
3348 ad_ubar_obc(i,ib,ir,lnew)=tl_ubar_obc(i,ib,ir,lnew)- &
3349 & ad_ubar_obc(i,ib,ir,lnew)
3350# ifdef MASKING
3351 ad_ubar_obc(i,ib,ir,lnew)=ad_ubar_obc(i,ib,ir,lnew)* &
3352 & umask(i,jstr-1)
3353# endif
3354 END DO
3355 END IF
3356 IF ((lobc(inorth,isubar,ng)).and. &
3357 & domain(ng)%Northern_Edge(tile)) THEN
3358 ib=inorth
3359 DO i=istru,iend
3360 ad_ubar_obc(i,ib,ir,lnew)=tl_ubar_obc(i,ib,ir,lnew)- &
3361 & ad_ubar_obc(i,ib,ir,lnew)
3362# ifdef MASKING
3363 ad_ubar_obc(i,ib,ir,lnew)=ad_ubar_obc(i,ib,ir,lnew)* &
3364 & umask(i,jend+1)
3365# endif
3366 END DO
3367 END IF
3368 END DO
3369 END IF
3370# endif
3371
3372# ifndef SOLVE3D
3373!
3374! 2D V-momentum.
3375!
3376 DO j=jstrp,jendt
3377 DO i=istrt,iendt
3378 ad_vbar(i,j,lnew)=tl_vbar(i,j,lnew)-ad_vbar(i,j,lnew)
3379# ifdef MASKING
3380 ad_vbar(i,j,lnew)=ad_vbar(i,j,lnew)*vmask(i,j)
3381# endif
3382 END DO
3383 END DO
3384# endif
3385
3386# ifdef ADJUST_BOUNDARY
3387!
3388! 2D V-momentum open boundaries.
3389!
3390 IF (any(lobc(:,isvbar,ng))) THEN
3391 DO ir=1,nbrec(ng)
3392 IF ((lobc(iwest,isvbar,ng)).and. &
3393 & domain(ng)%Western_Edge(tile)) THEN
3394 ib=iwest
3395 DO j=jstrv,jend
3396 ad_vbar_obc(j,ib,ir,lnew)=tl_vbar_obc(j,ib,ir,lnew)- &
3397 & ad_vbar_obc(j,ib,ir,lnew)
3398# ifdef MASKING
3399 ad_vbar_obc(j,ib,ir,lnew)=ad_vbar_obc(j,ib,ir,lnew)* &
3400 & vmask(istr-1,j)
3401# endif
3402 END DO
3403 END IF
3404 IF ((lobc(ieast,isvbar,ng)).and. &
3405 & domain(ng)%Eastern_Edge(tile)) THEN
3406 ib=ieast
3407 DO j=jstrv,jend
3408 ad_vbar_obc(j,ib,ir,lnew)=tl_vbar_obc(j,ib,ir,lnew)- &
3409 & ad_vbar_obc(j,ib,ir,lnew)
3410# ifdef MASKING
3411 ad_vbar_obc(j,ib,ir,lnew)=ad_vbar_obc(j,ib,ir,lnew)* &
3412 & vmask(iend+1,j)
3413# endif
3414 END DO
3415 END IF
3416 IF ((lobc(isouth,isvbar,ng)).and. &
3417 & domain(ng)%Southern_Edge(tile)) THEN
3418 ib=isouth
3419 DO i=istr,iend
3420 ad_vbar_obc(i,ib,ir,lnew)=tl_vbar_obc(i,ib,ir,lnew)- &
3421 & ad_vbar_obc(i,ib,ir,lnew)
3422# ifdef MASKING
3423 ad_vbar_obc(i,ib,ir,lnew)=ad_vbar_obc(i,ib,ir,lnew)* &
3424 & vmask(i,jstr)
3425# endif
3426 END DO
3427 END IF
3428 IF ((lobc(inorth,isvbar,ng)).and. &
3429 & domain(ng)%Northern_Edge(tile)) THEN
3430 ib=inorth
3431 DO i=istr,iend
3432 ad_vbar_obc(i,ib,ir,lnew)=tl_vbar_obc(i,ib,ir,lnew)- &
3433 & ad_vbar_obc(i,ib,ir,lnew)
3434# ifdef MASKING
3435 ad_vbar_obc(i,ib,ir,lnew)=ad_vbar_obc(i,ib,ir,lnew)* &
3436 & vmask(i,jend+1)
3437# endif
3438 END DO
3439 END IF
3440 END DO
3441 END IF
3442# endif
3443
3444# ifdef ADJUST_WSTRESS
3445!
3446! Surface momentum stress.
3447!
3448 DO ir=1,nfrec(ng)
3449 DO j=jstrt,jendt
3450 DO i=istrp,iendt
3451 ad_ustr(i,j,ir,lnew)=tl_ustr(i,j,ir,lnew)- &
3452 & ad_ustr(i,j,ir,lnew)
3453# ifdef MASKING
3454 ad_ustr(i,j,ir,lnew)=ad_ustr(i,j,ir,lnew)*umask(i,j)
3455# endif
3456 END DO
3457 END DO
3458 DO j=jstrp,jendt
3459 DO i=istrt,iendt
3460 ad_vstr(i,j,ir,lnew)=tl_vstr(i,j,ir,lnew)- &
3461 & ad_vstr(i,j,ir,lnew)
3462# ifdef MASKING
3463 ad_vstr(i,j,ir,lnew)=ad_vstr(i,j,ir,lnew)*vmask(i,j)
3464# endif
3465 END DO
3466 END DO
3467 END DO
3468# endif
3469
3470# ifdef SOLVE3D
3471!
3472! 3D U-momentum.
3473!
3474 DO k=1,n(ng)
3475 DO j=jstrt,jendt
3476 DO i=istrp,iendt
3477 ad_u(i,j,k,lnew)=tl_u(i,j,k,lnew)-ad_u(i,j,k,lnew)
3478# ifdef MASKING
3479 ad_u(i,j,k,lnew)=ad_u(i,j,k,lnew)*umask(i,j)
3480# endif
3481 END DO
3482 END DO
3483 END DO
3484
3485# ifdef ADJUST_BOUNDARY
3486!
3487! 3D U-momentum open boundaries.
3488!
3489 IF (any(lobc(:,isuvel,ng))) THEN
3490 DO ir=1,nbrec(ng)
3491 IF ((lobc(iwest,isuvel,ng)).and. &
3492 & domain(ng)%Western_Edge(tile)) THEN
3493 ib=iwest
3494 DO k=1,n(ng)
3495 DO j=jstr,jend
3496 ad_u_obc(j,k,ib,ir,lnew)=tl_u_obc(j,k,ib,ir,lnew)- &
3497 & ad_u_obc(j,k,ib,ir,lnew)
3498# ifdef MASKING
3499 ad_u_obc(j,k,ib,ir,lnew)=ad_u_obc(j,k,ib,ir,lnew)* &
3500 & umask(istr,j)
3501# endif
3502 END DO
3503 END DO
3504 END IF
3505 IF ((lobc(ieast,isuvel,ng)).and. &
3506 & domain(ng)%Eastern_Edge(tile)) THEN
3507 ib=ieast
3508 DO k=1,n(ng)
3509 DO j=jstr,jend
3510 ad_u_obc(j,k,ib,ir,lnew)=tl_u_obc(j,k,ib,ir,lnew)- &
3511 & ad_u_obc(j,k,ib,ir,lnew)
3512# ifdef MASKING
3513 ad_u_obc(j,k,ib,ir,lnew)=ad_u_obc(j,k,ib,ir,lnew)* &
3514 & umask(iend+1,j)
3515# endif
3516 END DO
3517 END DO
3518 END IF
3519 IF ((lobc(isouth,isuvel,ng)).and. &
3520 & domain(ng)%Southern_Edge(tile)) THEN
3521 ib=isouth
3522 DO k=1,n(ng)
3523 DO i=istru,iend
3524 ad_u_obc(i,k,ib,ir,lnew)=tl_u_obc(i,k,ib,ir,lnew)- &
3525 & ad_u_obc(i,k,ib,ir,lnew)
3526# ifdef MASKING
3527 ad_u_obc(i,k,ib,ir,lnew)=ad_u_obc(i,k,ib,ir,lnew)* &
3528 & umask(i,jstr-1)
3529# endif
3530 END DO
3531 END DO
3532 END IF
3533 IF ((lobc(inorth,isuvel,ng)).and. &
3534 & domain(ng)%Northern_Edge(tile)) THEN
3535 ib=inorth
3536 DO k=1,n(ng)
3537 DO i=istru,iend
3538 ad_u_obc(i,k,ib,ir,lnew)=tl_u_obc(i,k,ib,ir,lnew)- &
3539 & ad_u_obc(i,k,ib,ir,lnew)
3540# ifdef MASKING
3541 ad_u_obc(i,k,ib,ir,lnew)=ad_u_obc(i,k,ib,ir,lnew)* &
3542 & umask(i,jend+1)
3543# endif
3544 END DO
3545 END DO
3546 END IF
3547 END DO
3548 END IF
3549# endif
3550!
3551! 3D V-momentum.
3552!
3553 DO k=1,n(ng)
3554 DO j=jstrp,jendt
3555 DO i=istrt,iendt
3556 ad_v(i,j,k,lnew)=tl_v(i,j,k,lnew)-ad_v(i,j,k,lnew)
3557# ifdef MASKING
3558 ad_v(i,j,k,lnew)=ad_v(i,j,k,lnew)*vmask(i,j)
3559# endif
3560 END DO
3561 END DO
3562 END DO
3563
3564# ifdef ADJUST_BOUNDARY
3565!
3566! 3D V-momentum open boundaries.
3567!
3568 IF (any(lobc(:,isvvel,ng))) THEN
3569 DO ir=1,nbrec(ng)
3570 IF ((lobc(iwest,isvvel,ng)).and. &
3571 & domain(ng)%Western_Edge(tile)) THEN
3572 ib=iwest
3573 DO k=1,n(ng)
3574 DO j=jstrv,jend
3575 ad_v_obc(j,k,ib,ir,lnew)=tl_v_obc(j,k,ib,ir,lnew)- &
3576 & ad_v_obc(j,k,ib,ir,lnew)
3577# ifdef MASKING
3578 ad_v_obc(j,k,ib,ir,lnew)=ad_v_obc(j,k,ib,ir,lnew)* &
3579 & vmask(istr-1,j)
3580# endif
3581 END DO
3582 END DO
3583 END IF
3584 IF ((lobc(ieast,isvvel,ng)).and. &
3585 & domain(ng)%Eastern_Edge(tile)) THEN
3586 ib=ieast
3587 DO k=1,n(ng)
3588 DO j=jstrv,jend
3589 ad_v_obc(j,k,ib,ir,lnew)=tl_v_obc(j,k,ib,ir,lnew)- &
3590 & ad_v_obc(j,k,ib,ir,lnew)
3591# ifdef MASKING
3592 ad_v_obc(j,k,ib,ir,lnew)=ad_v_obc(j,k,ib,ir,lnew)* &
3593 & vmask(iend+1,j)
3594# endif
3595 END DO
3596 END DO
3597 END IF
3598 IF ((lobc(isouth,isvvel,ng)).and. &
3599 & domain(ng)%Southern_Edge(tile)) THEN
3600 ib=isouth
3601 DO k=1,n(ng)
3602 DO i=istr,iend
3603 ad_v_obc(i,k,ib,ir,lnew)=tl_v_obc(i,k,ib,ir,lnew)- &
3604 & ad_v_obc(i,k,ib,ir,lnew)
3605# ifdef MASKING
3606 ad_v_obc(i,k,ib,ir,lnew)=ad_v_obc(i,k,ib,ir,lnew)* &
3607 & vmask(i,jstr)
3608# endif
3609 END DO
3610 END DO
3611 END IF
3612 IF ((lobc(inorth,isvvel,ng)).and. &
3613 & domain(ng)%Northern_Edge(tile)) THEN
3614 ib=inorth
3615 DO k=1,n(ng)
3616 DO i=istr,iend
3617 ad_v_obc(i,k,ib,ir,lnew)=tl_v_obc(i,k,ib,ir,lnew)- &
3618 & ad_v_obc(i,k,ib,ir,lnew)
3619# ifdef MASKING
3620 ad_v_obc(i,k,ib,ir,lnew)=ad_v_obc(i,k,ib,ir,lnew)* &
3621 & vmask(i,jend+1)
3622# endif
3623 END DO
3624 END DO
3625 END IF
3626 END DO
3627 END IF
3628# endif
3629!
3630! Tracers.
3631!
3632 DO it=1,nt(ng)
3633 DO k=1,n(ng)
3634 DO j=jstrt,jendt
3635 DO i=istrt,iendt
3636 ad_t(i,j,k,lnew,it)=tl_t(i,j,k,lnew,it)- &
3637 & ad_t(i,j,k,lnew,it)
3638# ifdef MASKING
3639 ad_t(i,j,k,lnew,it)=ad_t(i,j,k,lnew,it)*rmask(i,j)
3640# endif
3641 END DO
3642 END DO
3643 END DO
3644 END DO
3645
3646# ifdef ADJUST_BOUNDARY
3647!
3648! Tracers open boundaries.
3649!
3650 DO it=1,nt(ng)
3651 IF (any(lobc(:,istvar(it),ng))) THEN
3652 DO ir=1,nbrec(ng)
3653 IF ((lobc(iwest,istvar(it),ng)).and. &
3654 & domain(ng)%Western_Edge(tile)) THEN
3655 ib=iwest
3656 DO k=1,n(ng)
3657 DO j=jstr,jend
3658 ad_t_obc(j,k,ib,ir,lnew,it)= &
3659 & tl_t_obc(j,k,ib,ir,lnew,it)- &
3660 & ad_t_obc(j,k,ib,ir,lnew,it)
3661# ifdef MASKING
3662 ad_t_obc(j,k,ib,ir,lnew,it)= &
3663 & ad_t_obc(j,k,ib,ir,lnew,it)*rmask(istr-1,j)
3664# endif
3665 END DO
3666 END DO
3667 END IF
3668 IF ((lobc(ieast,istvar(it),ng)).and. &
3669 & domain(ng)%Eastern_Edge(tile)) THEN
3670 ib=ieast
3671 DO k=1,n(ng)
3672 DO j=jstr,jend
3673 ad_t_obc(j,k,ib,ir,lnew,it)= &
3674 & tl_t_obc(j,k,ib,ir,lnew,it)- &
3675 & ad_t_obc(j,k,ib,ir,lnew,it)
3676# ifdef MASKING
3677 ad_t_obc(j,k,ib,ir,lnew,it)= &
3678 & ad_t_obc(j,k,ib,ir,lnew,it)*rmask(iend+1,j)
3679# endif
3680 END DO
3681 END DO
3682 END IF
3683 IF ((lobc(isouth,istvar(it),ng)).and. &
3684 & domain(ng)%Southern_Edge(tile)) THEN
3685 ib=isouth
3686 DO k=1,n(ng)
3687 DO i=istr,iend
3688 ad_t_obc(i,k,ib,ir,lnew,it)= &
3689 & tl_t_obc(i,k,ib,ir,lnew,it)- &
3690 & ad_t_obc(i,k,ib,ir,lnew,it)
3691# ifdef MASKING
3692 ad_t_obc(i,k,ib,ir,lnew,it)= &
3693 & ad_t_obc(i,k,ib,ir,lnew,it)*rmask(i,jstr-1)
3694# endif
3695 END DO
3696 END DO
3697 END IF
3698 IF ((lobc(inorth,istvar(it),ng)).and. &
3699 & domain(ng)%Northern_Edge(tile)) THEN
3700 ib=inorth
3701 DO k=1,n(ng)
3702 DO i=istr,iend
3703 ad_t_obc(i,k,ib,ir,lnew,it)= &
3704 & tl_t_obc(i,k,ib,ir,lnew,it)- &
3705 & ad_t_obc(i,k,ib,ir,lnew,it)
3706# ifdef MASKING
3707 ad_t_obc(i,k,ib,ir,lnew,it)= &
3708 & ad_t_obc(i,k,ib,ir,lnew,it)*rmask(i,jend+1)
3709# endif
3710 END DO
3711 END DO
3712 END IF
3713 END DO
3714 END IF
3715 END DO
3716# endif
3717
3718# ifdef ADJUST_STFLUX
3719!
3720! Surface tracers flux.
3721!
3722 DO it=1,nt(ng)
3723 IF (lstflux(it,ng)) THEN
3724 DO ir=1,nfrec(ng)
3725 DO j=jstrt,jendt
3726 DO i=istrt,iendt
3727 ad_tflux(i,j,ir,lnew,it)=tl_tflux(i,j,ir,lnew,it)- &
3728 & ad_tflux(i,j,ir,lnew,it)
3729# ifdef MASKING
3730 ad_tflux(i,j,ir,lnew,it)=ad_tflux(i,j,ir,lnew,it)* &
3731 & rmask(i,j)
3732# endif
3733 END DO
3734 END DO
3735 END DO
3736 END IF
3737 END DO
3738# endif
3739# endif
3740!
3741!-----------------------------------------------------------------------
3742! Compute norm Delta(k) as the dot-product between the new vector
3743! and previous Lanczos vector.
3744!-----------------------------------------------------------------------
3745!
3746 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3747!
3748! Compute current iteration norm Delta(k) used to compute tri-diagonal
3749! matrix T(k) in the Lanczos recurrence.
3750!
3751 CALL state_dotprod (ng, tile, model, &
3752 & lbi, ubi, lbj, ubj, lbij, ubij, &
3753 & nstatevar(ng), dot(0:), &
3754# ifdef MASKING
3755 & rmask, umask, vmask, &
3756# endif
3757# ifdef ADJUST_BOUNDARY
3758# ifdef SOLVE3D
3759 & ad_t_obc(:,:,:,:,lnew,:), &
3760 & tl_t_obc(:,:,:,:,lwrk,:), &
3761 & ad_u_obc(:,:,:,:,lnew), &
3762 & tl_u_obc(:,:,:,:,lwrk), &
3763 & ad_v_obc(:,:,:,:,lnew), &
3764 & tl_v_obc(:,:,:,:,lwrk), &
3765# endif
3766 & ad_ubar_obc(:,:,:,lnew), &
3767 & tl_ubar_obc(:,:,:,lwrk), &
3768 & ad_vbar_obc(:,:,:,lnew), &
3769 & tl_vbar_obc(:,:,:,lwrk), &
3770 & ad_zeta_obc(:,:,:,lnew), &
3771 & tl_zeta_obc(:,:,:,lwrk), &
3772# endif
3773# ifdef ADJUST_WSTRESS
3774 & ad_ustr(:,:,:,lnew), tl_ustr(:,:,:,lwrk), &
3775 & ad_vstr(:,:,:,lnew), tl_vstr(:,:,:,lwrk), &
3776# endif
3777# ifdef SOLVE3D
3778# ifdef ADJUST_STFLUX
3779 & ad_tflux(:,:,:,lnew,:), &
3780 & tl_tflux(:,:,:,lwrk,:), &
3781# endif
3782 & ad_t(:,:,:,lnew,:), tl_t(:,:,:,lwrk,:), &
3783 & ad_u(:,:,:,lnew), tl_u(:,:,:,lwrk), &
3784 & ad_v(:,:,:,lnew), tl_v(:,:,:,lwrk), &
3785# else
3786 & ad_ubar(:,:,lnew), tl_ubar(:,:,lwrk), &
3787 & ad_vbar(:,:,lnew), tl_vbar(:,:,lwrk), &
3788# endif
3789 & ad_zeta(:,:,lnew), tl_zeta(:,:,lwrk))
3790
3791 ae_delta(innloop,outloop)=dot(0)
3792
3793 RETURN
3794 END SUBROUTINE analysis_error
3795!
3796!***********************************************************************
3797 SUBROUTINE lanczos (ng, tile, model, &
3798 & LBi, UBi, LBj, UBj, LBij, UBij, &
3799 & IminS, ImaxS, JminS, JmaxS, &
3800 & Lold, Lnew, Lwrk, &
3801 & innLoop, outLoop, &
3802# ifdef MASKING
3803 & rmask, umask, vmask, &
3804# endif
3805# ifdef ADJUST_BOUNDARY
3806# ifdef SOLVE3D
3807 & tl_t_obc, tl_u_obc, tl_v_obc, &
3808# endif
3809 & tl_ubar_obc, tl_vbar_obc, &
3810 & tl_zeta_obc, &
3811# endif
3812# ifdef ADJUST_WSTRESS
3813 & tl_ustr, tl_vstr, &
3814# endif
3815# ifdef SOLVE3D
3816# ifdef ADJUST_STFLUX
3817 & tl_tflux, &
3818# endif
3819 & tl_t, tl_u, tl_v, &
3820# if defined WEAK_CONSTRAINT && defined TIME_CONV
3821 & tl_ubar, tl_vbar, &
3822# endif
3823# else
3824 & tl_ubar, tl_vbar, &
3825# endif
3826 & tl_zeta, &
3827# ifdef ADJUST_BOUNDARY
3828# ifdef SOLVE3D
3829 & ad_t_obc, ad_u_obc, ad_v_obc, &
3830# endif
3831 & ad_ubar_obc, ad_vbar_obc, &
3832 & ad_zeta_obc, &
3833# endif
3834# ifdef ADJUST_WSTRESS
3835 & ad_ustr, ad_vstr, &
3836# endif
3837# ifdef SOLVE3D
3838# ifdef ADJUST_STFLUX
3839 & ad_tflux, &
3840# endif
3841 & ad_t, ad_u, ad_v, &
3842# if defined WEAK_CONSTRAINT && defined TIME_CONV
3843 & ad_ubar, ad_vbar, &
3844# endif
3845# else
3846 & ad_ubar, ad_vbar, &
3847# endif
3848 & ad_zeta)
3849!***********************************************************************
3850!
3851! Imported variable declarations.
3852!
3853 integer, intent(in) :: ng, tile, model
3854 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
3855 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
3856 integer, intent(in) :: Lold, Lnew, Lwrk
3857 integer, intent(in) :: innLoop, outLoop
3858!
3859# ifdef ASSUMED_SHAPE
3860# ifdef MASKING
3861 real(r8), intent(in) :: rmask(LBi:,LBj:)
3862 real(r8), intent(in) :: umask(LBi:,LBj:)
3863 real(r8), intent(in) :: vmask(LBi:,LBj:)
3864# endif
3865# ifdef ADJUST_BOUNDARY
3866# ifdef SOLVE3D
3867 real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
3868 real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
3869 real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
3870# endif
3871 real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
3872 real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
3873 real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
3874# endif
3875# ifdef ADJUST_WSTRESS
3876 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
3877 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
3878# endif
3879# ifdef SOLVE3D
3880# ifdef ADJUST_STFLUX
3881 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
3882# endif
3883 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
3884 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
3885 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
3886# if defined WEAK_CONSTRAINT && defined TIME_CONV
3887 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
3888 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
3889# endif
3890# else
3891 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
3892 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
3893# endif
3894 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
3895# ifdef ADJUST_BOUNDARY
3896# ifdef SOLVE3D
3897 real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
3898 real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
3899 real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
3900# endif
3901 real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
3902 real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
3903 real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
3904# endif
3905# ifdef ADJUST_WSTRESS
3906 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
3907 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
3908# endif
3909# ifdef SOLVE3D
3910# ifdef ADJUST_STFLUX
3911 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
3912# endif
3913 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
3914 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
3915 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
3916# if defined WEAK_CONSTRAINT && defined TIME_CONV
3917 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
3918 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
3919# endif
3920# else
3921 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
3922 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
3923# endif
3924 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
3925
3926# else
3927
3928# ifdef MASKING
3929 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
3930 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
3931 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
3932# endif
3933# ifdef ADJUST_BOUNDARY
3934# ifdef SOLVE3D
3935 real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, &
3936 & Nbrec(ng),2,NT(ng))
3937 real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
3938 real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
3939# endif
3940 real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
3941 real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
3942 real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
3943# endif
3944# ifdef ADJUST_WSTRESS
3945 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
3946 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
3947# endif
3948# ifdef SOLVE3D
3949# ifdef ADJUST_STFLUX
3950 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, &
3951 & Nfrec(ng),2,NT(ng))
3952# endif
3953 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
3954 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
3955 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
3956# if defined WEAK_CONSTRAINT && defined TIME_CONV
3957 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
3958 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
3959# endif
3960# else
3961 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
3962 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
3963# endif
3964 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
3965# ifdef ADJUST_BOUNDARY
3966# ifdef SOLVE3D
3967 real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, &
3968 & Nbrec(ng),2,NT(ng))
3969 real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
3970 real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
3971# endif
3972 real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
3973 real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
3974 real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
3975# endif
3976# ifdef ADJUST_WSTRESS
3977 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
3978 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
3979# endif
3980# ifdef SOLVE3D
3981# ifdef ADJUST_STFLUX
3982 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, &
3983 & Nfrec(ng),2,NT(ng))
3984# endif
3985 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
3986 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
3987 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
3988# if defined WEAK_CONSTRAINT && defined TIME_CONV
3989 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
3990 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
3991# endif
3992# else
3993 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
3994 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
3995# endif
3996 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
3997# endif
3998!
3999! Local variable declarations.
4000!
4001 integer :: i, j, lstr, rec
4002!
4003 real(r8) :: fac, fac1, fac2
4004!
4005 real(r8), dimension(0:NstateVar(ng)) :: dot
4006 real(r8), dimension(0:NpostI) :: DotProd, dot_new, dot_old
4007!
4008 character (len=256) :: ncname
4009
4010 character (len=*), parameter :: MyFile = &
4011 & __FILE__//", lanczos"
4012
4013# include "set_bounds.h"
4014!
4015!-----------------------------------------------------------------------
4016! Calculate the new Lanczos vector, q(k+1) using reccurence equation
4017! for the gradient vectors:
4018!
4019! H q(k+1) = Gamma(k+1) q(k+2) + Delta(k+1) q(k+1) + Gamma(k) q(k)
4020!
4021! where Gamma(k) = - SQRT ( Beta(k+1) ) / Alpha(k)
4022!-----------------------------------------------------------------------
4023!
4024! NOTE: In the case of weak constraint ("WEAK_CONSTRAINT") and
4025! and time convolutions ("TIME_CONV"), the state arrays
4026! tl_ubar, tl_vbar, ad_ubar, and ad_vbar are only passed
4027! as required by the "state" operators routines but they
4028! are not used in subsequent calculations.
4029!
4030! At this point, the previous orthonormal Lanczos vector is still in
4031! tangent linear state arrays (index Lwrk) - it was read in the
4032! routine anerror.
4033!
4034 IF (innloop.gt.0) THEN
4035!
4036! Compute new Lanczos vector: (See NOTE above)
4037!
4038! ad_var(Lnew) = fac1 * ad_var(Lnew) + fac2 * tl_var(Lwrk)
4039!
4040 fac1=1.0_r8
4041 fac2=-ae_delta(innloop,outloop)
4042!
4043 CALL state_addition (ng, tile, &
4044 & lbi, ubi, lbj, ubj, lbij, ubij, &
4045 & lnew, lwrk, lnew, fac1, fac2, &
4046# ifdef MASKING
4047 & rmask, umask, vmask, &
4048# endif
4049# ifdef ADJUST_BOUNDARY
4050# ifdef SOLVE3D
4051 & ad_t_obc, tl_t_obc, &
4052 & ad_u_obc, tl_u_obc, &
4053 & ad_v_obc, tl_v_obc, &
4054# endif
4055 & ad_ubar_obc, tl_ubar_obc, &
4056 & ad_vbar_obc, tl_vbar_obc, &
4057 & ad_zeta_obc, tl_zeta_obc, &
4058# endif
4059# ifdef ADJUST_WSTRESS
4060 & ad_ustr, tl_ustr, &
4061 & ad_vstr, tl_vstr, &
4062# endif
4063# ifdef SOLVE3D
4064# ifdef ADJUST_STFLUX
4065 & ad_tflux, tl_tflux, &
4066# endif
4067 & ad_t, tl_t, &
4068 & ad_u, tl_u, &
4069 & ad_v, tl_v, &
4070# if defined WEAK_CONSTRAINT && defined TIME_CONV
4071 & ad_ubar, tl_ubar, &
4072 & ad_vbar, tl_vbar, &
4073# endif
4074# else
4075 & ad_ubar, tl_ubar, &
4076 & ad_vbar, tl_vbar, &
4077# endif
4078 & ad_zeta, tl_zeta)
4079 END IF
4080!
4081! Substract previous orthonormal Lanczos vector.
4082!
4083 IF (innloop.gt.1) THEN
4084!
4085! Determine adjoint file to process.
4086!
4087 ncname=adm(ng)%name
4088!
4089! Read in the previous (innLoop-1) orthonormal Lanczos vector.
4090!
4091 CALL state_read (ng, tile, model, adm(ng)%IOtype, &
4092 & lbi, ubi, lbj, ubj, lbij, ubij, &
4093 & lwrk, innloop-1, &
4094 & ndefadj(ng), adm(ng)%ncid, ncname, &
4095# ifdef MASKING
4096 & rmask, umask, vmask, &
4097# endif
4098# ifdef ADJUST_BOUNDARY
4099# ifdef SOLVE3D
4100 & tl_t_obc, tl_u_obc, tl_v_obc, &
4101# endif
4102 & tl_ubar_obc, tl_vbar_obc, &
4103 & tl_zeta_obc, &
4104# endif
4105# ifdef ADJUST_WSTRESS
4106 & tl_ustr, tl_vstr, &
4107# endif
4108# ifdef SOLVE3D
4109# ifdef ADJUST_STFLUX
4110 & tl_tflux, &
4111# endif
4112 & tl_t, tl_u, tl_v, &
4113# else
4114 & tl_ubar, tl_vbar, &
4115# endif
4116 & tl_zeta)
4117 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4118!
4119! Substract previous orthonormal Lanczos vector: (See NOTE above)
4120!
4121! ad_var(Lnew) = fac1 * ad_var(Lnew) + fac2 * tl_var(Lwrk)
4122!
4123 fac1=1.0_r8
4124 fac2=-ae_beta(innloop,outloop)
4125
4126 CALL state_addition (ng, tile, &
4127 & lbi, ubi, lbj, ubj, lbij, ubij, &
4128 & lnew, lwrk, lnew, fac1, fac2, &
4129# ifdef MASKING
4130 & rmask, umask, vmask, &
4131# endif
4132# ifdef ADJUST_BOUNDARY
4133# ifdef SOLVE3D
4134 & ad_t_obc, tl_t_obc, &
4135 & ad_u_obc, tl_u_obc, &
4136 & ad_v_obc, tl_v_obc, &
4137# endif
4138 & ad_ubar_obc, tl_ubar_obc, &
4139 & ad_vbar_obc, tl_vbar_obc, &
4140 & ad_zeta_obc, tl_zeta_obc, &
4141# endif
4142# ifdef ADJUST_WSTRESS
4143 & ad_ustr, tl_ustr, &
4144 & ad_vstr, tl_vstr, &
4145# endif
4146# ifdef SOLVE3D
4147# ifdef ADJUST_STFLUX
4148 & ad_tflux, tl_tflux, &
4149# endif
4150 & ad_t, tl_t, &
4151 & ad_u, tl_u, &
4152 & ad_v, tl_v, &
4153# if defined WEAK_CONSTRAINT && defined TIME_CONV
4154 & ad_ubar, tl_ubar, &
4155 & ad_vbar, tl_vbar, &
4156# endif
4157# else
4158 & ad_ubar, tl_ubar, &
4159 & ad_vbar, tl_vbar, &
4160# endif
4161 & ad_zeta, tl_zeta)
4162 END IF
4163!
4164!-----------------------------------------------------------------------
4165! Orthogonalize current gradient, q(k+1), against all previous
4166! gradients (reverse order) using Gramm-Schmidt procedure.
4167!-----------------------------------------------------------------------
4168!
4169! We can overwrite adjoint arrays at index Lnew each time around the
4170! the following loop because the preceding gradient vectors that we
4171! read are orthogonal to each other. The reversed order of the loop
4172! is important for the Lanczos vector calculations.
4173!
4174 ncname=adm(ng)%name
4175!
4176 DO rec=innloop,1,-1
4177!
4178! Read in each previous gradient state solutions, G(0) to G(k), and
4179! compute its associated dot angaint curret G(k+1). Each gradient
4180! solution is loaded into TANGENT LINEAR STATE ARRAYS at index Lwrk.
4181!
4182 CALL state_read (ng, tile, model, adm(ng)%IOtype, &
4183 & lbi, ubi, lbj, ubj, lbij, ubij, &
4184 & lwrk, rec, &
4185 & ndefadj(ng), adm(ng)%ncid, ncname, &
4186# ifdef MASKING
4187 & rmask, umask, vmask, &
4188# endif
4189# ifdef ADJUST_BOUNDARY
4190# ifdef SOLVE3D
4191 & tl_t_obc, tl_u_obc, tl_v_obc, &
4192# endif
4193 & tl_ubar_obc, tl_vbar_obc, &
4194 & tl_zeta_obc, &
4195# endif
4196# ifdef ADJUST_WSTRESS
4197 & tl_ustr, tl_vstr, &
4198# endif
4199# ifdef SOLVE3D
4200# ifdef ADJUST_STFLUX
4201 & tl_tflux, &
4202# endif
4203 & tl_t, tl_u, tl_v, &
4204# else
4205 & tl_ubar, tl_vbar, &
4206# endif
4207 & tl_zeta)
4208 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4209!
4210! Compute dot product <q(k+1), q(rec)>.
4211!
4212 CALL state_dotprod (ng, tile, model, &
4213 & lbi, ubi, lbj, ubj, lbij, ubij, &
4214 & nstatevar(ng), dot(0:), &
4215# ifdef MASKING
4216 & rmask, umask, vmask, &
4217# endif
4218# ifdef ADJUST_BOUNDARY
4219# ifdef SOLVE3D
4220 & ad_t_obc(:,:,:,:,lnew,:), &
4221 & tl_t_obc(:,:,:,:,lwrk,:), &
4222 & ad_u_obc(:,:,:,:,lnew), &
4223 & tl_u_obc(:,:,:,:,lwrk), &
4224 & ad_v_obc(:,:,:,:,lnew), &
4225 & tl_v_obc(:,:,:,:,lwrk), &
4226# endif
4227 & ad_ubar_obc(:,:,:,lnew), &
4228 & tl_ubar_obc(:,:,:,lwrk), &
4229 & ad_vbar_obc(:,:,:,lnew), &
4230 & tl_vbar_obc(:,:,:,lwrk), &
4231 & ad_zeta_obc(:,:,:,lnew), &
4232 & tl_zeta_obc(:,:,:,lwrk), &
4233# endif
4234# ifdef ADJUST_WSTRESS
4235 & ad_ustr(:,:,:,lnew), tl_ustr(:,:,:,lwrk), &
4236 & ad_vstr(:,:,:,lnew), tl_vstr(:,:,:,lwrk), &
4237# endif
4238# ifdef SOLVE3D
4239# ifdef ADJUST_STFLUX
4240 & ad_tflux(:,:,:,lnew,:), &
4241 & tl_tflux(:,:,:,lwrk,:), &
4242# endif
4243 & ad_t(:,:,:,lnew,:), tl_t(:,:,:,lwrk,:), &
4244 & ad_u(:,:,:,lnew), tl_u(:,:,:,lwrk), &
4245 & ad_v(:,:,:,lnew), tl_v(:,:,:,lwrk), &
4246# else
4247 & ad_ubar(:,:,lnew), tl_ubar(:,:,lwrk), &
4248 & ad_vbar(:,:,lnew), tl_vbar(:,:,lwrk), &
4249# endif
4250 & ad_zeta(:,:,lnew), tl_zeta(:,:,lwrk))
4251!
4252! Compute Gramm-Schmidt scaling coefficient.
4253!
4254 dotprod(rec)=dot(0)
4255!
4256! Gramm-Schmidt orthonormalization, free-surface. See NOTE above.
4257!
4258! ad_var(Lnew) = fac1 * ad_var(Lnew) + fac2 * tl_var(Lwrk)
4259!
4260 fac1=1.0_r8
4261 fac2=-dotprod(rec)
4262
4263 CALL state_addition (ng, tile, &
4264 & lbi, ubi, lbj, ubj, lbij, ubij, &
4265 & lnew, lwrk, lnew, fac1, fac2, &
4266# ifdef MASKING
4267 & rmask, umask, vmask, &
4268# endif
4269# ifdef ADJUST_BOUNDARY
4270# ifdef SOLVE3D
4271 & ad_t_obc, tl_t_obc, &
4272 & ad_u_obc, tl_u_obc, &
4273 & ad_v_obc, tl_v_obc, &
4274# endif
4275 & ad_ubar_obc, tl_ubar_obc, &
4276 & ad_vbar_obc, tl_vbar_obc, &
4277 & ad_zeta_obc, tl_zeta_obc, &
4278# endif
4279# ifdef ADJUST_WSTRESS
4280 & ad_ustr, tl_ustr, &
4281 & ad_vstr, tl_vstr, &
4282# endif
4283# ifdef SOLVE3D
4284# ifdef ADJUST_STFLUX
4285 & ad_tflux, tl_tflux, &
4286# endif
4287 & ad_t, tl_t, &
4288 & ad_u, tl_u, &
4289 & ad_v, tl_v, &
4290# if defined WEAK_CONSTRAINT && defined TIME_CONV
4291 & ad_ubar, tl_ubar, &
4292 & ad_vbar, tl_vbar, &
4293# endif
4294# else
4295 & ad_ubar, tl_ubar, &
4296 & ad_vbar, tl_vbar, &
4297# endif
4298 & ad_zeta, tl_zeta)
4299 END DO
4300!
4301!-----------------------------------------------------------------------
4302! Normalize current orthogonal Lanczos vector.
4303!-----------------------------------------------------------------------
4304!
4305 CALL state_dotprod (ng, tile, model, &
4306 & lbi, ubi, lbj, ubj, lbij, ubij, &
4307 & nstatevar(ng), dot(0:), &
4308# ifdef MASKING
4309 & rmask, umask, vmask, &
4310# endif
4311# ifdef ADJUST_BOUNDARY
4312# ifdef SOLVE3D
4313 & ad_t_obc(:,:,:,:,lnew,:), &
4314 & ad_t_obc(:,:,:,:,lnew,:), &
4315 & ad_u_obc(:,:,:,:,lnew), &
4316 & ad_u_obc(:,:,:,:,lnew), &
4317 & ad_v_obc(:,:,:,:,lnew), &
4318 & ad_v_obc(:,:,:,:,lnew), &
4319# endif
4320 & ad_ubar_obc(:,:,:,lnew), &
4321 & ad_ubar_obc(:,:,:,lnew), &
4322 & ad_vbar_obc(:,:,:,lnew), &
4323 & ad_vbar_obc(:,:,:,lnew), &
4324 & ad_zeta_obc(:,:,:,lnew), &
4325 & ad_zeta_obc(:,:,:,lnew), &
4326# endif
4327# ifdef ADJUST_WSTRESS
4328 & ad_ustr(:,:,:,lnew), ad_ustr(:,:,:,lnew), &
4329 & ad_vstr(:,:,:,lnew), ad_vstr(:,:,:,lnew), &
4330# endif
4331# ifdef SOLVE3D
4332# ifdef ADJUST_STFLUX
4333 & ad_tflux(:,:,:,lnew,:), &
4334 & ad_tflux(:,:,:,lnew,:), &
4335# endif
4336 & ad_t(:,:,:,lnew,:), ad_t(:,:,:,lnew,:), &
4337 & ad_u(:,:,:,lnew), ad_u(:,:,:,lnew), &
4338 & ad_v(:,:,:,lnew), ad_v(:,:,:,lnew), &
4339# else
4340 & ad_ubar(:,:,lnew), ad_ubar(:,:,lnew), &
4341 & ad_vbar(:,:,lnew), ad_vbar(:,:,lnew), &
4342# endif
4343 & ad_zeta(:,:,lnew), ad_zeta(:,:,lnew))
4344!
4345! Compute normalization factor.
4346!
4347 IF (innloop.eq.0) THEN
4348 ae_gnorm(outloop)=sqrt(dot(0))
4349 ELSE
4350 ae_beta(innloop+1,outloop)=sqrt(dot(0))
4351 END IF
4352!
4353! Normalize the vector: ad_var(Lnew) = fac * ad_var(Lnew)
4354!
4355 fac=1.0_r8/sqrt(dot(0))
4356
4357 CALL state_scale (ng, tile, &
4358 & lbi, ubi, lbj, ubj, lbij, ubij, &
4359 & lnew, lnew, fac, &
4360# ifdef MASKING
4361 & rmask, umask, vmask, &
4362# endif
4363# ifdef ADJUST_BOUNDARY
4364# ifdef SOLVE3D
4365 & ad_t_obc, ad_u_obc, ad_v_obc, &
4366# endif
4367 & ad_ubar_obc, ad_vbar_obc, &
4368 & ad_zeta_obc, &
4369# endif
4370# ifdef ADJUST_WSTRESS
4371 & ad_ustr, ad_vstr, &
4372# endif
4373# ifdef SOLVE3D
4374# ifdef ADJUST_STFLUX
4375 & ad_tflux, &
4376# endif
4377 & ad_t, ad_u, ad_v, &
4378# else
4379 & ad_ubar, ad_vbar, &
4380# endif
4381 & ad_zeta)
4382!
4383# ifdef TEST_ORTHOGONALIZATION
4384!
4385!-----------------------------------------------------------------------
4386! Test orthogonal properties of the new gradient.
4387!-----------------------------------------------------------------------
4388!
4389! Determine adjoint file to process.
4390!
4391 ncname=adm(ng)%name
4392!
4393 DO rec=innloop,1,-1
4394!
4395! Read in each previous gradient state solutions, q(0) to q(k), and
4396! compute its associated dot angaint orthogonalized q(k+1). Again,
4397! each gradient solution is loaded into TANGENT LINEAR STATE ARRAYS
4398! at index Lwrk.
4399!
4400 CALL state_read (ng, tile, model, adm(ng)%IOtype, &
4401 & lbi, ubi, lbj, ubj, lbij, ubij, &
4402 & lwrk, rec, &
4403 & ndefadj(ng), adm(ng)%ncid, ncname, &
4404# ifdef MASKING
4405 & rmask, umask, vmask, &
4406# endif
4407# ifdef ADJUST_BOUNDARY
4408# ifdef SOLVE3D
4409 & tl_t_obc, tl_u_obc, tl_v_obc, &
4410# endif
4411 & tl_ubar_obc, tl_vbar_obc, &
4412 & tl_zeta_obc, &
4413# endif
4414# ifdef ADJUST_WSTRESS
4415 & tl_ustr, tl_vstr, &
4416# endif
4417# ifdef SOLVE3D
4418# ifdef ADJUST_STFLUX
4419 & tl_tflux, &
4420# endif
4421 & tl_t, tl_u, tl_v, &
4422# else
4423 & tl_ubar, tl_vbar, &
4424# endif
4425 & tl_zeta)
4426 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4427!
4428 CALL state_dotprod (ng, tile, model, &
4429 & lbi, ubi, lbj, ubj, lbij, ubij, &
4430 & nstatevar(ng), dot(0:), &
4431# ifdef MASKING
4432 & rmask, umask, vmask, &
4433# endif
4434# ifdef ADJUST_BOUNDARY
4435# ifdef SOLVE3D
4436 & ad_t_obc(:,:,:,:,lnew,:), &
4437 & tl_t_obc(:,:,:,:,lwrk,:), &
4438 & ad_u_obc(:,:,:,:,lnew), &
4439 & tl_u_obc(:,:,:,:,lwrk), &
4440 & ad_v_obc(:,:,:,:,lnew), &
4441 & tl_v_obc(:,:,:,:,lwrk), &
4442# endif
4443 & ad_ubar_obc(:,:,:,lnew), &
4444 & tl_ubar_obc(:,:,:,lwrk), &
4445 & ad_vbar_obc(:,:,:,lnew), &
4446 & tl_vbar_obc(:,:,:,lwrk), &
4447 & ad_zeta_obc(:,:,:,lnew), &
4448 & tl_zeta_obc(:,:,:,lwrk), &
4449# endif
4450# ifdef ADJUST_WSTRESS
4451 & ad_ustr(:,:,:,lnew), tl_ustr(:,:,:,lwrk), &
4452 & ad_vstr(:,:,:,lnew), tl_vstr(:,:,:,lwrk), &
4453# endif
4454# ifdef SOLVE3D
4455# ifdef ADJUST_STFLUX
4456 & ad_tflux(:,:,:,lnew,:), &
4457 & tl_tflux(:,:,:,lwrk,:), &
4458# endif
4459 & ad_t(:,:,:,lnew,:), tl_t(:,:,:,lwrk,:), &
4460 & ad_u(:,:,:,lnew), tl_u(:,:,:,lwrk), &
4461 & ad_v(:,:,:,lnew), tl_v(:,:,:,lwrk), &
4462# else
4463 & ad_ubar(:,:,lnew), tl_ubar(:,:,lwrk), &
4464 & ad_vbar(:,:,lnew), tl_vbar(:,:,lwrk), &
4465# endif
4466 & ad_zeta(:,:,lnew), tl_zeta(:,:,lwrk))
4467 dot_new(rec)=dot(0)
4468 END DO
4469!
4470! Report dot products. If everything is working correctly, at the
4471! end of the orthogonalization dot_new(rec) << dot_old(rec).
4472!
4473 IF (master) THEN
4474 WRITE (stdout,20) outloop, innloop
4475 DO rec=innloop,1,-1
4476 WRITE (stdout,30) dotprod(rec), rec-1
4477 END DO
4478 WRITE (stdout,*) ' '
4479 DO rec=innloop,1,-1
4480 WRITE (stdout,40) innloop, rec-1, dot_new(rec), &
4481 & rec-1, rec-1, dot_old(rec)
4482 END DO
4483 20 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
4484 & 'Gramm-Schmidt Orthogonalization:',/)
4485 30 FORMAT (12x,'Orthogonalization Factor = ',1p,e19.12,3x, &
4486 & '(Iter=',i3.3,')')
4487 40 FORMAT (2x,'Ortho Test: ', &
4488 & '<G(',i3.3,'),G(',i3.3,')> = ',1p,e15.8,1x, &
4489 & '<G(',i3.3,'),G(',i3.3,')> = ',1p,e15.8)
4490 END IF
4491# endif
4492
4493 RETURN
4494 END SUBROUTINE lanczos
4495!
4496!***********************************************************************
4497 SUBROUTINE posterior_eofs (ng, tile, model, &
4498 & LBi, UBi, LBj, UBj, LBij, UBij, &
4499 & IminS, ImaxS, JminS, JmaxS, &
4500 & Lold, Lnew, Lwrk, &
4501 & innLoop, outLoop, &
4502# ifdef MASKING
4503 & rmask, umask, vmask, &
4504# endif
4505# ifdef ADJUST_BOUNDARY
4506# ifdef SOLVE3D
4507 & nl_t_obc, nl_u_obc, nl_v_obc, &
4508# endif
4509 & nl_ubar_obc, nl_vbar_obc, &
4510 & nl_zeta_obc, &
4511# endif
4512# ifdef ADJUST_WSTRESS
4513 & nl_ustr, nl_vstr, &
4514# endif
4515# ifdef SOLVE3D
4516# ifdef ADJUST_STFLUX
4517 & nl_tflux, &
4518# endif
4519 & nl_t, nl_u, nl_v, &
4520# if defined WEAK_CONSTRAINT && defined TIME_CONV
4521 & nl_ubar, nl_vbar, &
4522# endif
4523# else
4524 & nl_ubar, nl_vbar, &
4525# endif
4526 & nl_zeta, &
4527# ifdef ADJUST_BOUNDARY
4528# ifdef SOLVE3D
4529 & tl_t_obc, tl_u_obc, tl_v_obc, &
4530# endif
4531 & tl_ubar_obc, tl_vbar_obc, &
4532 & tl_zeta_obc, &
4533# endif
4534# ifdef ADJUST_WSTRESS
4535 & tl_ustr, tl_vstr, &
4536# endif
4537# ifdef SOLVE3D
4538# ifdef ADJUST_STFLUX
4539 & tl_tflux, &
4540# endif
4541 & tl_t, tl_u, tl_v, &
4542# if defined WEAK_CONSTRAINT && defined TIME_CONV
4543 & tl_ubar, tl_vbar, &
4544# endif
4545# else
4546 & tl_ubar, tl_vbar, &
4547# endif
4548 & tl_zeta, &
4549# ifdef ADJUST_BOUNDARY
4550# ifdef SOLVE3D
4551 & ad_t_obc, ad_u_obc, ad_v_obc, &
4552# endif
4553 & ad_ubar_obc, ad_vbar_obc, &
4554 & ad_zeta_obc, &
4555# endif
4556# ifdef ADJUST_WSTRESS
4557 & ad_ustr, ad_vstr, &
4558# endif
4559# ifdef SOLVE3D
4560# ifdef ADJUST_STFLUX
4561 & ad_tflux, &
4562# endif
4563 & ad_t, ad_u, ad_v, &
4564# if defined WEAK_CONSTRAINT && defined TIME_CONV
4565 & ad_ubar, ad_vbar, &
4566# endif
4567# else
4568 & ad_ubar, ad_vbar, &
4569# endif
4570 & ad_zeta)
4571!***********************************************************************
4572!
4573! Imported variable declarations.
4574!
4575 integer, intent(in) :: ng, tile, model
4576 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
4577 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
4578 integer, intent(in) :: Lold, Lnew, Lwrk
4579 integer, intent(in) :: innLoop, outLoop
4580!
4581# ifdef ASSUMED_SHAPE
4582# ifdef MASKING
4583 real(r8), intent(in) :: rmask(LBi:,LBj:)
4584 real(r8), intent(in) :: umask(LBi:,LBj:)
4585 real(r8), intent(in) :: vmask(LBi:,LBj:)
4586# endif
4587# ifdef ADJUST_BOUNDARY
4588# ifdef SOLVE3D
4589 real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
4590 real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
4591 real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
4592# endif
4593 real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
4594 real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
4595 real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
4596# endif
4597# ifdef ADJUST_WSTRESS
4598 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
4599 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
4600# endif
4601# ifdef SOLVE3D
4602# ifdef ADJUST_STFLUX
4603 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
4604# endif
4605 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
4606 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
4607 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
4608# else
4609 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
4610 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
4611# endif
4612 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
4613# ifdef ADJUST_BOUNDARY
4614# ifdef SOLVE3D
4615 real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
4616 real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
4617 real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
4618# endif
4619 real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
4620 real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
4621 real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
4622# endif
4623# ifdef ADJUST_WSTRESS
4624 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
4625 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
4626# endif
4627# ifdef SOLVE3D
4628# ifdef ADJUST_STFLUX
4629 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
4630# endif
4631 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
4632 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
4633 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
4634# else
4635 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
4636 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
4637# endif
4638 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
4639# ifdef ADJUST_BOUNDARY
4640# ifdef SOLVE3D
4641 real(r8), intent(inout) :: nl_t_obc(LBij:,:,:,:,:,:)
4642 real(r8), intent(inout) :: nl_u_obc(LBij:,:,:,:,:)
4643 real(r8), intent(inout) :: nl_v_obc(LBij:,:,:,:,:)
4644# endif
4645 real(r8), intent(inout) :: nl_ubar_obc(LBij:,:,:,:)
4646 real(r8), intent(inout) :: nl_vbar_obc(LBij:,:,:,:)
4647 real(r8), intent(inout) :: nl_zeta_obc(LBij:,:,:,:)
4648# endif
4649# ifdef ADJUST_WSTRESS
4650 real(r8), intent(inout) :: nl_ustr(LBi:,LBj:,:,:)
4651 real(r8), intent(inout) :: nl_vstr(LBi:,LBj:,:,:)
4652# endif
4653# ifdef SOLVE3D
4654# ifdef ADJUST_STFLUX
4655 real(r8), intent(inout) :: nl_tflux(LBi:,LBj:,:,:,:)
4656# endif
4657 real(r8), intent(inout) :: nl_t(LBi:,LBj:,:,:,:)
4658 real(r8), intent(inout) :: nl_u(LBi:,LBj:,:,:)
4659 real(r8), intent(inout) :: nl_v(LBi:,LBj:,:,:)
4660# else
4661 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
4662 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
4663# endif
4664 real(r8), intent(inout) :: nl_zeta(LBi:,LBj:,:)
4665
4666# else
4667
4668# ifdef MASKING
4669 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
4670 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
4671 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
4672# endif
4673# ifdef ADJUST_BOUNDARY
4674# ifdef SOLVE3D
4675 real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, &
4676 & Nbrec(ng),2,NT(ng))
4677 real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
4678 real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
4679# endif
4680 real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
4681 real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
4682 real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
4683# endif
4684# ifdef ADJUST_WSTRESS
4685 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
4686 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
4687# endif
4688# ifdef SOLVE3D
4689# ifdef ADJUST_STFLUX
4690 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, &
4691 & Nfrec(ng),2,NT(ng))
4692# endif
4693 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
4694 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
4695 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
4696# else
4697 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
4698 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
4699# endif
4700 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
4701# ifdef ADJUST_BOUNDARY
4702# ifdef SOLVE3D
4703 real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, &
4704 & Nbrec(ng),2,NT(ng))
4705 real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
4706 real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
4707# endif
4708 real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
4709 real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
4710 real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
4711# endif
4712# ifdef ADJUST_WSTRESS
4713 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
4714 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
4715# endif
4716# ifdef SOLVE3D
4717# ifdef ADJUST_STFLUX
4718 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, &
4719 & Nfrec(ng),2,NT(ng))
4720# endif
4721 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
4722 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
4723 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
4724# else
4725 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
4726 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
4727# endif
4728 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
4729# ifdef ADJUST_BOUNDARY
4730# ifdef SOLVE3D
4731 real(r8), intent(inout) :: nl_t_obc(LBij:UBij,N(ng),4, &
4732 & Nbrec(ng),2,NT(ng))
4733 real(r8), intent(inout) :: nl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
4734 real(r8), intent(inout) :: nl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
4735# endif
4736 real(r8), intent(inout) :: nl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
4737 real(r8), intent(inout) :: nl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
4738 real(r8), intent(inout) :: nl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
4739# endif
4740# ifdef ADJUST_WSTRESS
4741 real(r8), intent(inout) :: nl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
4742 real(r8), intent(inout) :: nl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
4743# endif
4744# ifdef SOLVE3D
4745# ifdef ADJUST_STFLUX
4746 real(r8), intent(inout) :: nl_tflux(LBi:UBi,LBj:UBj, &
4747 & Nfrec(ng),2,NT(ng))
4748# endif
4749 real(r8), intent(inout) :: nl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
4750 real(r8), intent(inout) :: nl_u(LBi:UBi,LBj:UBj,N(ng),2)
4751 real(r8), intent(inout) :: nl_v(LBi:UBi,LBj:UBj,N(ng),2)
4752# else
4753 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,:)
4754 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,:)
4755# endif
4756 real(r8), intent(inout) :: nl_zeta(LBi:UBi,LBj:UBj,:)
4757# endif
4758!
4759! Local variable declarations.
4760!
4761 integer :: i, ingood, lstr, nvec, rec, status, varid
4762 integer :: L1
4763 integer :: start(4), total(4)
4764!
4765 real(r8) :: fac, fac1, fac2
4766
4767 real(r8), dimension(NpostI) :: RitzErr
4768
4769 real(r8), dimension(0:NstateVar(ng)) :: dot
4770!
4771 character (len=256) :: ncname
4772
4773 character (len=*), parameter :: MyFile = &
4774 & __FILE__//", posterior_eofs"
4775
4776# include "set_bounds.h"
4777!
4778 calledfrom=myfile
4779 sourcefile=myfile
4780!
4781!-----------------------------------------------------------------------
4782! Calculate converged eigenvectors of the Hessian.
4783!-----------------------------------------------------------------------
4784!
4785! NOTE: In the case of weak constraint ("WEAK_CONSTRAINT") and
4786! and time convolutions ("TIME_CONV"), the state arrays
4787! tl_ubar, tl_vbar, ad_ubar, and ad_vbar are only passed
4788! as required by the "state" operators routines but they
4789! are not used in subsequent calculations.
4790!
4791! Count and collect the converged eigenvalues.
4792!
4793 ingood=0
4794 DO i=innloop,1,-1
4795 ingood=ingood+1
4796 ritz(ingood)=ae_ritz(i,outloop)
4797 ritzerr(ingood)=ae_ritzerr(i,outloop)
4798 END DO
4799 nconvritz=ingood
4800!
4801! Write out number of converged eigenvalues.
4802!
4803 SELECT CASE (hss(ng)%IOtype)
4804 CASE (io_nf90)
4805 CALL netcdf_put_ivar (ng, model, hss(ng)%name, &
4806 & 'nConvRitz', nconvritz, &
4807 & (/0/), (/0/), &
4808 & ncid = hss(ng)%ncid)
4809
4810# if defined PIO_LIB && defined DISTRIBUTE
4811 CASE (io_pio)
4812 CALL pio_netcdf_put_ivar (ng, model, hss(ng)%name, &
4813 & 'nConvRitz', nconvritz, &
4814 & (/0/), (/0/), &
4815 & piofile = hss(ng)%pioFile)
4816# endif
4817 END SELECT
4818 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4819!
4820!-----------------------------------------------------------------------
4821! First, premultiply the eigenvectors of the tridiagonal
4822! matrix T(k) by the matrix of Lanczos vectors Q(k). Use tangent
4823! linear (index Lwrk) and adjoint (index Lold) state arrays as
4824! temporary storage.
4825!-----------------------------------------------------------------------
4826!
4827 IF (master) WRITE (stdout,10)
4828!
4829 hss(ng)%Rindex=0
4830!
4831 columns : DO nvec=innloop,1,-1
4832!
4833! Initialize adjoint state arrays: ad_var(Lold) = fac
4834!
4835 fac=0.0_r8
4836
4837 CALL state_initialize (ng, tile, &
4838 & lbi, ubi, lbj, ubj, lbij, ubij, &
4839 & lold, fac, &
4840# ifdef MASKING
4841 & rmask, umask, vmask, &
4842# endif
4843# ifdef ADJUST_BOUNDARY
4844# ifdef SOLVE3D
4845 & ad_t_obc, ad_u_obc, ad_v_obc, &
4846# endif
4847 & ad_ubar_obc, ad_vbar_obc, &
4848 & ad_zeta_obc, &
4849# endif
4850# ifdef ADJUST_WSTRESS
4851 & ad_ustr, ad_vstr, &
4852# endif
4853# ifdef SOLVE3D
4854# ifdef ADJUST_STFLUX
4855 & ad_tflux, &
4856# endif
4857 & ad_t, ad_u, ad_v, &
4858# else
4859 & ad_ubar, ad_vbar, &
4860# endif
4861 & ad_zeta)
4862!
4863! Compute Hessian eigenvectors.
4864!
4865 ncname=adm(ng)%name
4866!
4867 rows : DO rec=1,innloop
4868!
4869! Read gradient solution and load it into TANGENT LINEAR STATE ARRAYS
4870! at index Lwrk.
4871!
4872 CALL state_read (ng, tile, model, adm(ng)%IOtype, &
4873 & lbi, ubi, lbj, ubj, lbij, ubij, &
4874 & lwrk, rec, &
4875 & ndefadj(ng), adm(ng)%ncid, ncname, &
4876# ifdef MASKING
4877 & rmask, umask, vmask, &
4878# endif
4879# ifdef ADJUST_BOUNDARY
4880# ifdef SOLVE3D
4881 & tl_t_obc, tl_u_obc, tl_v_obc, &
4882# endif
4883 & tl_ubar_obc, tl_vbar_obc, &
4884 & tl_zeta_obc, &
4885# endif
4886# ifdef ADJUST_WSTRESS
4887 & tl_ustr, tl_vstr, &
4888# endif
4889# ifdef SOLVE3D
4890# ifdef ADJUST_STFLUX
4891 & tl_tflux, &
4892# endif
4893 & tl_t, tl_u, tl_v, &
4894# else
4895 & tl_ubar, tl_vbar, &
4896# endif
4897 & tl_zeta)
4898 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4899!
4900! Compute Hessian eigenvectors: (See NOTE above)
4901!
4902! ad_var(Lold) = fac1 * ad_var(Lold) + fac2 * tl_var(Lwrk)
4903!
4904 fac1=1.0_r8
4905 fac2=ae_zv(rec,nvec)
4906
4907 CALL state_addition (ng, tile, &
4908 & lbi, ubi, lbj, ubj, lbij, ubij, &
4909 & lold, lwrk, lold, fac1, fac2, &
4910# ifdef MASKING
4911 & rmask, umask, vmask, &
4912# endif
4913# ifdef ADJUST_BOUNDARY
4914# ifdef SOLVE3D
4915 & ad_t_obc, tl_t_obc, &
4916 & ad_u_obc, tl_u_obc, &
4917 & ad_v_obc, tl_v_obc, &
4918# endif
4919 & ad_ubar_obc, tl_ubar_obc, &
4920 & ad_vbar_obc, tl_vbar_obc, &
4921 & ad_zeta_obc, tl_zeta_obc, &
4922# endif
4923# ifdef ADJUST_WSTRESS
4924 & ad_ustr, tl_ustr, &
4925 & ad_vstr, tl_vstr, &
4926# endif
4927# ifdef SOLVE3D
4928# ifdef ADJUST_STFLUX
4929 & ad_tflux, tl_tflux, &
4930# endif
4931 & ad_t, tl_t, &
4932 & ad_u, tl_u, &
4933 & ad_v, tl_v, &
4934# if defined WEAK_CONSTRAINT && defined TIME_CONV
4935 & ad_ubar, tl_ubar, &
4936 & ad_vbar, tl_vbar, &
4937# endif
4938# else
4939 & ad_ubar, tl_ubar, &
4940 & ad_vbar, tl_vbar, &
4941# endif
4942 & ad_zeta, tl_zeta)
4943 END DO rows
4944!
4945! Write eigenvectors into Hessian NetCDF.
4946!
4947 lwrtstate2d(ng)=.true.
4948# ifdef DISTRIBUTE
4949 CALL wrt_hessian (ng, myrank, lold, lold)
4950# else
4951 CALL wrt_hessian (ng, -1, lold, lold)
4952# endif
4953 lwrtstate2d(ng)=.false.
4954 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
4955
4956 END DO columns
4957!
4958!-----------------------------------------------------------------------
4959! Second, orthonormalize the converged Hessian vectors against each
4960! other. Use tangent linear state arrays (index Lwrk) as temporary
4961! storage.
4962!-----------------------------------------------------------------------
4963!
4964! Use nl_var(1) as temporary storage since we need to preserve
4965! ad_var(Lnew).
4966!
4967 ncname=hss(ng)%name
4968 IF (master) WRITE (stdout,30)
4969!
4970 DO nvec=1,innloop
4971!
4972! Read in just computed Hessian eigenvectors into adjoint state array
4973! index Lold.
4974!
4975 CALL state_read (ng, tile, model, hss(ng)%IOtype, &
4976 & lbi, ubi, lbj, ubj, lbij, ubij, &
4977 & lold, nvec, &
4978 & 0, hss(ng)%ncid, ncname, &
4979# ifdef MASKING
4980 & rmask, umask, vmask, &
4981# endif
4982# ifdef ADJUST_BOUNDARY
4983# ifdef SOLVE3D
4984 & ad_t_obc, ad_u_obc, ad_v_obc, &
4985# endif
4986 & ad_ubar_obc, ad_vbar_obc, &
4987 & ad_zeta_obc, &
4988# endif
4989# ifdef ADJUST_WSTRESS
4990 & ad_ustr, ad_vstr, &
4991# endif
4992# ifdef SOLVE3D
4993# ifdef ADJUST_STFLUX
4994 & ad_tflux, &
4995# endif
4996 & ad_t, ad_u, ad_v, &
4997# else
4998 & ad_ubar, ad_vbar, &
4999# endif
5000 & ad_zeta)
5001 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
5002!
5003! Initialize nonlinear state arrays index L1 with just read Hessian
5004! vector in index Lold (initialize the summation): (See NOTE above)
5005!
5006! Copy ad_var(Lold) into nl_var(L1).
5007!
5008 l1=1
5009
5010 CALL state_copy (ng, tile, &
5011 & lbi, ubi, lbj, ubj, lbij, ubij, &
5012 & lold, l1, &
5013# ifdef ADJUST_BOUNDARY
5014# ifdef SOLVE3D
5015 & nl_t_obc, ad_t_obc, &
5016 & nl_u_obc, ad_u_obc, &
5017 & nl_v_obc, ad_v_obc, &
5018# endif
5019 & nl_ubar_obc, ad_ubar_obc, &
5020 & nl_vbar_obc, ad_vbar_obc, &
5021 & nl_zeta_obc, ad_zeta_obc, &
5022# endif
5023# ifdef ADJUST_WSTRESS
5024 & nl_ustr, ad_ustr, &
5025 & nl_vstr, ad_vstr, &
5026# endif
5027# ifdef SOLVE3D
5028# ifdef ADJUST_STFLUX
5029 & nl_tflux, ad_tflux, &
5030# endif
5031 & nl_t, ad_t, &
5032 & nl_u, ad_u, &
5033 & nl_v, ad_v, &
5034# if defined WEAK_CONSTRAINT && defined TIME_CONV
5035 & nl_ubar, ad_ubar, &
5036 & nl_vbar, ad_vbar, &
5037# endif
5038# else
5039 & nl_ubar, ad_ubar, &
5040 & nl_vbar, ad_vbar, &
5041# endif
5042 & nl_zeta, ad_zeta)
5043!
5044! Orthogonalize Hessian eigenvectors against each other.
5045!
5046 DO rec=1,nvec-1
5047!
5048! Read in gradient just computed Hessian eigenvectors into tangent
5049! linear state array index Lwrk.
5050!
5051 CALL state_read (ng, tile, model, hss(ng)%IOtype, &
5052 & lbi, ubi, lbj, ubj, lbij, ubij, &
5053 & lwrk, rec, &
5054 & 0, hss(ng)%ncid, ncname, &
5055# ifdef MASKING
5056 & rmask, umask, vmask, &
5057# endif
5058# ifdef ADJUST_BOUNDARY
5059# ifdef SOLVE3D
5060 & tl_t_obc, tl_u_obc, tl_v_obc, &
5061# endif
5062 & tl_ubar_obc, tl_vbar_obc, &
5063 & tl_zeta_obc, &
5064# endif
5065# ifdef ADJUST_WSTRESS
5066 & tl_ustr, tl_vstr, &
5067# endif
5068# ifdef SOLVE3D
5069# ifdef ADJUST_STFLUX
5070 & tl_tflux, &
5071# endif
5072 & tl_t, tl_u, tl_v, &
5073# else
5074 & tl_ubar, tl_vbar, &
5075# endif
5076 & tl_zeta)
5077 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
5078!
5079! Compute dot product.
5080!
5081 CALL state_dotprod (ng, tile, model, &
5082 & lbi, ubi, lbj, ubj, lbij, ubij, &
5083 & nstatevar(ng), dot(0:), &
5084# ifdef MASKING
5085 & rmask, umask, vmask, &
5086# endif
5087# ifdef ADJUST_BOUNDARY
5088# ifdef SOLVE3D
5089 & ad_t_obc(:,:,:,:,lold,:), &
5090 & tl_t_obc(:,:,:,:,lwrk,:), &
5091 & ad_u_obc(:,:,:,:,lold), &
5092 & tl_u_obc(:,:,:,:,lwrk), &
5093 & ad_v_obc(:,:,:,:,lold), &
5094 & tl_v_obc(:,:,:,:,lwrk), &
5095# endif
5096 & ad_ubar_obc(:,:,:,lold), &
5097 & tl_ubar_obc(:,:,:,lwrk), &
5098 & ad_vbar_obc(:,:,:,lold), &
5099 & tl_vbar_obc(:,:,:,lwrk), &
5100 & ad_zeta_obc(:,:,:,lold), &
5101 & tl_zeta_obc(:,:,:,lwrk), &
5102# endif
5103# ifdef ADJUST_WSTRESS
5104 & ad_ustr(:,:,:,lold), tl_ustr(:,:,:,lwrk), &
5105 & ad_vstr(:,:,:,lold), tl_vstr(:,:,:,lwrk), &
5106# endif
5107# ifdef SOLVE3D
5108# ifdef ADJUST_STFLUX
5109 & ad_tflux(:,:,:,lold,:), &
5110 & tl_tflux(:,:,:,lwrk,:), &
5111# endif
5112 & ad_t(:,:,:,lold,:), tl_t(:,:,:,lwrk,:), &
5113 & ad_u(:,:,:,lold), tl_u(:,:,:,lwrk), &
5114 & ad_v(:,:,:,lold), tl_v(:,:,:,lwrk), &
5115# else
5116 & ad_ubar(:,:,lold), tl_ubar(:,:,lwrk), &
5117 & ad_vbar(:,:,lold), tl_vbar(:,:,lwrk), &
5118# endif
5119 & ad_zeta(:,:,lold), tl_zeta(:,:,lwrk))
5120!
5121! Orthogonalize Hessian eigenvectors: (See NOTE above)
5122!
5123! nl_var(L1) = fac1 * nl_var(L1) + fac2 * tl_var(Lwrk)
5124!
5125 fac1=1.0_r8
5126 fac2=-dot(0)
5127
5128 CALL state_addition (ng, tile, &
5129 & lbi, ubi, lbj, ubj, lbij, ubij, &
5130 & l1, lwrk, l1, fac1, fac2, &
5131# ifdef MASKING
5132 & rmask, umask, vmask, &
5133# endif
5134# ifdef ADJUST_BOUNDARY
5135# ifdef SOLVE3D
5136 & nl_t_obc, tl_t_obc, &
5137 & nl_u_obc, tl_u_obc, &
5138 & nl_v_obc, tl_v_obc, &
5139# endif
5140 & nl_ubar_obc, tl_ubar_obc, &
5141 & nl_vbar_obc, tl_vbar_obc, &
5142 & nl_zeta_obc, tl_zeta_obc, &
5143# endif
5144# ifdef ADJUST_WSTRESS
5145 & nl_ustr, tl_ustr, &
5146 & nl_vstr, tl_vstr, &
5147# endif
5148# ifdef SOLVE3D
5149# ifdef ADJUST_STFLUX
5150 & nl_tflux, tl_tflux, &
5151# endif
5152 & nl_t, tl_t, &
5153 & nl_u, tl_u, &
5154 & nl_v, tl_v, &
5155# if defined WEAK_CONSTRAINT && defined TIME_CONV
5156 & nl_ubar, tl_ubar, &
5157 & nl_vbar, tl_vbar, &
5158# endif
5159# else
5160 & nl_ubar, tl_ubar, &
5161 & nl_vbar, tl_vbar, &
5162# endif
5163 & nl_zeta, tl_zeta)
5164 END DO
5165!
5166! Compute normalization factor.
5167!
5168 CALL state_dotprod (ng, tile, model, &
5169 & lbi, ubi, lbj, ubj, lbij, ubij, &
5170 & nstatevar(ng), dot(0:), &
5171# ifdef MASKING
5172 & rmask, umask, vmask, &
5173# endif
5174# ifdef ADJUST_BOUNDARY
5175# ifdef SOLVE3D
5176 & nl_t_obc(:,:,:,:,l1,:), &
5177 & nl_t_obc(:,:,:,:,l1,:), &
5178 & nl_u_obc(:,:,:,:,l1), &
5179 & nl_u_obc(:,:,:,:,l1), &
5180 & nl_v_obc(:,:,:,:,l1), &
5181 & nl_v_obc(:,:,:,:,l1), &
5182# endif
5183 & nl_ubar_obc(:,:,:,l1), &
5184 & nl_ubar_obc(:,:,:,l1), &
5185 & nl_vbar_obc(:,:,:,l1), &
5186 & nl_vbar_obc(:,:,:,l1), &
5187 & nl_zeta_obc(:,:,:,l1), &
5188 & nl_zeta_obc(:,:,:,l1), &
5189# endif
5190# ifdef ADJUST_WSTRESS
5191 & nl_ustr(:,:,:,l1), nl_ustr(:,:,:,l1), &
5192 & nl_vstr(:,:,:,l1), nl_vstr(:,:,:,l1), &
5193# endif
5194# ifdef SOLVE3D
5195# ifdef ADJUST_STFLUX
5196 & nl_tflux(:,:,:,l1,:), &
5197 & nl_tflux(:,:,:,l1,:), &
5198# endif
5199 & nl_t(:,:,:,l1,:), nl_t(:,:,:,l1,:), &
5200 & nl_u(:,:,:,l1), nl_u(:,:,:,l1), &
5201 & nl_v(:,:,:,l1), nl_v(:,:,:,l1), &
5202# else
5203 & nl_ubar(:,:,l1), nl_ubar(:,:,l1), &
5204 & nl_vbar(:,:,l1), nl_vbar(:,:,l1), &
5205# endif
5206 & nl_zeta(:,:,l1), nl_zeta(:,:,l1))
5207!
5208! Normalize Hessian eigenvectors:
5209!
5210! nl_var(L1) = fac * nl_var(L1)
5211!
5212 fac=1.0_r8/sqrt(dot(0))
5213
5214 CALL state_scale (ng, tile, &
5215 & lbi, ubi, lbj, ubj, lbij, ubij, &
5216 & l1, l1, fac, &
5217# ifdef MASKING
5218 & rmask, umask, vmask, &
5219# endif
5220# ifdef ADJUST_BOUNDARY
5221# ifdef SOLVE3D
5222 & nl_t_obc, nl_u_obc, nl_v_obc, &
5223# endif
5224 & nl_ubar_obc, nl_vbar_obc, &
5225 & nl_zeta_obc, &
5226# endif
5227# ifdef ADJUST_WSTRESS
5228 & nl_ustr, nl_vstr, &
5229# endif
5230# ifdef SOLVE3D
5231# ifdef ADJUST_STFLUX
5232 & nl_tflux, &
5233# endif
5234 & nl_t, nl_u, nl_v, &
5235# else
5236 & nl_ubar, nl_vbar, &
5237# endif
5238 & nl_zeta)
5239!
5240! Copy nl_var(L1) into ad_var(Lold). See NOTE above.
5241!
5242 CALL state_copy (ng, tile, &
5243 & lbi, ubi, lbj, ubj, lbij, ubij, &
5244 & l1, lold, &
5245# ifdef ADJUST_BOUNDARY
5246# ifdef SOLVE3D
5247 & ad_t_obc, nl_t_obc, &
5248 & ad_u_obc, nl_u_obc, &
5249 & ad_v_obc, nl_v_obc, &
5250# endif
5251 & ad_ubar_obc, nl_ubar_obc, &
5252 & ad_vbar_obc, nl_vbar_obc, &
5253 & ad_zeta_obc, nl_zeta_obc, &
5254# endif
5255# ifdef ADJUST_WSTRESS
5256 & ad_ustr, nl_ustr, &
5257 & ad_vstr, nl_vstr, &
5258# endif
5259# ifdef SOLVE3D
5260# ifdef ADJUST_STFLUX
5261 & ad_tflux, nl_tflux, &
5262# endif
5263 & ad_t, nl_t, &
5264 & ad_u, nl_u, &
5265 & ad_v, nl_v, &
5266# if defined WEAK_CONSTRAINT && defined TIME_CONV
5267 & ad_ubar, nl_ubar, &
5268 & ad_vbar, nl_vbar, &
5269# endif
5270# else
5271 & ad_ubar, nl_ubar, &
5272 & ad_vbar, nl_vbar, &
5273# endif
5274 & ad_zeta, nl_zeta)
5275!
5276! Write out converged Ritz eigenvalues and is associated accuracy.
5277!
5278 SELECT CASE (hss(ng)%IOtype)
5279 CASE (io_nf90)
5280 CALL netcdf_put_fvar (ng, model, hss(ng)%name, &
5281 & 'Ritz', ritz(nvec:), &
5282 & (/nvec/), (/1/), &
5283 & ncid = hss(ng)%ncid)
5284 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
5285
5286 CALL netcdf_put_fvar (ng, model, hss(ng)%name, &
5287 & 'Ritz_error', ritzerr(nvec:), &
5288 & (/nvec/), (/1/), &
5289 & ncid = hss(ng)%ncid)
5290 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
5291
5292# if defined PIO_LIB && defined DISTRIBUTE
5293 CASE (io_pio)
5294 CALL pio_netcdf_put_fvar (ng, model, hss(ng)%name, &
5295 & 'Ritz', ritz(nvec:), &
5296 & (/nvec/), (/1/), &
5297 & piofile = hss(ng)%pioFile)
5298 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
5299
5300 CALL pio_netcdf_put_fvar (ng, model, hss(ng)%name, &
5301 & 'Ritz_error', ritzerr(nvec:), &
5302 & (/nvec/), (/1/), &
5303 & piofile = hss(ng)%pioFile)
5304 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
5305# endif
5306 END SELECT
5307!
5308! Replace record "nvec" of Hessian eigenvectors NetCDF with the
5309! normalized value in adjoint state arrays at index Lold.
5310!
5311 hss(ng)%Rindex=nvec-1
5312 lwrtstate2d(ng)=.true.
5313# ifdef DISTRIBUTE
5314 CALL wrt_hessian (ng, myrank, lold, lold)
5315# else
5316 CALL wrt_hessian (ng, -1, lold, lold)
5317# endif
5318 lwrtstate2d(ng)=.false.
5319 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
5320
5321 END DO
5322!
5323 10 FORMAT (/,' Computing converged analysis error eofs...',/)
5324 20 FORMAT (a,'_',i3.3,'.nc')
5325 30 FORMAT (/,' Orthonormalizing converged analysis error eofs...',/)
5326!
5327 RETURN
5328 END SUBROUTINE posterior_eofs
5329#endif
5330
5331 END MODULE posterior_mod
distribute_mod::mp_bcastf
Definition distribute.F:75
distribute_mod::mp_bcasti
Definition distribute.F:95
mod_netcdf::netcdf_put_fvar
Definition mod_netcdf.F:82
mod_netcdf::netcdf_put_ivar
Definition mod_netcdf.F:96
mod_pio_netcdf::pio_netcdf_put_fvar
Definition mod_pio_netcdf.F:94
mod_pio_netcdf::pio_netcdf_put_ivar
Definition mod_pio_netcdf.F:108
distribute_mod
Definition distribute.F:3
lapack_mod
Definition lapack_mod.F:2
lapack_mod::dsteqr
subroutine, public dsteqr(compz, n, d, e, z, ldz, work, info)
Definition lapack_mod.F:66
mod_boundary
Definition mod_boundary.F:2
mod_boundary::boundary
type(t_boundary), dimension(:), allocatable boundary
Definition mod_boundary.F:746
mod_coupling
Definition mod_coupling.F:2
mod_forces
Definition mod_forces.F:2
mod_forces::forces
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
mod_fourdvar
Definition mod_fourdvar.F:2
mod_fourdvar::cg_beta
real(dp), dimension(:,:), allocatable cg_beta
Definition mod_fourdvar.F:654
mod_fourdvar::ritz
real(dp), dimension(:), allocatable ritz
Definition mod_fourdvar.F:664
mod_fourdvar::ae_gnorm
real(r8), dimension(:), allocatable ae_gnorm
Definition mod_fourdvar.F:745
mod_fourdvar::hevecerr
real(r8) hevecerr
Definition mod_fourdvar.F:858
mod_fourdvar::ae_beta
real(r8), dimension(:,:), allocatable ae_beta
Definition mod_fourdvar.F:739
mod_fourdvar::ae_delta
real(r8), dimension(:,:), allocatable ae_delta
Definition mod_fourdvar.F:738
mod_fourdvar::nposti
integer nposti
Definition mod_fourdvar.F:452
mod_fourdvar::ae_trace
real(r8), dimension(:), allocatable ae_trace
Definition mod_fourdvar.F:741
mod_fourdvar::ritzmaxerr
real(dp) ritzmaxerr
Definition mod_fourdvar.F:659
mod_fourdvar::ae_tmatrix
real(r8), dimension(:,:), allocatable ae_tmatrix
Definition mod_fourdvar.F:749
mod_fourdvar::cg_delta
real(dp), dimension(:,:), allocatable cg_delta
Definition mod_fourdvar.F:691
mod_fourdvar::ae_ritz
real(r8), dimension(:,:), allocatable ae_ritz
Definition mod_fourdvar.F:755
mod_fourdvar::ae_ritzerr
real(r8), dimension(:,:), allocatable ae_ritzerr
Definition mod_fourdvar.F:756
mod_fourdvar::ae_zv
real(r8), dimension(:,:), allocatable ae_zv
Definition mod_fourdvar.F:740
mod_fourdvar::nstatevar
integer, dimension(:), allocatable nstatevar
Definition mod_fourdvar.F:362
mod_fourdvar::nconvritz
integer nconvritz
Definition mod_fourdvar.F:842
mod_grid
Definition mod_grid.F:2
mod_grid::grid
type(t_grid), dimension(:), allocatable grid
Definition mod_grid.F:365
mod_iounits
Definition mod_iounits.F:2
mod_iounits::calledfrom
character(len=256) calledfrom
Definition mod_iounits.F:315
mod_iounits::adm
type(t_io), dimension(:), allocatable adm
Definition mod_iounits.F:200
mod_iounits::hss
type(t_io), dimension(:), allocatable hss
Definition mod_iounits.F:223
mod_iounits::stdout
integer stdout
Definition mod_iounits.F:195
mod_iounits::sourcefile
character(len=256) sourcefile
Definition mod_iounits.F:317
mod_ncparam
Definition mod_ncparam.F:2
mod_ncparam::io_nf90
integer, parameter io_nf90
Definition mod_ncparam.F:95
mod_ncparam::isvvel
integer isvvel
Definition mod_ncparam.F:507
mod_ncparam::io_pio
integer, parameter io_pio
Definition mod_ncparam.F:96
mod_ncparam::isvbar
integer isvbar
Definition mod_ncparam.F:505
mod_ncparam::istvar
integer, dimension(:), allocatable istvar
Definition mod_ncparam.F:520
mod_ncparam::isuvel
integer isuvel
Definition mod_ncparam.F:506
mod_ncparam::isfsur
integer isfsur
Definition mod_ncparam.F:503
mod_ncparam::isubar
integer isubar
Definition mod_ncparam.F:504
mod_netcdf
Definition mod_netcdf.F:2
mod_ocean
Definition mod_ocean.F:2
mod_ocean::ocean
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
mod_parallel
Definition mod_parallel.F:2
mod_parallel::master
logical master
Definition mod_parallel.F:40
mod_parallel::myrank
integer myrank
Definition mod_parallel.F:157
mod_param
Definition mod_param.F:2
mod_param::domain
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
mod_pio_netcdf
Definition mod_pio_netcdf.F:2
mod_scalars
Definition mod_scalars.F:2
mod_scalars::ninner
integer ninner
Definition mod_scalars.F:353
mod_scalars::lobc
logical, dimension(:,:,:), allocatable lobc
Definition mod_scalars.F:1127
mod_scalars::iwest
integer, parameter iwest
Definition mod_scalars.F:1432
mod_scalars::lstflux
logical, dimension(:,:), allocatable lstflux
Definition mod_scalars.F:1156
mod_scalars::exit_flag
integer exit_flag
Definition mod_scalars.F:559
mod_scalars::lwrtstate2d
logical, dimension(:), allocatable lwrtstate2d
Definition mod_scalars.F:1020
mod_scalars::isouth
integer, parameter isouth
Definition mod_scalars.F:1433
mod_scalars::ieast
integer, parameter ieast
Definition mod_scalars.F:1434
mod_scalars::inorth
integer, parameter inorth
Definition mod_scalars.F:1435
mod_scalars::ndefadj
integer, dimension(:), allocatable ndefadj
Definition mod_scalars.F:1057
mod_scalars::inner
integer inner
Definition mod_scalars.F:360
mod_scalars::noerror
integer noerror
Definition mod_scalars.F:561
mod_stepping
Definition mod_stepping.F:2
mod_stepping::lold
integer, dimension(:), allocatable lold
Definition mod_stepping.F:107
mod_stepping::lnew
integer, dimension(:), allocatable lnew
Definition mod_stepping.F:106
posterior_mod
Definition posterior.F:2
posterior_mod::posterior_eofs
subroutine posterior_eofs(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, lold, lnew, lwrk, innloop, outloop, rmask, umask, vmask, nl_t_obc, nl_u_obc, nl_v_obc, nl_ubar_obc, nl_vbar_obc, nl_zeta_obc, nl_ustr, nl_vstr, nl_tflux, nl_t, nl_u, nl_v, nl_ubar, nl_vbar, nl_zeta, tl_t_obc, tl_u_obc, tl_v_obc, tl_ubar_obc, tl_vbar_obc, tl_zeta_obc, tl_ustr, tl_vstr, tl_tflux, tl_t, tl_u, tl_v, tl_ubar, tl_vbar, tl_zeta, ad_t_obc, ad_u_obc, ad_v_obc, ad_ubar_obc, ad_vbar_obc, ad_zeta_obc, ad_ustr, ad_vstr, ad_tflux, ad_t, ad_u, ad_v, ad_ubar, ad_vbar, ad_zeta)
Definition posterior.F:4571
posterior_mod::posterior
subroutine, public posterior(ng, tile, model, innloop, outloop, ltrace)
Definition posterior.F:65
posterior_mod::new_direction
subroutine new_direction(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, lold, lnew, rmask, umask, vmask, ad_t_obc, ad_u_obc, ad_v_obc, ad_ubar_obc, ad_vbar_obc, ad_zeta_obc, ad_ustr, ad_vstr, ad_tflux, ad_t, ad_u, ad_v, ad_zeta, d_t_obc, d_u_obc, d_v_obc, d_ubar_obc, d_vbar_obc, d_zeta_obc, d_sustr, d_svstr, d_stflx, d_t, d_u, d_v, d_zeta)
Definition posterior.F:2033
posterior_mod::lanczos
subroutine lanczos(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, lold, lnew, lwrk, innloop, outloop, rmask, umask, vmask, tl_t_obc, tl_u_obc, tl_v_obc, tl_ubar_obc, tl_vbar_obc, tl_zeta_obc, tl_ustr, tl_vstr, tl_tflux, tl_t, tl_u, tl_v, tl_ubar, tl_vbar, tl_zeta, ad_t_obc, ad_u_obc, ad_v_obc, ad_ubar_obc, ad_vbar_obc, ad_zeta_obc, ad_ustr, ad_vstr, ad_tflux, ad_t, ad_u, ad_v, ad_ubar, ad_vbar, ad_zeta)
Definition posterior.F:3849
posterior_mod::tl_new_vector
subroutine tl_new_vector(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, linp, lout, innloop, outloop, rmask, umask, vmask, d_t_obc, d_u_obc, d_v_obc, d_ubar_obc, d_vbar_obc, d_zeta_obc, d_sustr, d_svstr, d_stflx, d_t, d_u, d_v, d_zeta, tl_t_obc, tl_u_obc, tl_v_obc, tl_ubar_obc, tl_vbar_obc, tl_zeta_obc, tl_ustr, tl_vstr, tl_tflux, tl_t, tl_u, tl_v, tl_zeta, ad_t_obc, ad_u_obc, ad_v_obc, ad_ubar_obc, ad_vbar_obc, ad_zeta_obc, ad_ustr, ad_vstr, ad_tflux, ad_t, ad_u, ad_v, ad_zeta)
Definition posterior.F:1304
posterior_mod::posterior_tile
subroutine posterior_tile(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, lold, lnew, innloop, outloop, ltrace, rmask, umask, vmask, nl_t_obc, nl_u_obc, nl_v_obc, nl_ubar_obc, nl_vbar_obc, nl_zeta_obc, nl_ustr, nl_vstr, nl_tflux, nl_t, nl_u, nl_v, nl_ubar, nl_vbar, nl_zeta, tl_t_obc, tl_u_obc, tl_v_obc, tl_ubar_obc, tl_vbar_obc, tl_zeta_obc, tl_ustr, tl_vstr, tl_tflux, tl_t, tl_u, tl_v, tl_ubar, tl_vbar, tl_zeta, d_t_obc, d_u_obc, d_v_obc, d_ubar_obc, d_vbar_obc, d_zeta_obc, d_sustr, d_svstr, d_stflx, d_t, d_u, d_v, d_ubar, d_vbar, d_zeta, ad_t_obc, ad_u_obc, ad_v_obc, ad_ubar_obc, ad_vbar_obc, ad_zeta_obc, ad_ustr, ad_vstr, ad_tflux, ad_t, ad_u, ad_v, ad_ubar, ad_vbar, ad_zeta)
Definition posterior.F:317
posterior_mod::analysis_error
subroutine analysis_error(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, lold, lnew, lwrk, innloop, outloop, rmask, umask, vmask, ad_t_obc, ad_u_obc, ad_v_obc, ad_ubar_obc, ad_vbar_obc, ad_zeta_obc, ad_ustr, ad_vstr, ad_tflux, ad_t, ad_u, ad_v, ad_ubar, ad_vbar, ad_zeta, tl_t_obc, tl_u_obc, tl_v_obc, tl_ubar_obc, tl_vbar_obc, tl_zeta_obc, tl_ustr, tl_vstr, tl_tflux, tl_t, tl_u, tl_v, tl_ubar, tl_vbar, tl_zeta, nl_t_obc, nl_u_obc, nl_v_obc, nl_ubar_obc, nl_vbar_obc, nl_zeta_obc, nl_ustr, nl_vstr, nl_tflux, nl_t, nl_u, nl_v, nl_ubar, nl_vbar, nl_zeta)
Definition posterior.F:2733
state_addition_mod
Definition state_addition.F:2
state_addition_mod::state_addition
subroutine, public state_addition(ng, tile, lbi, ubi, lbj, ubj, lbij, ubij, lin1, lin2, lout, fac1, fac2, rmask, umask, vmask, s1_t_obc, s2_t_obc, s1_u_obc, s2_u_obc, s1_v_obc, s2_v_obc, s1_ubar_obc, s2_ubar_obc, s1_vbar_obc, s2_vbar_obc, s1_zeta_obc, s2_zeta_obc, s1_sustr, s2_sustr, s1_svstr, s2_svstr, s1_tflux, s2_tflux, s1_t, s2_t, s1_u, s2_u, s1_v, s2_v, s1_ubar, s2_ubar, s1_vbar, s2_vbar, s1_zeta, s2_zeta)
Definition state_addition.F:64
state_copy_mod
Definition state_copy.F:2
state_copy_mod::state_copy
subroutine, public state_copy(ng, tile, lbi, ubi, lbj, ubj, lbij, ubij, linp, lout, s1_t_obc, s2_t_obc, s1_u_obc, s2_u_obc, s1_v_obc, s2_v_obc, s1_ubar_obc, s2_ubar_obc, s1_vbar_obc, s2_vbar_obc, s1_zeta_obc, s2_zeta_obc, s1_sustr, s2_sustr, s1_svstr, s2_svstr, s1_tflux, s2_tflux, s1_t, s2_t, s1_u, s2_u, s1_v, s2_v, s1_ubar, s2_ubar, s1_vbar, s2_vbar, s1_zeta, s2_zeta)
Definition state_copy.F:57
state_dotprod_mod
Definition state_dotprod.F:2
state_dotprod_mod::state_dotprod
subroutine, public state_dotprod(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, nstatevars, dotprod, rmask, umask, vmask, s1_t_obc, s2_t_obc, s1_u_obc, s2_u_obc, s1_v_obc, s2_v_obc, s1_ubar_obc, s2_ubar_obc, s1_vbar_obc, s2_vbar_obc, s1_zeta_obc, s2_zeta_obc, s1_sustr, s2_sustr, s1_svstr, s2_svstr, s1_tflux, s2_tflux, s1_t, s2_t, s1_u, s2_u, s1_v, s2_v, s1_zeta, s2_zeta)
Definition state_dotprod.F:85
state_initialize_mod
Definition state_initialize.F:2
state_initialize_mod::state_initialize
subroutine, public state_initialize(ng, tile, lbi, ubi, lbj, ubj, lbij, ubij, lout, fac, rmask, umask, vmask, s_t_obc, s_u_obc, s_v_obc, s_ubar_obc, s_vbar_obc, s_zeta_obc, s_sustr, s_svstr, s_tflux, s_t, s_u, s_v, s_zeta)
Definition state_initialize.F:52
state_read_mod
Definition state_read.F:2
state_read_mod::state_read
subroutine, public state_read(ng, tile, model, iotype, lbi, ubi, lbj, ubj, lbij, ubij, lout, rec, nopen, ncid, piofile, ncname, rmask, umask, vmask, s_t_obc, s_u_obc, s_v_obc, s_ubar_obc, s_vbar_obc, s_zeta_obc, s_ustr, s_vstr, s_tflux, s_t, s_u, s_v, s_zeta)
Definition state_read.F:137
state_scale_mod
Definition state_scale.F:2
state_scale_mod::state_scale
subroutine, public state_scale(ng, tile, lbi, ubi, lbj, ubj, lbij, ubij, linp, lout, fac, rmask, umask, vmask, s_t_obc, s_u_obc, s_v_obc, s_ubar_obc, s_vbar_obc, s_zeta_obc, s_sustr, s_svstr, s_tflux, s_t, s_u, s_v, s_zeta)
Definition state_scale.F:50
strings_mod
Definition strings.F:2
strings_mod::founderror
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
wrt_hessian_mod
Definition wrt_hessian.F:2
wrt_hessian_mod::wrt_hessian
subroutine, public wrt_hessian(ng, tile, kout, nout)
Definition wrt_hessian.F:62
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3