Ticket #25: cgradient.h

File cgradient.h, 113.7 KB (added by m.hadfield, 19 years ago)
Line 
1!
2!svn $Id: cgradient.h 48 2007-05-09 16:15:44Z arango $
3!================================================== Hernan G. Arango ===
4! Copyright (c) 2002-2007 The ROMS/TOMS Group Andrew M. Moore !
5! Licensed under a MIT/X style license !
6! See License_ROMS.txt !
7!=======================================================================
8! !
9! This module minimizes incremental 4Dvar quadratic cost function !
10! using a preconditioned version of the conjugate gradient algorithm !
11! proposed by Mike Fisher (ECMWF). !
12! !
13! In the following, M represents the preconditioner. Specifically, !
14! !
15! M = I + SUM_i [ (mu_i-1) h_i (h_i)'], !
16! !
17! where mu_i can take the following values: !
18! !
19! Lscale= 1: mu_i = lambda_i !
20! Lscale=-1: mu_i = 1 / lambda_i !
21! Lscale= 2: mu_i = SQRT (lambda_i) !
22! Lscale=-2: mu_i = 1 / SQRT(lambda_i) !
23! !
24! where lambda_i are the Hessian eigenvalues and h_i are the Hessian !
25! eigenvectors. For Lscale=-1 the preconditioner is an approximation !
26! of the inverse Hessian matrix constructed from the leading Hessian !
27! eigenvectors. A full description is given in Fisher and Courtier !
28! (1995), ECMWF Tech Memo 220. !
29! !
30! Given an initial model state X(0), gradient G(0), descent direction !
31! d(0), and trial step size tau(1), the minimization algorithm at the !
32! k-iteration is : !
33! !
34! (1) Run tangent linear model initialized with trial step, Xhat(k): !
35! !
36! Xhat(k) = X(k) + tau(k) * d(k) (Eq 5a) !
37! !
38! (2) Run adjoint model to compute gradient at trial point, Ghat(k): !
39! !
40! Ghat(k) = GRAD[ f(Xhat(k)) ] (Eq 5b) !
41! !
42! (3) Compute optimum step size, alpha(k): !
43! !
44! alpha(k) = tau(k) * <d(k), G(k)> / !
45! !
46! (<d(k),G(k)> - <d(k), Ghat(k)>) (Eq 5c) !
47! !
48! here <...> denotes dot product !
49! !
50! (4) Compute new starting point (TLM increments), X(k+1): !
51! !
52! X(k+1) = X(k) + alpha(k) * d(k) (Eq 5d) !
53! !
54! (5) Compute gradient at new point, G(k+1): !
55! !
56! G(k+1) = G(k) + (alpha(k) / tau(k)) * (Ghat(k) - G(k)) (Eq 5e) !
57! !
58! overwrite G(k+1) in the NetCDF for latter use. !
59! !
60! (6) Orthogonalize new gradient, G(k+1), against all previous !
61! gradients [G(k), ..., G(0)], in reverse order, using the !
62! modified Gramm-Schmidt algorithm. Notice that we need to !
63! save all inner loop gradient solutions. !
64! !
65! For the preconditioned case the appropriate inner-product !
66! for the orthonormalizatio is <G,MG>. !
67! !
68! (7) Compute new descent direction, d(k+1): !
69! !
70! beta(k+1) = <G(k+1),M G(k+1)> / <G(k),M G(k)> (Eq 5g) !
71! !
72! d(k+1) = - MG(k+1) + beta(k+1) * d(k) (Eq 5f) !
73! !
74! After the first iteration, the trial step size is: !
75! !
76! tau(k) = alpha(k-1) !
77! !
78! NOTE: In all of the following computations we are using the NLM !
79! state variable arrays as temporary arrays. !
80! !
81! Reference: !
82! !
83! Fisher, M., 1997: Efficient Minimization of Quadratic Penalty !
84! funtions, unpublish manuscript, 1-14. !
85! !
86!=======================================================================
87!
88 implicit none
89
90 PRIVATE
91 PUBLIC :: cgradient
92
93 CONTAINS
94!
95!***********************************************************************
96 SUBROUTINE cgradient (ng, tile, model, innLoop, outLoop)
97!***********************************************************************
98!
99 USE mod_param
100#ifdef SOLVE3D
101 USE mod_coupling
102#endif
103#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
104 USE mod_forces
105#endif
106 USE mod_grid
107 USE mod_ocean
108 USE mod_stepping
109!
110! Imported variable declarations.
111!
112 integer, intent(in) :: ng, tile, model, innLoop, outLoop
113!
114! Local variable declarations.
115!
116#include "tile.h"
117!
118#ifdef PROFILE
119 CALL wclock_on (ng, model, 36)
120#endif
121 CALL cgradient_tile (ng, model, Istr, Iend, Jstr, Jend, &
122 & LBi, UBi, LBj, UBj, &
123 & Lold(ng), Lnew(ng), &
124 & innLoop, outLoop, &
125#ifdef MASKING
126 & GRID(ng) % rmask, &
127 & GRID(ng) % umask, &
128 & GRID(ng) % vmask, &
129#endif
130#ifdef ADJUST_WSTRESS
131 & FORCES(ng) % sustrG, &
132 & FORCES(ng) % svstrG, &
133#endif
134#ifdef SOLVE3D
135# ifdef ADJUST_STFLUX
136 & FORCES(ng) % tflux, &
137# endif
138 & OCEAN(ng) % t, &
139 & OCEAN(ng) % u, &
140 & OCEAN(ng) % v, &
141#else
142 & OCEAN(ng) % ubar, &
143 & OCEAN(ng) % vbar, &
144#endif
145 & OCEAN(ng) % zeta, &
146#ifdef ADJUST_WSTRESS
147 & OCEAN(ng) % d_sustr, &
148 & OCEAN(ng) % d_svstr, &
149#endif
150#ifdef SOLVE3D
151# ifdef ADJUST_STFLUX
152 & OCEAN(ng) % d_stflx, &
153# endif
154 & OCEAN(ng) % d_t, &
155 & OCEAN(ng) % d_u, &
156 & OCEAN(ng) % d_v, &
157#else
158 & OCEAN(ng) % d_ubar, &
159 & OCEAN(ng) % d_vbar, &
160#endif
161 & OCEAN(ng) % d_zeta, &
162#ifdef ADJUST_WSTRESS
163 & FORCES(ng) % ad_ustr, &
164 & FORCES(ng) % ad_vstr, &
165#endif
166#ifdef SOLVE3D
167# ifdef ADJUST_STFLUX
168 & FORCES(ng) % ad_tflux, &
169# endif
170 & OCEAN(ng) % ad_t, &
171 & OCEAN(ng) % ad_u, &
172 & OCEAN(ng) % ad_v, &
173#else
174 & OCEAN(ng) % ad_ubar, &
175 & OCEAN(ng) % ad_vbar, &
176#endif
177 & OCEAN(ng) % ad_zeta, &
178#ifdef ADJUST_WSTRESS
179 & FORCES(ng) % tl_ustr, &
180 & FORCES(ng) % tl_vstr, &
181#endif
182#ifdef SOLVE3D
183# ifdef ADJUST_STFLUX
184 & FORCES(ng) % tl_tflux, &
185# endif
186 & OCEAN(ng) % tl_t, &
187 & OCEAN(ng) % tl_u, &
188 & OCEAN(ng) % tl_v, &
189#else
190 & OCEAN(ng) % tl_ubar, &
191 & OCEAN(ng) % tl_vbar, &
192#endif
193 & OCEAN(ng) % tl_zeta)
194#ifdef PROFILE
195 CALL wclock_on (ng, model, 36)
196#endif
197 RETURN
198 END SUBROUTINE cgradient
199!
200!***********************************************************************
201 SUBROUTINE cgradient_tile (ng, model, Istr, Iend, Jstr, Jend, &
202 & LBi, UBi, LBj, UBj, &
203 & Lold, Lnew, &
204 & innLoop, outLoop, &
205#ifdef MASKING
206 & rmask, umask, vmask, &
207#endif
208#ifdef ADJUST_WSTRESS
209 & nl_ustr, nl_vstr, &
210#endif
211#ifdef SOLVE3D
212# ifdef ADJUST_STFLUX
213 & nl_tflux, &
214# endif
215 & nl_t, nl_u, nl_v, &
216#else
217 & nl_ubar, nl_vbar, &
218#endif
219 & nl_zeta, &
220#ifdef ADJUST_WSTRESS
221 & d_sustr, d_svstr, &
222#endif
223#ifdef SOLVE3D
224# ifdef ADJUST_STFLUX
225 & d_stflx, &
226# endif
227 & d_t, d_u, d_v, &
228#else
229 & d_ubar, d_vbar, &
230#endif
231 & d_zeta, &
232#ifdef ADJUST_WSTRESS
233 & ad_ustr, ad_vstr, &
234#endif
235#ifdef SOLVE3D
236# ifdef ADJUST_STFLUX
237 & ad_tflux, &
238# endif
239 & ad_t, ad_u, ad_v, &
240#else
241 & ad_ubar, ad_vbar, &
242#endif
243 & ad_zeta, &
244#ifdef ADJUST_WSTRESS
245 & tl_ustr, tl_vstr, &
246#endif
247#ifdef SOLVE3D
248# ifdef ADJUST_STFLUX
249 & tl_tflux, &
250# endif
251 & tl_t, tl_u, tl_v, &
252#else
253 & tl_ubar, tl_vbar, &
254#endif
255 & tl_zeta)
256!***********************************************************************
257!
258 USE mod_param
259 USE mod_parallel
260 USE mod_fourdvar
261 USE mod_iounits
262 USE mod_ncparam
263 USE mod_netcdf
264 USE mod_scalars
265
266#ifdef DISTRIBUTE
267!
268 USE distribute_mod, ONLY : mp_bcastf, mp_bcasti
269#endif
270!
271! Imported variable declarations.
272!
273 integer, intent(in) :: ng, model, Iend, Istr, Jend, Jstr
274 integer, intent(in) :: LBi, UBi, LBj, UBj
275 integer, intent(in) :: Lold, Lnew
276 integer, intent(in) :: innLoop, outLoop
277!
278#ifdef ASSUMED_SHAPE
279# ifdef MASKING
280 real(r8), intent(in) :: rmask(LBi:,LBj:)
281 real(r8), intent(in) :: umask(LBi:,LBj:)
282 real(r8), intent(in) :: vmask(LBi:,LBj:)
283# endif
284# ifdef ADJUST_WSTRESS
285 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:)
286 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:)
287# endif
288# ifdef SOLVE3D
289# ifdef ADJUST_STFLUX
290 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:)
291# endif
292 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
293 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
294 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
295# else
296 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
297 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
298# endif
299 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
300# ifdef ADJUST_WSTRESS
301 real(r8), intent(inout) :: d_sustr(LBi:,LBj:)
302 real(r8), intent(inout) :: d_svstr(LBi:,LBj:)
303# endif
304# ifdef SOLVE3D
305# ifdef ADJUST_STFLUX
306 real(r8), intent(inout) :: d_stflx(LBi:,LBj:,:)
307# endif
308 real(r8), intent(inout) :: d_t(LBi:,LBj:,:,:)
309 real(r8), intent(inout) :: d_u(LBi:,LBj:,:)
310 real(r8), intent(inout) :: d_v(LBi:,LBj:,:)
311# else
312 real(r8), intent(inout) :: d_ubar(LBi:,LBj:)
313 real(r8), intent(inout) :: d_vbar(LBi:,LBj:)
314# endif
315 real(r8), intent(inout) :: d_zeta(LBi:,LBj:)
316# ifdef ADJUST_WSTRESS
317 real(r8), intent(inout) :: nl_ustr(LBi:,LBj:,:)
318 real(r8), intent(inout) :: nl_vstr(LBi:,LBj:,:)
319# endif
320# ifdef SOLVE3D
321# ifdef ADJUST_STFLUX
322 real(r8), intent(inout) :: nl_tflux(LBi:,LBj:,:,:)
323# endif
324 real(r8), intent(inout) :: nl_t(LBi:,LBj:,:,:,:)
325 real(r8), intent(inout) :: nl_u(LBi:,LBj:,:,:)
326 real(r8), intent(inout) :: nl_v(LBi:,LBj:,:,:)
327# else
328 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
329 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
330# endif
331 real(r8), intent(inout) :: nl_zeta(LBi:,LBj:,:)
332# ifdef ADJUST_WSTRESS
333 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:)
334 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:)
335# endif
336# ifdef SOLVE3D
337# ifdef ADJUST_STFLUX
338 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:)
339# endif
340 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
341 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
342 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
343# else
344 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
345 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
346# endif
347 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
348
349#else
350
351# ifdef MASKING
352 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
353 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
354 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
355# endif
356# ifdef ADJUST_WSTRESS
357 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,2)
358 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,2)
359# endif
360# ifdef SOLVE3D
361# ifdef ADJUST_STFLUX
362 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
363# endif
364 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
365 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
366 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
367# else
368 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,3)
369 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,3)
370# endif
371 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,3)
372# ifdef ADJUST_WSTRESS
373 real(r8), intent(inout) :: d_sustr(LBi:UBi,LBj:UBj)
374 real(r8), intent(inout) :: d_svstr(LBi:UBi,LBj:UBj)
375# endif
376# ifdef SOLVE3D
377# ifdef ADJUST_STFLUX
378 real(r8), intent(inout) :: d_stflx(LBi:UBi,LBj:UBj,NT(ng))
379# endif
380 real(r8), intent(inout) :: d_t(LBi:UBi,LBj:UBj,N(ng),NT(ng))
381 real(r8), intent(inout) :: d_u(LBi:UBi,LBj:UBj,N(ng))
382 real(r8), intent(inout) :: d_v(LBi:UBi,LBj:UBj,N(ng))
383# else
384 real(r8), intent(inout) :: d_ubar(LBi:UBi,LBj:UBj)
385 real(r8), intent(inout) :: d_vbar(LBi:UBi,LBj:UBj)
386# endif
387 real(r8), intent(inout) :: d_zeta(LBi:UBi,LBj:UBj)
388# ifdef ADJUST_WSTRESS
389 real(r8), intent(inout) :: nl_ustr(LBi:UBi,LBj:UBj,2)
390 real(r8), intent(inout) :: nl_vstr(LBi:UBi,LBj:UBj,2)
391# endif
392# ifdef SOLVE3D
393# ifdef ADJUST_STFLUX
394 real(r8), intent(inout) :: nl_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
395# endif
396 real(r8), intent(inout) :: nl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
397 real(r8), intent(inout) :: nl_u(LBi:UBi,LBj:UBj,N(ng),2)
398 real(r8), intent(inout) :: nl_v(LBi:UBi,LBj:UBj,N(ng),2)
399# else
400 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,3)
401 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,3)
402# endif
403 real(r8), intent(inout) :: nl_zeta(LBi:UBi,LBj:UBj,3)
404# ifdef ADJUST_WSTRESS
405 real(r8), intent(inout) :: nl_tl_ustr(LBi:UBi,LBj:UBj,2)
406 real(r8), intent(inout) :: nl_tl_vstr(LBi:UBi,LBj:UBj,2)
407# endif
408# ifdef SOLVE3D
409# ifdef ADJUST_STFLUX
410 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
411# endif
412 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
413 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
414 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
415# else
416 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,3)
417 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,3)
418# endif
419 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,3)
420#endif
421!
422! Local variable declarations.
423!
424 integer :: Linp, Lout, Lwrk, L1, L2, i, Lscale
425 integer :: status, varid
426
427 real(r8) :: norm
428
429 real(r8), dimension(0:NstateVar(ng)) :: Adjust
430 real(r8), dimension(0:NstateVar(ng)) :: dot_old, dot_new
431 real(r8), dimension(0:NstateVar(ng)) :: old_dot, new_dot
432!
433!-----------------------------------------------------------------------
434! Initialize trial step size.
435!-----------------------------------------------------------------------
436!
437 IF (innLoop.eq.0) THEN
438 cg_tau(innLoop,outLoop)=CGstepI
439 cg_alpha(innLoop,outLoop)=cg_tau(innLoop,outLoop)
440 DO i=0,NstateVar(ng)
441 dot_old(i)=0.0_r8
442 dot_new(i)=0.0_r8
443 old_dot(i)=0.0_r8
444 new_dot(i)=0.0_r8
445 FOURDVAR(ng)%CostGradDot(i)=0.0_r8
446 END DO
447 END IF
448 IF (Master) THEN
449 WRITE (stdout,10)
450 10 FORMAT (/,' <<<< Descent Algorithm >>>>')
451 END IF
452!
453! If preconditioning, read in number of converged eigenvectors and their
454! associated eigenvalues.
455!
456 IF (Lprecond.and.((innLoop.eq.0).and.(outLoop.eq.1))) THEN
457 IF (InpThread) THEN
458 status=nf_inq_varid(ncHSSid(ng),'nConvRitz',varid)
459 status=nf_get_var1_int(ncHSSid(ng), varid, 1, nConvRitz)
460 IF (status.ne.nf_noerr) THEN
461 WRITE (stdout,20) 'nConvRitz', TRIM(HSSname(ng))
462 20 FORMAT (/,' CGRADIENT - error while reading variable: ', &
463 & a,/,12x,'from NetCDF file: ',a)
464 exit_flag=2
465 ioerror=status
466 RETURN
467 END IF
468 status=nf_inq_varid(ncHSSid(ng),'Ritz',varid)
469 status=nf_get_vara_TYPE(ncHSSid(ng), varid, 1, nConvRitz, Ritz)
470 IF (status.ne.nf_noerr) THEN
471 WRITE (stdout,20) 'Ritz', TRIM(HSSname(ng))
472 exit_flag=2
473 ioerror=status
474 RETURN
475 END IF
476 END IF
477#ifdef DISTRIBUTE
478 CALL mp_bcasti (ng, iADM, nConvRitz, 1)
479 CALL mp_bcastf (ng, iADM, Ritz, nConvRitz)
480#endif
481 END IF
482!
483!-----------------------------------------------------------------------
484! Compute conjugate gradient optimum step size, alpha(k).
485!-----------------------------------------------------------------------
486!
487 IF (innLoop.gt.0) THEN
488 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
489 & LBi, UBi, LBj, UBj, &
490 & NstateVar(ng), dot_old(0:), &
491#ifdef MASKING
492 & rmask, umask, vmask, &
493#endif
494#ifdef ADJUST_WSTRESS
495 & d_sustr, ad_ustr(:,:,Lold), &
496 & d_svstr, ad_vstr(:,:,Lold), &
497#endif
498#ifdef SOLVE3D
499# ifdef ADJUST_STFLUX
500 & d_stflx, ad_tflux(:,:,:,Lold,:), &
501# endif
502 & d_t, ad_t(:,:,:,Lold,:), &
503 & d_u, ad_u(:,:,:,Lold), &
504 & d_v, ad_v(:,:,:,Lold), &
505#else
506 & d_ubar, ad_ubar(:,:,Lold), &
507 & d_vbar, ad_vbar(:,:,Lold), &
508#endif
509 & d_zeta, ad_zeta(:,:,Lold))
510!
511! If preconditioning, compute new dot product, <d(k), H^-1 * Ghat(k)>.
512!
513 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
514 & LBi, UBi, LBj, UBj, &
515 & NstateVar(ng), dot_new(0:), &
516#ifdef MASKING
517 & rmask, umask, vmask, &
518#endif
519#ifdef ADJUST_WSTRESS
520 & d_sustr, ad_ustr(:,:,Lnew), &
521 & d_svstr, ad_vstr(:,:,Lnew), &
522#endif
523#ifdef SOLVE3D
524# ifdef ADJUST_STFLUX
525 & d_stflx, ad_tflux(:,:,Lnew,:), &
526# endif
527 & d_t, ad_t(:,:,:,Lnew,:), &
528 & d_u, ad_u(:,:,:,Lnew), &
529 & d_v, ad_v(:,:,:,Lnew), &
530#else
531 & d_ubar, ad_ubar(:,:,Lnew), &
532 & d_vbar, ad_vbar(:,:,Lnew), &
533#endif
534 & d_zeta, ad_zeta(:,:,Lnew))
535!
536! Compute new optimal step size.
537!
538 cg_tau(innLoop,outLoop)=cg_alpha(innLoop-1,outLoop)
539 cg_alpha(innLoop,outLoop)=cg_tau(innLoop,outLoop)* &
540 & dot_old(0)/(dot_old(0)-dot_new(0))
541 END IF
542!
543! Adjust the cost function for the previous inner-loop iteration.
544! This is based on a first-order Taylor expansion of the cost function.
545! Let vhat=v+tau*d. During each inner-loop the tangent linear
546! model provides J(vhat). What we require is J(v). Using a 1st-order
547! Taylor expansion we have: J(vhat)=J(v)+tau*<d,grad> where grad is
548! the cost function gradient computed during the last inner-loop
549! immediately prior to the orthogonalization. Rearranging this
550! equation we have: J(v)=J(vhat)-tau*<d,grad>. In the code
551! J(vhat)=CostFun(:) and <d,grad>=CostFunDot(:). Remember though
552! that J(v) is the cost function associated with v from the previous
553! inner-loop.
554!
555 DO i=0,NstateVar(ng)
556 Adjust(i)=cg_tau(innLoop,outLoop)*FOURDVAR(ng)%CostGradDot(i)
557 FOURDVAR(ng)%CostFun(i)=FOURDVAR(ng)%CostFun(i)-Adjust(i)
558 END DO
559!
560!-----------------------------------------------------------------------
561! Estimate the gradient for the new state vector, G(k+1).
562!-----------------------------------------------------------------------
563!
564! If preconditioning, compute old dot product, <G(k), H^-1 * G(k)>.
565! The ADM arrays, index Lold, will be used a as temporary storage
566! after this.
567!
568 IF (Lprecond) THEN
569 Lscale=-1
570 Lwrk=2
571 CALL precond (ng, model, Istr, Iend, Jstr, Jend, &
572 & LBi, UBi, LBj, UBj, &
573 & NstateVar(ng), Lold, Lwrk, Lscale, &
574 & nConvRitz, Ritz, &
575#ifdef MASKING
576 & rmask, umask, vmask, &
577#endif
578#ifdef ADJUST_WSTRESS
579 & ad_ustr, ustr, tl_ustr, &
580 & ad_vstr, vstr, tl_vstr, &
581#endif
582#ifdef SOLVE3D
583# ifdef ADJUST_STFLUX
584 & ad_tflux, nl_tflux, tl_tflx, &
585# endif
586 & ad_t, nl_t, tl_t, &
587 & ad_u, nl_u, tl_u, &
588 & ad_v, nl_v, tl_v, &
589#else
590 & ad_ubar, nl_ubar, tl_ubar, &
591 & ad_vbar, nl_vbar, tl_vbar, &
592#endif
593 & ad_zeta, nl_zeta, tl_zeta)
594!
595 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
596 & LBi, UBi, LBj, UBj, &
597 & NstateVar(ng), old_dot(0:), &
598#ifdef MASKING
599 & rmask, umask, vmask, &
600#endif
601#ifdef ADJUST_WSTRESS
602 & ad_ustr(:,:,Lold), tl_ustr(:,:,Lwrk), &
603 & ad_vstr(:,:,Lold), tl_vstr(:,:,Lwrk), &
604#endif
605#ifdef SOLVE3D
606# ifdef ADJUST_STFLUX
607 & ad_tflux(:,:,Lold,:), tl_tflux(:,:,Lwrk,:), &
608# endif
609 & ad_t(:,:,:,Lold,:), tl_t(:,:,:,Lwrk,:), &
610 & ad_u(:,:,:,Lold), tl_u(:,:,:,Lwrk), &
611 & ad_v(:,:,:,Lold), tl_v(:,:,:,Lwrk), &
612#else
613 & ad_ubar(:,:,Lold), tl_ubar(:,:,Lwrk), &
614 & ad_vbar(:,:,Lold), tl_vbar(:,:,Lwrk), &
615#endif
616 & ad_zeta(:,:,Lold), tl_zeta(:,:,Lwrk))
617!
618! If not preconditioning, compute old dot product, <G(k), G(k)>.
619! The ADM arrays, index Lold, will be used a as temporary storage
620! after this.
621!
622 ELSE
623 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
624 & LBi, UBi, LBj, UBj, &
625 & NstateVar(ng), old_dot(0:), &
626#ifdef MASKING
627 & rmask, umask, vmask, &
628#endif
629#ifdef ADJUST_WSTRESS
630 & ad_ustr(:,:,Lold), ad_ustr(:,:,Lold), &
631 & ad_vstr(:,:,Lold), ad_vstr(:,:,Lold), &
632#endif
633#ifdef SOLVE3D
634# ifdef ADJUST_STFLUX
635 & ad_tflux(:,:,Lold,:), ad_tflux(:,:,Lold,:), &
636# endif
637 & ad_t(:,:,:,Lold,:), ad_t(:,:,:,Lold,:), &
638 & ad_u(:,:,:,Lold), ad_u(:,:,:,Lold), &
639 & ad_v(:,:,:,Lold), ad_v(:,:,:,Lold), &
640#else
641 & ad_ubar(:,:,Lold), ad_ubar(:,:,Lold), &
642 & ad_vbar(:,:,Lold), ad_vbar(:,:,Lold), &
643#endif
644 & ad_zeta(:,:,Lold), ad_zeta(:,:,Lold))
645 END IF
646!
647! Notice that the current gradient Ghat(k) in time index Lnew is
648! overwritten with the new gradient G(k+1).
649!
650! G(k+1) = G(k) + (alpha(k) / tau(k)) * (Ghat(k) - G(k))
651! Lnew Lold Lnew Lold index
652!
653! Also save G(k+1) in time index Lold as a non-orthogonalized new
654! gradient.
655!
656 CALL ad_new_state (ng, Istr, Iend, Jstr, Jend, &
657 & LBi, UBi, LBj, UBj, &
658 & Lold, Lnew, &
659 & cg_alpha(innLoop,outLoop), &
660 & cg_tau(innLoop,outLoop), &
661#ifdef MASKING
662 & rmask, umask, vmask, &
663#endif
664#ifdef ADJUST_WSTRESS
665 & ad_ustr, ad_vstr, &
666#endif
667#ifdef SOLVE3D
668# ifdef ADJUST_STFLUX
669 & ad_tflux, &
670# endif
671 & ad_t, ad_u, ad_v, &
672#else
673 & ad_ubar, ad_vbar, &
674#endif
675 & ad_zeta)
676
677#ifdef ORTHOGONALIZATION
678!
679! Orthogonalize new gradient, G(k+1), against all previous gradients
680! G(0) to G(k). Use TLM state arrays at time index Lwrk=2, to load
681! each of the previous gradients.
682!
683 IF (innLoop.gt.0) THEN
684 Lwrk=2
685 CALL orthogonalize (ng, model, Istr, Iend, Jstr, Jend, &
686 & LBi, UBi, LBj, UBj, &
687 & Lold, Lnew, Lwrk, &
688 & innLoop, outLoop, &
689# ifdef MASKING
690 & rmask, umask, vmask, &
691# endif
692# ifdef SOLVE3D
693# ifdef ADJUST_STFLUX
694 & nl_tflux, &
695# endif
696 & nl_t, nl_u, nl_v, &
697# else
698 & nl_ubar, nl_vbar, &
699# endif
700 & nl_zeta, &
701# ifdef ADJUST_WSTRESS
702 & tl_ustr, tl_vstr, &
703# endif
704# ifdef SOLVE3D
705# ifdef ADJUST_STFLUX
706 & tl_tflux, &
707# endif
708 & tl_t, tl_u, tl_v, &
709# else
710 & tl_ubar, tl_vbar, &
711# endif
712 & tl_zeta, &
713# ifdef ADJUST_WSTRESS
714 & ustr, vstr, &
715# endif
716# ifdef ADJUST_WSTRESS
717 & ad_ustr, ad_vstr, &
718# endif
719# ifdef SOLVE3D
720# ifdef ADJUST_STFLUX
721 & ad_tflux, &
722# endif
723 & ad_t, ad_u, ad_v, &
724# else
725 & ad_ubar, ad_vbar, &
726# endif
727 & ad_zeta)
728 END IF
729#endif
730!
731!-----------------------------------------------------------------------
732! Compute new starting tangent linear state vector, X(k+1).
733!-----------------------------------------------------------------------
734!
735! Here we are doing step (4), equation 5d, the new TLM increment for
736! the initial conditions are always saved at time level Lout=1.
737!
738! X(k+1) = X(k) + alpha(k) * d(k)
739! Lout Linp index
740!
741 IF (innLoop.gt.0) THEN
742 Linp=1
743 Lout=1
744 CALL tl_new_state (ng, Istr, Iend, Jstr, Jend, &
745 & LBi, UBi, LBj, UBj, &
746 & Linp, Lout, &
747 & cg_alpha(innLoop,outLoop), &
748#ifdef MASKING
749 & rmask, umask, vmask, &
750#endif
751#ifdef ADJUST_WSTRESS
752 & d_sustr, d_svstr, &
753#endif
754#ifdef SOLVE3D
755# ifdef ADJUST_STFLUX
756 & d_stflx, &
757# endif
758 & d_t, d_u, d_v, &
759#else
760 & d_ubar, d_vbar, &
761#endif
762 & d_zeta, &
763#ifdef ADJUST_WSTRESS
764 & tl_ustr, tl_vstr, &
765#endif
766#ifdef SOLVE3D
767# ifdef ADJUST_STFLUX
768 & tl_tflux, &
769# endif
770 & tl_t, tl_u, tl_v, &
771#else
772 & tl_ubar, tl_vbar, &
773#endif
774 & tl_zeta)
775!
776! If last iteration of inner loop, skip remaining computations. The
777! TLM increments computed here are the ones that are needed update
778! the NLM model initial conditions.
779!
780!! IF (innLoop.eq.Ninner) RETURN
781 END IF
782!
783!-----------------------------------------------------------------------
784! Compute new conjugate descent direction, d(k+1).
785!-----------------------------------------------------------------------
786!
787! If preconditioning, multiply the new gradient by H^-1 and save in
788! nl_var(Lwrk).
789!
790 IF (Lprecond) THEN
791 Lscale=-1
792 Lwrk=2
793 CALL precond (ng, model, Istr, Iend, Jstr, Jend, &
794 & LBi, UBi, LBj, UBj, &
795 & NstateVar(ng), Lnew, Lwrk, Lscale, &
796 & nConvRitz, Ritz, &
797#ifdef MASKING
798 & rmask, umask, vmask, &
799#endif
800#ifdef ADJUST_WSTRESS
801 & ad_ustr, nl_ustr, nl_ustr, &
802 & ad_vstr, nl_vstr, nl_vstr, &
803#endif
804#ifdef SOLVE3D
805# ifdef ADJUST_STFLUX
806 & ad_tflux, tflux, tflux &
807# endif
808 & ad_t, nl_t, nl_t, &
809 & ad_u, nl_u, nl_u, &
810 & ad_v, nl_v, nl_v, &
811#else
812 & ad_ubar, nl_ubar, nl_ubar, &
813 & ad_vbar, nl_vbar, nl_vbar, &
814#endif
815 & ad_zeta, nl_zeta, nl_zeta)
816 END IF
817!
818! If preconditioning, compute new dot product, <G(k+1), H^-1 * G(k+1)>.
819!
820 IF (innLoop.gt.0) THEN
821 IF (Lprecond) THEN
822 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
823 & LBi, UBi, LBj, UBj, &
824 & NstateVar(ng), new_dot(0:), &
825#ifdef MASKING
826 & rmask, umask, vmask, &
827#endif
828#ifdef ADJUST_WSTRESS
829 & ad_ustr(:,:,Lnew), nl_ustr(:,:,Lwrk), &
830 & ad_vstr(:,:,Lnew), nl_vstr(:,:,Lwrk), &
831#endif
832#ifdef SOLVE3D
833# ifdef ADJUST_STFLUX
834 & ad_tflux(:,:,Lnew,:),nl_tflux(:,:,Lwrk,:),&
835# endif
836 & ad_t(:,:,:,Lnew,:), nl_t(:,:,:,Lwrk,:), &
837 & ad_u(:,:,:,Lnew), nl_u(:,:,:,Lwrk), &
838 & ad_v(:,:,:,Lnew), nl_v(:,:,:,Lwrk), &
839#else
840 & ad_ubar(:,:,Lnew), nl_ubar(:,:,Lwrk), &
841 & ad_vbar(:,:,Lnew), nl_vbar(:,:,Lwrk), &
842#endif
843 & ad_zeta(:,:,Lnew), nl_zeta(:,:,Lwrk))
844 ELSE
845!
846! If not preconditioning, compute new dot product, <G(k+1), G(k+1)>.
847!
848 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
849 & LBi, UBi, LBj, UBj, &
850 & NstateVar(ng), new_dot(0:), &
851#ifdef MASKING
852 & rmask, umask, vmask, &
853#endif
854#ifdef ADJUST_WSTRESS
855 & ad_ustr(:,:,Lnew), ad_ustr(:,:,Lnew), &
856 & ad_vstr(:,:,Lnew), ad_vstr(:,:,Lnew), &
857#endif
858#ifdef SOLVE3D
859# ifdef ADJUST_STFLUX
860 & ad_tflux(:,:,Lnew,:), &
861 & ad_tflux(:,:,Lnew,:), &
862# endif
863 & ad_t(:,:,:,Lnew,:), ad_t(:,:,:,Lnew,:), &
864 & ad_u(:,:,:,Lnew), ad_u(:,:,:,Lnew), &
865 & ad_v(:,:,:,Lnew), ad_v(:,:,:,Lnew), &
866#else
867 & ad_ubar(:,:,Lnew), ad_ubar(:,:,Lnew), &
868 & ad_vbar(:,:,Lnew), ad_vbar(:,:,Lnew), &
869#endif
870 & ad_zeta(:,:,Lnew), ad_zeta(:,:,Lnew))
871 END IF
872!
873! Compute conjugate direction coefficient, beta(k+1).
874!
875 cg_beta(innLoop,outLoop)=new_dot(0)/old_dot(0)
876 ELSE
877 cg_beta(innLoop,outLoop)=0.0_r8
878 END IF
879!
880! If preconditioning, compute new conjugate direction, d(k+1). Notice
881! that the preconditined gradient is in NLM (index Lwrk) state arrays.
882!
883 IF (Lprecond) THEN
884 CALL new_direction (ng, model, Istr, Iend, Jstr, Jend, &
885 & LBi, UBi, LBj, UBj, &
886 & Lold, Lwrk, &
887 & cg_beta(innLoop,outLoop), &
888#ifdef MASKING
889 & rmask, umask, vmask, &
890#endif
891#ifdef ADJUST_WSTRESS
892 & nl_ustr, nl_vstr, &
893#endif
894#ifdef SOLVE3D
895# ifdef ADJUST_STFLUX
896 & nl_tflux, &
897# endif
898 & nl_t, nl_u, nl_v, &
899#else
900 & nl_ubar, nl_vbar, &
901#endif
902 & nl_zeta, &
903#ifdef ADJUST_WSTRESS
904 & d_sustr, d_svstr, &
905#endif
906#ifdef SOLVE3D
907# ifdef ADJUST_STFLUX
908 & d_stflx, &
909# endif
910 & d_t, d_u, d_v, &
911#else
912 & d_ubar, d_vbar, &
913#endif
914 & d_zeta)
915!
916! If not preconditioning, compute new conjugate direction, d(k+1).
917!
918 ELSE
919 CALL new_direction (ng, model, Istr, Iend, Jstr, Jend, &
920 & LBi, UBi, LBj, UBj, &
921 & Lold, Lnew, &
922 & cg_beta(innLoop,outLoop), &
923#ifdef MASKING
924 & rmask, umask, vmask, &
925#endif
926#ifdef ADJUST_WSTRESS
927 & ad_ustr, ad_vstr, &
928#endif
929#ifdef SOLVE3D
930# ifdef ADJUST_STFLUX
931 & ad_tflux, &
932# endif
933 & ad_t, ad_u, ad_v, &
934#else
935 & ad_ubar, ad_vbar, &
936#endif
937 & ad_zeta, &
938#ifdef ADJUST_WSTRESS
939 & d_sustr, d_svstr, &
940#endif
941#ifdef SOLVE3D
942# ifdef ADJUST_STFLUX
943 & d_stflx, &
944# endif
945 & d_t, d_u, d_v, &
946#else
947 & d_ubar, d_vbar, &
948#endif
949 & d_zeta)
950 END IF
951!
952! Compute next iteration dot product, <d(k), G(k)>, using new d(k+1)
953! and non-orthogonalized G(k+1) used to adjust cost function.
954!
955 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
956 & LBi, UBi, LBj, UBj, &
957 & NstateVar(ng), &
958 & FOURDVAR(ng)%CostGradDot(0:), &
959#ifdef MASKING
960 & rmask, umask, vmask, &
961#endif
962#ifdef ADJUST_WSTRESS
963 & d_sustr, ad_ustr(:,:,Lold), &
964 & d_svstr, ad_vstr(:,:,Lold), &
965#endif
966#ifdef SOLVE3D
967# ifdef ADJUST_STFLUX
968 & d_stflx, ad_tflux(:,:,:,Lold,:), &
969# endif
970 & d_t, ad_t(:,:,:,Lold,:), &
971 & d_u, ad_u(:,:,:,Lold), &
972 & d_v, ad_v(:,:,:,Lold), &
973#else
974 & d_ubar, ad_ubar(:,:,Lold), &
975 & d_vbar, ad_vbar(:,:,Lold), &
976#endif
977 & d_zeta, ad_zeta(:,:,Lold))
978!
979!-----------------------------------------------------------------------
980! Set TLM initial conditions for next inner loop, Xhat(k+1).
981!-----------------------------------------------------------------------
982!
983! Here we are doing step (1), equation 5a, the new TLM initial
984! conditions for the next inner loop are always saved at Lout=2.
985!
986! Xhat(k+1) = X(k+1) + tau(k+1) * d(k+1), where tau(k+1)=alpha(k)
987! Lout Linp index
988!
989 Linp=1
990 Lout=2
991 CALL tl_new_state (ng, Istr, Iend, Jstr, Jend, &
992 & LBi, UBi, LBj, UBj, &
993 & Linp, Lout, &
994 & cg_alpha(innLoop,outLoop), &
995#ifdef MASKING
996 & rmask, umask, vmask, &
997#endif
998#ifdef ADJUST_WSTRESS
999 & d_sustr, d_svstr, &
1000#endif
1001#ifdef SOLVE3D
1002# ifdef ADJUST_STFLUX
1003 & d_stflx, &
1004# endif
1005 & d_t, d_u, d_v, &
1006#else
1007 & d_ubar, d_vbar, &
1008#endif
1009 & d_zeta, &
1010#ifdef ADJUST_WSTRESS
1011 & tl_ustr, tl_vstr, &
1012#endif
1013#ifdef SOLVE3D
1014# ifdef ADJUST_STFLUX
1015 & tl_tflux, &
1016# endif
1017 & tl_t, tl_u, tl_v, &
1018#else
1019 & tl_ubar, tl_vbar, &
1020#endif
1021 & tl_zeta)
1022!
1023!-----------------------------------------------------------------------
1024! Write out conjugate gradient information into NetCDF file.
1025!-----------------------------------------------------------------------
1026!
1027 CALL cg_write (ng, innLoop, outLoop)
1028!
1029! Report algorithm parameters.
1030!
1031 IF (Master) THEN
1032 WRITE (stdout,30) outLoop, innLoop, &
1033 & cg_tau(innLoop,outLoop), &
1034 & cg_alpha(innLoop,outLoop), &
1035 & cg_beta(innLoop,outLoop), &
1036 & outLoop, MAX(0,innLoop-1), Adjust(0), &
1037 & outLoop, innLoop, &
1038 & 'dot product', innLoop, innLoop, &
1039 & dot_old(0), 'alpha', &
1040 & 'dot product', innLoop, innLoop, &
1041 & dot_new(0), 'alpha', &
1042 & 'dot product', innLoop, innLoop, &
1043 & old_dot(0), 'beta', &
1044 & 'dot product', innLoop+1, innLoop+1, &
1045 & new_dot(0), 'beta'
1046 30 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1047 & 'tau = ',1p,e14.7, &
1048 & ', alpha = ',1p,e14.7, &
1049 & ', Beta = ',1p,e14.7, &
1050 & /,1x,'(',i3.3,',',i3.3,'): ', &
1051 & 'Total COST Function Adjustment = ',1p,e19.12, &
1052 & /,1x,'(',i3.3,',',i3.3,'): ', &
1053 & a,' <d(',i3.3,'),G(',i3.3,')> = ',1p,e19.12,3x,a,/,12x, &
1054 & a,' <d(',i3.3,'),g(',i3.3,')> = ',1p,e19.12,3x,a,/,12x, &
1055 & a,' <G(',i3.3,'),G(',i3.3,')> = ',1p,e19.12,3x,a,/,12x, &
1056 & a,' <G(',i3.3,'),G(',i3.3,')> = ',1p,e19.12,3x,a,/)
1057 END IF
1058
1059 RETURN
1060 END SUBROUTINE cgradient_tile
1061
1062!
1063!***********************************************************************
1064 SUBROUTINE tl_new_state (ng, Istr, Iend, Jstr, Jend, &
1065 & LBi, UBi, LBj, UBj, &
1066 & Linp, Lout, alphaK, &
1067#ifdef MASKING
1068 & rmask, umask, vmask, &
1069#endif
1070#ifdef ADJUST_WSTRESS
1071 & d_sustr, d_svstr, &
1072#endif
1073#ifdef SOLVE3D
1074# ifdef ADJUST_STFLUX
1075 & d_stflx, &
1076# endif
1077 & d_t, d_u, d_v, &
1078#else
1079 & d_ubar, d_vbar, &
1080#endif
1081 & d_zeta, &
1082#ifdef ADJUST_WSTRESS
1083 & tl_ustr, tl_vstr, &
1084#endif
1085#ifdef SOLVE3D
1086# ifdef ADJUST_STFLUX
1087 & tl_tflux, &
1088# endif
1089 & tl_t, tl_u, tl_v, &
1090#else
1091 & tl_ubar, tl_vbar, &
1092#endif
1093 & tl_zeta)
1094!***********************************************************************
1095!
1096 USE mod_param
1097!
1098! Imported variable declarations.
1099!
1100 integer, intent(in) :: ng, Iend, Istr, Jend, Jstr
1101 integer, intent(in) :: LBi, UBi, LBj, UBj
1102 integer, intent(in) :: Linp, Lout
1103
1104 real(r8), intent(in) :: alphaK
1105!
1106#ifdef ASSUMED_SHAPE
1107# ifdef MASKING
1108 real(r8), intent(in) :: rmask(LBi:,LBj:)
1109 real(r8), intent(in) :: umask(LBi:,LBj:)
1110 real(r8), intent(in) :: vmask(LBi:,LBj:)
1111# endif
1112# ifdef ADJUST_WSTRESS
1113 real(r8), intent(inout) :: d_sustr(LBi:,LBj:)
1114 real(r8), intent(inout) :: d_svstr(LBi:,LBj:)
1115# endif
1116# ifdef SOLVE3D
1117# ifdef ADJUST_STFLUX
1118 real(r8), intent(inout) :: d_stflx(LBi:,LBj:,:)
1119# endif
1120 real(r8), intent(inout) :: d_t(LBi:,LBj:,:,:)
1121 real(r8), intent(inout) :: d_u(LBi:,LBj:,:)
1122 real(r8), intent(inout) :: d_v(LBi:,LBj:,:)
1123# else
1124 real(r8), intent(inout) :: d_ubar(LBi:,LBj:)
1125 real(r8), intent(inout) :: d_vbar(LBi:,LBj:)
1126# endif
1127 real(r8), intent(inout) :: d_zeta(LBi:,LBj:)
1128# ifdef ADJUST_WSTRESS
1129 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:)
1130 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:)
1131# endif
1132# ifdef SOLVE3D
1133# ifdef ADJUST_STFLUX
1134 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:)
1135# endif
1136 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
1137 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
1138 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
1139# else
1140 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
1141 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
1142# endif
1143 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
1144#else
1145# ifdef MASKING
1146 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1147 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1148 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1149# endif
1150# ifdef ADJUST_WSTRESS
1151 real(r8), intent(inout) :: d_sustr(LBi:UBi,LBj:UBj)
1152 real(r8), intent(inout) :: d_svstr(LBi:UBi,LBj:UBj)
1153# endif
1154# ifdef SOLVE3D
1155# ifdef ADJUST_STFLUX
1156 real(r8), intent(inout) :: d_stflx(LBi:UBi,LBj:UBj,NT(ng))
1157# endif
1158 real(r8), intent(inout) :: d_t(LBi:UBi,LBj:UBj,N(ng),NT(ng))
1159 real(r8), intent(inout) :: d_u(LBi:UBi,LBj:UBj,N(ng))
1160 real(r8), intent(inout) :: d_v(LBi:UBi,LBj:UBj,N(ng))
1161# else
1162 real(r8), intent(inout) :: d_ubar(LBi:UBi,LBj:UBj)
1163 real(r8), intent(inout) :: d_vbar(LBi:UBi,LBj:UBj)
1164# endif
1165 real(r8), intent(inout) :: d_zeta(LBi:UBi,LBj:UBj)
1166# ifdef ADJUST_WSTRESS
1167 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,2)
1168 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,2)
1169# endif
1170# ifdef SOLVE3D
1171# ifdef ADJUST_STFLUX
1172 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
1173# endif
1174 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1175 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
1176 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
1177# else
1178 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,3)
1179 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,3)
1180# endif
1181 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,3)
1182#endif
1183!
1184! Local variable declarations.
1185!
1186 integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
1187 integer :: i, j
1188#ifdef SOLVE3D
1189 integer :: itrc, k
1190#endif
1191
1192#include "set_bounds.h"
1193!
1194!-----------------------------------------------------------------------
1195! Compute new starting tangent linear state vector, X(k+1).
1196!-----------------------------------------------------------------------
1197!
1198! Free-surface.
1199!
1200 DO j=JstrR,JendR
1201 DO i=IstrR,IendR
1202 tl_zeta(i,j,Lout)=tl_zeta(i,j,Linp)+ &
1203 & alphaK*d_zeta(i,j)
1204#ifdef MASKING
1205 tl_zeta(i,j,Lout)=tl_zeta(i,j,Lout)*rmask(i,j)
1206#endif
1207 END DO
1208 END DO
1209
1210#ifndef SOLVE3D
1211!
1212! 2D momentum.
1213!
1214 DO j=JstrR,JendR
1215 DO i=Istr,IendR
1216 tl_ubar(i,j,Lout)=tl_ubar(i,j,Linp)+ &
1217 & alphaK*d_ubar(i,j)
1218# ifdef MASKING
1219 tl_ubar(i,j,Lout)=tl_ubar(i,j,Lout)*umask(i,j)
1220# endif
1221 END DO
1222 END DO
1223 DO j=Jstr,JendR
1224 DO i=IstrR,IendR
1225 tl_vbar(i,j,Lout)=tl_vbar(i,j,Linp)+ &
1226 & alphaK*d_vbar(i,j)
1227# ifdef MASKING
1228 tl_vbar(i,j,Lout)=tl_vbar(i,j,Lout)*vmask(i,j)
1229# endif
1230 END DO
1231 END DO
1232#endif
1233#ifdef ADJUST_WSTRESS
1234!
1235! Surface momentum stress.
1236!
1237 DO j=JstrR,JendR
1238 DO i=Istr,IendR
1239 tl_ustr(i,j,Lout)=tl_ustr(i,j,Linp)+ &
1240 & alphaK*d_sustr(i,j)
1241# ifdef MASKING
1242 tl_ustr(i,j,Lout)=tl_ustr(i,j,Lout)*umask(i,j)
1243# endif
1244 END DO
1245 END DO
1246 DO j=Jstr,JendR
1247 DO i=IstrR,IendR
1248 tl_vstr(i,j,Lout)=tl_vstr(i,j,Linp)+ &
1249 & alphaK*d_vbar(i,j)
1250# ifdef MASKING
1251 tl_vstr(i,j,Lout)=tl_vstr(i,j,Lout)*vmask(i,j)
1252# endif
1253 END DO
1254 END DO
1255#endif
1256
1257#ifdef SOLVE3D
1258!
1259! 3D momentum.
1260!
1261 DO k=1,N(ng)
1262 DO j=JstrR,JendR
1263 DO i=Istr,IendR
1264 tl_u(i,j,k,Lout)=tl_u(i,j,k,Linp)+ &
1265 & alphaK*d_u(i,j,k)
1266# ifdef MASKING
1267 tl_u(i,j,k,Lout)=tl_u(i,j,k,Lout)*umask(i,j)
1268# endif
1269 END DO
1270 END DO
1271 DO j=Jstr,JendR
1272 DO i=IstrR,IendR
1273 tl_v(i,j,k,Lout)=tl_v(i,j,k,Linp)+ &
1274 & alphaK*d_v(i,j,k)
1275# ifdef MASKING
1276 tl_v(i,j,k,Lout)=tl_v(i,j,k,Lout)*vmask(i,j)
1277# endif
1278 END DO
1279 END DO
1280 END DO
1281!
1282! Tracers.
1283!
1284 DO itrc=1,NT(ng)
1285 DO k=1,N(ng)
1286 DO j=JstrR,JendR
1287 DO i=IstrR,IendR
1288 tl_t(i,j,k,Lout,itrc)=tl_t(i,j,k,Linp,itrc)+ &
1289 & alphaK*d_t(i,j,k,itrc)
1290# ifdef MASKING
1291 tl_t(i,j,k,Lout,itrc)=tl_t(i,j,k,Lout,itrc)*rmask(i,j)
1292# endif
1293 END DO
1294 END DO
1295 END DO
1296 END DO
1297
1298# ifdef ADJUST_STFLUX
1299!
1300! Surface tracers flux.
1301!
1302 DO itrc=1,NT(ng)
1303 DO j=JstrR,JendR
1304 DO i=IstrR,IendR
1305 tl_tflux(i,j,Lout,itrc)=tl_tflux(i,j,Linp,itrc)+ &
1306 & alphaK*d_t(i,j,k,itrc)
1307# ifdef MASKING
1308 tl_tflux(i,j,Lout,itrc)=tl_tflux(i,j,Lout,itrc)*rmask(i,j)
1309# endif
1310 END DO
1311 END DO
1312 END DO
1313# endif
1314#endif
1315
1316 RETURN
1317 END SUBROUTINE tl_new_state
1318!
1319!***********************************************************************
1320 SUBROUTINE ad_new_state (ng, Istr, Iend, Jstr, Jend, &
1321 & LBi, UBi, LBj, UBj, &
1322 & Lold, Lnew, alphaK, tauK, &
1323#ifdef MASKING
1324 & rmask, umask, vmask, &
1325#endif
1326#ifdef ADJUST_WSTRESS
1327 & ad_ustr, ad_vstr, &
1328#endif
1329#ifdef SOLVE3D
1330# ifdef ADJUST_STFLUX
1331 & ad_tflux, &
1332# endif
1333 & ad_t, ad_u, ad_v, &
1334#else
1335 & ad_ubar, ad_vbar, &
1336#endif
1337 & ad_zeta)
1338!***********************************************************************
1339!
1340 USE mod_param
1341!
1342! Imported variable declarations.
1343!
1344 integer, intent(in) :: ng, Iend, Istr, Jend, Jstr
1345 integer, intent(in) :: LBi, UBi, LBj, UBj
1346 integer, intent(in) :: Lold, Lnew
1347
1348 real(r8), intent(in) :: alphaK, tauK
1349!
1350#ifdef ASSUMED_SHAPE
1351# ifdef MASKING
1352 real(r8), intent(in) :: rmask(LBi:,LBj:)
1353 real(r8), intent(in) :: umask(LBi:,LBj:)
1354 real(r8), intent(in) :: vmask(LBi:,LBj:)
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#else
1373# ifdef MASKING
1374 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1375 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1376 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1377# endif
1378# ifdef ADJUST_WSTRESS
1379 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,2)
1380 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,2)
1381# endif
1382# ifdef SOLVE3D
1383# ifdef ADJUST_STFLUX
1384 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
1385# endif
1386 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1387 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
1388 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
1389# else
1390 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,3)
1391 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,3)
1392# endif
1393 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,3)
1394#endif
1395!
1396! Local variable declarations.
1397!
1398 integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
1399 integer :: i, j
1400#ifdef SOLVE3D
1401 integer :: itrc, k
1402#endif
1403 real(r8) :: fac
1404
1405#include "set_bounds.h"
1406!
1407!-----------------------------------------------------------------------
1408! Estimate the gradient for the new state vector, G(k+1). Notice that
1409! the Lnew record is overwritten.
1410!-----------------------------------------------------------------------
1411!
1412 fac=alphaK/tauK
1413!
1414! Free-surface.
1415!
1416 DO j=JstrR,JendR
1417 DO i=IstrR,IendR
1418 ad_zeta(i,j,Lnew)=ad_zeta(i,j,Lold)+ &
1419 & fac*(ad_zeta(i,j,Lnew)- &
1420 & ad_zeta(i,j,Lold))
1421#ifdef MASKING
1422 ad_zeta(i,j,Lnew)=ad_zeta(i,j,Lnew)*rmask(i,j)
1423#endif
1424 ad_zeta(i,j,Lold)=ad_zeta(i,j,Lnew)
1425 END DO
1426 END DO
1427
1428#ifndef SOLVE3D
1429!
1430! 2D momentum.
1431!
1432 DO j=JstrR,JendR
1433 DO i=Istr,IendR
1434 ad_ubar(i,j,Lnew)=ad_ubar(i,j,Lold)+ &
1435 & fac*(ad_ubar(i,j,Lnew)- &
1436 & ad_ubar(i,j,Lold))
1437# ifdef MASKING
1438 ad_ubar(i,j,Lnew)=ad_ubar(i,j,Lnew)*umask(i,j)
1439# endif
1440 ad_ubar(i,j,Lold)=ad_ubar(i,j,Lnew)
1441 END DO
1442 END DO
1443 DO j=Jstr,JendR
1444 DO i=IstrR,IendR
1445 ad_vbar(i,j,Lnew)=ad_vbar(i,j,Lold)+ &
1446 & fac*(ad_vbar(i,j,Lnew)- &
1447 & ad_vbar(i,j,Lold))
1448# ifdef MASKING
1449 ad_vbar(i,j,Lnew)=ad_vbar(i,j,Lnew)*vmask(i,j)
1450# endif
1451 ad_vbar(i,j,Lold)=ad_vbar(i,j,Lnew)
1452 END DO
1453 END DO
1454#endif
1455
1456#ifdef ADJUST_WSTRESS
1457!
1458! Surface momentum stress.
1459!
1460 DO j=JstrR,JendR
1461 DO i=Istr,IendR
1462 ad_ustr(i,j,Lnew)=ad_ustr(i,j,Lold)+ &
1463 & fac*(ad_ustr(i,j,Lnew)- &
1464 & ad_ustr(i,j,Lold))
1465# ifdef MASKING
1466 ad_ustr(i,j,Lnew)=ad_ustr(i,j,Lnew)*umask(i,j)
1467# endif
1468 ad_ustr(i,j,Lold)=ad_ustr(i,j,Lnew)
1469 END DO
1470 END DO
1471 DO j=Jstr,JendR
1472 DO i=IstrR,IendR
1473 ad_vstr(i,j,Lnew)=ad_vstr(i,j,Lold)+ &
1474 & fac*(ad_vstr(i,j,Lnew)- &
1475 & ad_vstr(i,j,Lold))
1476# ifdef MASKING
1477 ad_vstr(i,j,Lnew)=ad_vstr(i,j,Lnew)*vmask(i,j)
1478# endif
1479 ad_vstr(i,j,Lold)=ad_vstr(i,j,Lnew)
1480 END DO
1481 END DO
1482#endif
1483#ifdef SOLVE3D
1484!
1485! 3D state variables.
1486!
1487 DO k=1,N(ng)
1488 DO j=JstrR,JendR
1489 DO i=Istr,IendR
1490 ad_u(i,j,k,Lnew)=ad_u(i,j,k,Lold)+ &
1491 & fac*(ad_u(i,j,k,Lnew)- &
1492 & ad_u(i,j,k,Lold))
1493# ifdef MASKING
1494 ad_u(i,j,k,Lnew)=ad_u(i,j,k,Lnew)*umask(i,j)
1495# endif
1496 ad_u(i,j,k,Lold)=ad_u(i,j,k,Lnew)
1497 END DO
1498 END DO
1499 DO j=Jstr,JendR
1500 DO i=IstrR,IendR
1501 ad_v(i,j,k,Lnew)=ad_v(i,j,k,Lold)+ &
1502 & fac*(ad_v(i,j,k,Lnew)- &
1503 & ad_v(i,j,k,Lold))
1504# ifdef MASKING
1505 ad_v(i,j,k,Lnew)=ad_v(i,j,k,Lnew)*vmask(i,j)
1506# endif
1507 ad_v(i,j,k,Lold)=ad_v(i,j,k,Lnew)
1508 END DO
1509 END DO
1510 END DO
1511!
1512! Tracers.
1513!
1514 DO itrc=1,NT(ng)
1515 DO k=1,N(ng)
1516 DO j=JstrR,JendR
1517 DO i=IstrR,IendR
1518 ad_t(i,j,k,Lnew,itrc)=ad_t(i,j,k,Lold,itrc)+ &
1519 & fac*(ad_t(i,j,k,Lnew,itrc)- &
1520 & ad_t(i,j,k,Lold,itrc))
1521# ifdef MASKING
1522 ad_t(i,j,k,Lnew,itrc)=ad_t(i,j,k,Lnew,itrc)*rmask(i,j)
1523# endif
1524 ad_t(i,j,k,Lold,itrc)=ad_t(i,j,k,Lnew,itrc)
1525 END DO
1526 END DO
1527 END DO
1528 END DO
1529
1530# ifdef ADJUST_STFLUX
1531!
1532! Tracers.
1533!
1534 DO itrc=1,NT(ng)
1535 DO j=JstrR,JendR
1536 DO i=IstrR,IendR
1537 ad_tflux(i,j,Lnew,itrc)=ad_tflux(i,j,Lold,itrc)+ &
1538 & fac*(ad_tflux(i,j,Lnew,itrc)- &
1539 & ad_tflux(i,j,Lold,itrc))
1540# ifdef MASKING
1541 ad_tflux(i,j,Lnew,itrc)=ad_tflux(i,j,Lnew,itrc)*rmask(i,j)
1542# endif
1543 ad_tflux(i,j,Lold,itrc)=ad_tflux(i,j,Lnew,itrc)
1544 END DO
1545 END DO
1546 END DO
1547# endif
1548#endif
1549
1550 RETURN
1551 END SUBROUTINE ad_new_state
1552!
1553!***********************************************************************
1554 SUBROUTINE orthogonalize (ng, model, Istr, Iend, Jstr, Jend, &
1555 & LBi, UBi, LBj, UBj, &
1556 & Lold, Lnew, Lwrk, &
1557 & innLoop, outLoop, &
1558#ifdef MASKING
1559 & rmask, umask, vmask, &
1560#endif
1561#ifdef ADJUST_WSTRESS
1562 & nl_ustr, nl_vstr, &
1563#endif
1564#ifdef SOLVE3D
1565# ifdef ADJUST_STFLUX
1566 & nl_tflux, &
1567# endif
1568 & nl_t, nl_u, nl_v, &
1569#else
1570 & nl_ubar, nl_vbar, &
1571#endif
1572 & nl_zeta, &
1573#ifdef ADJUST_WSTRESS
1574 & tl_ustr, tl_vstr, &
1575#endif
1576#ifdef SOLVE3D
1577# ifdef ADJUST_STFLUX
1578 & tl_tflux, &
1579# endif
1580 & tl_t, tl_u, tl_v, &
1581#else
1582 & tl_ubar, tl_vbar, &
1583#endif
1584 & tl_zeta, &
1585#ifdef ADJUST_WSTRESS
1586 & ad_ustr, ad_vstr, &
1587#endif
1588#ifdef SOLVE3D
1589# ifdef ADJUST_STFLUX
1590 & ad_tflux, &
1591# endif
1592 & ad_t, ad_u, ad_v, &
1593#else
1594 & ad_ubar, ad_vbar, &
1595#endif
1596 & ad_zeta)
1597!***********************************************************************
1598!
1599 USE mod_param
1600 USE mod_parallel
1601 USE mod_fourdvar
1602 USE mod_iounits
1603 USE mod_ncparam
1604 USE mod_scalars
1605!
1606! Imported variable declarations.
1607!
1608 integer, intent(in) :: ng, model, Iend, Istr, Jend, Jstr
1609 integer, intent(in) :: LBi, UBi, LBj, UBj
1610 integer, intent(in) :: Lold, Lnew, Lwrk
1611 integer, intent(in) :: innLoop, outLoop
1612!
1613#ifdef ASSUMED_SHAPE
1614# ifdef MASKING
1615 real(r8), intent(in) :: rmask(LBi:,LBj:)
1616 real(r8), intent(in) :: umask(LBi:,LBj:)
1617 real(r8), intent(in) :: vmask(LBi:,LBj:)
1618# endif
1619# ifdef ADJUST_WSTRESS
1620 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:)
1621 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:)
1622# endif
1623# ifdef SOLVE3D
1624# ifdef ADJUST_STFLUX
1625 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:)
1626# endif
1627 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
1628 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
1629 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
1630# else
1631 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
1632 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
1633# endif
1634 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
1635# ifdef ADJUST_WSTRESS
1636 real(r8), intent(inout) :: nl_ustr(LBi:,LBj:,:)
1637 real(r8), intent(inout) :: nl_vstr(LBi:,LBj:,:)
1638# endif
1639# ifdef SOLVE3D
1640# ifdef ADJUST_STFLUX
1641 real(r8), intent(inout) :: nl_tflux(LBi:,LBj:,:,:)
1642# endif
1643 real(r8), intent(inout) :: nl_t(LBi:,LBj:,:,:,:)
1644 real(r8), intent(inout) :: nl_u(LBi:,LBj:,:,:)
1645 real(r8), intent(inout) :: nl_v(LBi:,LBj:,:,:)
1646# else
1647 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
1648 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
1649# endif
1650 real(r8), intent(inout) :: nl_zeta(LBi:,LBj:,:)
1651# ifdef ADJUST_WSTRESS
1652 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:)
1653 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:)
1654# endif
1655# ifdef SOLVE3D
1656# ifdef ADJUST_STFLUX
1657 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:)
1658# endif
1659 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
1660 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
1661 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
1662# else
1663 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
1664 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
1665# endif
1666 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
1667
1668#else
1669
1670# ifdef MASKING
1671 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
1672 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
1673 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
1674# endif
1675# ifdef ADJUST_WSTRESS
1676 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,2)
1677 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,2)
1678# endif
1679# ifdef SOLVE3D
1680# ifdef ADJUST_STFLUX
1681 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
1682# endif
1683 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1684 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
1685 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
1686# else
1687 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,3)
1688 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,3)
1689# endif
1690 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,3)
1691# ifdef ADJUST_WSTRESS
1692 real(r8), intent(inout) :: nl_ustr(LBi:UBi,LBj:UBj,2)
1693 real(r8), intent(inout) :: nl_vstr(LBi:UBi,LBj:UBj,2)
1694# endif
1695# ifdef SOLVE3D
1696# ifdef ADJUST_STFLUX
1697 real(r8), intent(inout) :: nl_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
1698# endif
1699 real(r8), intent(inout) :: nl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1700 real(r8), intent(inout) :: nl_u(LBi:UBi,LBj:UBj,N(ng),2)
1701 real(r8), intent(inout) :: nl_v(LBi:UBi,LBj:UBj,N(ng),2)
1702# else
1703 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,3)
1704 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,3)
1705# endif
1706 real(r8), intent(inout) :: nl_zeta(LBi:UBi,LBj:UBj,3)
1707# ifdef ADJUST_WSTRESS
1708 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,2)
1709 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,2)
1710# endif
1711# ifdef SOLVE3D
1712# ifdef ADJUST_STFLUX
1713 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
1714# endif
1715 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
1716 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
1717 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
1718# else
1719 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,3)
1720 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,3)
1721# endif
1722 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,3)
1723#endif
1724!
1725! Local variable declarations.
1726!
1727 integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
1728 integer :: i, j, lstr, rec, Lscale
1729#ifdef SOLVE3D
1730 integer :: itrc, k
1731#endif
1732 integer :: L1 = 1
1733 integer :: L2 = 2
1734
1735 real(r8) :: fac, fac1, fac2
1736
1737 real(r8), dimension(0:NstateVar(ng)) :: dot
1738 real(r8), dimension(0:Ninner) :: DotProd, dot_new, dot_old
1739
1740 character (len=80) :: ncname
1741
1742#include "set_bounds.h"
1743!
1744!-----------------------------------------------------------------------
1745! Orthogonalize current gradient, G(k+1), against all previous
1746! gradients (reverse order) using Gramm-Schmidt procedure.
1747!-----------------------------------------------------------------------
1748!
1749! We can overwrite adjoint arrays at index Lnew each time around the
1750! the following loop because the preceding gradient vectors that we
1751! read are orthogonal to each other. The reversed order of the loop
1752! is important for the Lanczos vector calculations.
1753!
1754 DO rec=innLoop,1,-1
1755!
1756! Determine adjoint file to process.
1757!
1758 IF (ndefADJ(ng).gt.0) THEN
1759 lstr=LEN_TRIM(ADJbase(ng))
1760 WRITE (ncname,10) ADJbase(ng)(1:lstr-3), rec
1761 10 FORMAT (a,'_',i3.3,'.nc')
1762 ELSE
1763 ncname=ADJname(ng)
1764 END IF
1765!
1766! Read in each previous gradient state solutions, G(0) to G(k), and
1767! compute its associated dot against current G(k+1). Each gradient
1768! solution is loaded NLM (index L2, if preconditioning) or
1769! TLM (index Lwrk, if not preconditioning) state arrays.
1770!
1771 IF (Lprecond) THEN
1772 CALL read_state (ng, model, Istr, Iend, Jstr, Jend, &
1773 & LBi, UBi, LBj, UBj, &
1774 & L2, rec, &
1775 & ndefADJ(ng), ncADJid(ng), ncname, &
1776#ifdef MASKING
1777 & rmask, umask, vmask, &
1778#endif
1779#ifdef ADJUST_WSTRESS
1780 & nl_ustr, nl_vstr, &
1781#endif
1782#ifdef SOLVE3D
1783# ifdef ADJUST_STFLUX
1784 & nl_tflux, &
1785# endif
1786 & nl_t, nl_u, nl_v, &
1787#else
1788 & nl_ubar, nl_vbar, &
1789#endif
1790 & nl_zeta)
1791!
1792 ELSE
1793 CALL read_state (ng, model, Istr, Iend, Jstr, Jend, &
1794 & LBi, UBi, LBj, UBj, &
1795 & Lwrk, rec, &
1796 & ndefADJ(ng), ncADJid(ng), ncname, &
1797#ifdef MASKING
1798 & rmask, umask, vmask, &
1799#endif
1800#ifdef ADJUST_WSTRESS
1801 & tl_ustr, tl_vstr, &
1802#endif
1803#ifdef SOLVE3D
1804# ifdef ADJUST_STFLUX
1805 & tl_tflux, &
1806# endif
1807 & tl_t, tl_u, tl_v, &
1808#else
1809 & tl_ubar, tl_vbar, &
1810#endif
1811 & tl_zeta)
1812 END IF
1813!
1814! If preconditioning, compute H^-1 * G(rec) and store it TLM state
1815! arrays (index Lwrk).
1816!
1817 IF (Lprecond) THEN
1818 Lscale=-1
1819 CALL precond (ng, model, Istr, Iend, Jstr, Jend, &
1820 & LBi, UBi, LBj, UBj, &
1821 & NstateVar(ng), L2, Lwrk, Lscale, &
1822 & nConvRitz, Ritz, &
1823#ifdef MASKING
1824 & rmask, umask, vmask, &
1825#endif
1826#ifdef ADJUST_WSTRESS
1827 & nl_ustr, nl_ustr, tl_ustr, &
1828 & nl_vstr, nl_vstr, tl_vstr, &
1829#endif
1830#ifdef SOLVE3D
1831# ifdef ADJUST_STFLUX
1832 & nl_tflux, nl_tflux, tl_tflux, &
1833# endif
1834 & nl_t, nl_t, tl_t, &
1835 & nl_u, nl_u, tl_u, &
1836 & nl_v, nl_v, tl_v, &
1837#else
1838 & nl_ubar, nl_ubar, tl_ubar, &
1839 & nl_vbar, nl_vbar, tl_vbar, &
1840#endif
1841 & nl_zeta, nl_zeta, tl_zeta)
1842 END IF
1843!
1844! If preconditioning, compute dot product <G(k+1), H^-1 G(rec)>.
1845! Otherwise, compute <G(k+1), G(rec)>. Recall that the TLM
1846! (index Lwrk) contains either H^-1 G(rec) or G(rec).
1847!
1848 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
1849 & LBi, UBi, LBj, UBj, &
1850 & NstateVar(ng), dot(0:), &
1851#ifdef MASKING
1852 & rmask, umask, vmask, &
1853#endif
1854#ifdef ADJUST_WSTRESS
1855 & ad_ustr(:,:,Lnew), tl_ustr(:,:,Lwrk), &
1856 & ad_vstr(:,:,Lnew), tl_vstr(:,:,Lwrk), &
1857#endif
1858#ifdef SOLVE3D
1859# ifdef ADJUST_STFLUX
1860 & ad_tflux(:,:,Lnew,:), tl_tflux(:,:,Lwrk,:), &
1861# endif
1862 & ad_t(:,:,:,Lnew,:), tl_t(:,:,:,Lwrk,:), &
1863 & ad_u(:,:,:,Lnew), tl_u(:,:,:,Lwrk), &
1864 & ad_v(:,:,:,Lnew), tl_v(:,:,:,Lwrk), &
1865#else
1866 & ad_ubar(:,:,Lnew), tl_ubar(:,:,Lwrk), &
1867 & ad_vbar(:,:,Lnew), tl_vbar(:,:,Lwrk), &
1868#endif
1869 & ad_zeta(:,:,Lnew), tl_zeta(:,:,Lwrk))
1870 dot_new(rec)=dot(0)
1871!
1872! If preconditioning, compute dot product <G(rec), H^-1 * G(rec)>.
1873!
1874 IF (Lprecond) THEN
1875 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
1876 & LBi, UBi, LBj, UBj, &
1877 & NstateVar(ng), dot(0:), &
1878#ifdef MASKING
1879 & rmask, umask, vmask, &
1880#endif
1881#ifdef ADJUST_WSTRESS
1882 & nl_ustr(:,:,L2), tl_ustr(:,:,Lwrk), &
1883 & nl_vstr(:,:,L2), tl_vstr(:,:,Lwrk), &
1884#endif
1885#ifdef SOLVE3D
1886# ifdef ADJUST_STFLUX
1887 & nl_tflux(:,:,L2,:), tl_tflux(:,:,Lwrk,:), &
1888# endif
1889 & nl_t(:,:,:,L2,:), tl_t(:,:,:,Lwrk,:), &
1890 & nl_u(:,:,:,L2), tl_u(:,:,:,Lwrk), &
1891 & nl_v(:,:,:,L2), tl_v(:,:,:,Lwrk), &
1892#else
1893 & nl_ubar(:,:,L2), tl_ubar(:,:,Lwrk), &
1894 & nl_vbar(:,:,L2), tl_vbar(:,:,Lwrk), &
1895#endif
1896 & nl_zeta(:,:,L2), tl_zeta(:,:,Lwrk))
1897!
1898! Otherwise, compute dot product <G(rec), G(rec)>.
1899!
1900 ELSE
1901 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
1902 & LBi, UBi, LBj, UBj, &
1903 & NstateVar(ng), dot(0:), &
1904#ifdef MASKING
1905 & rmask, umask, vmask, &
1906#endif
1907#ifdef ADJUST_WSTRESS
1908 & tl_ustr(:,:,Lwrk), tl_ustr(:,:,Lwrk), &
1909 & tl_vstr(:,:,Lwrk), tl_vstr(:,:,Lwrk), &
1910#endif
1911#ifdef SOLVE3D
1912# ifdef ADJUST_STFLUX
1913 & tl_tflux(:,:,Lwrk,:),tl_tflux(:,:,Lwrk,:),&
1914# endif
1915 & tl_t(:,:,:,Lwrk,:), tl_t(:,:,:,Lwrk,:), &
1916 & tl_u(:,:,:,Lwrk), tl_u(:,:,:,Lwrk), &
1917 & tl_v(:,:,:,Lwrk), tl_v(:,:,:,Lwrk), &
1918#else
1919 & tl_ubar(:,:,Lwrk), tl_ubar(:,:,Lwrk), &
1920 & tl_vbar(:,:,Lwrk), tl_vbar(:,:,Lwrk), &
1921#endif
1922 & tl_zeta(:,:,Lwrk), tl_zeta(:,:,Lwrk))
1923 END IF
1924 dot_old(rec)=dot(0)
1925!
1926! Compute Gramm-Schmidt scaling coefficient.
1927!
1928 DotProd(rec)=dot_new(rec)/dot_old(rec)
1929
1930 fac1=1.0_r8
1931 fac2=-DotProd(rec)
1932!
1933! If preconditioning, perform Gramm-Schmidt orthonormalization as:
1934!
1935! ad_var(Lnew) = fac1 * ad_var(Lnew) + fac2 * nl_var(L2)
1936!
1937 IF (Lprecond) THEN
1938 CALL state_addition (ng, Istr, Iend, Jstr, Jend, &
1939 & LBi, UBi, LBj, UBj, &
1940 & Lnew, L2, Lnew, fac1, fac2, &
1941#ifdef MASKING
1942 & rmask, umask, vmask, &
1943#endif
1944#ifdef ADJUST_WSTRESS
1945 & ad_ustr, nl_ustr, &
1946 & ad_vstr, nl_vstr, &
1947#endif
1948#ifdef SOLVE3D
1949# ifdef ADJUST_STFLUX
1950 & ad_tflux, nl_tflux, &
1951# endif
1952 & ad_t, nl_t, &
1953 & ad_u, nl_u, &
1954 & ad_v, nl_v, &
1955#else
1956 & ad_ubar, nl_ubar, &
1957 & ad_vbar, nl_vbar, &
1958#endif
1959 & ad_zeta, nl_zeta)
1960!
1961! If not preconditioning, perform Gramm-Schmidt orthonormalization as:
1962!
1963! ad_var(Lnew) = fac1 * ad_var(Lnew) + fac2 * tl_var(Lwrk)
1964!
1965 ELSE
1966 CALL state_addition (ng, Istr, Iend, Jstr, Jend, &
1967 & LBi, UBi, LBj, UBj, &
1968 & Lnew, Lwrk, Lnew, fac1, fac2, &
1969#ifdef MASKING
1970 & rmask, umask, vmask, &
1971#endif
1972#ifdef ADJUST_WSTRESS
1973 & ad_ustr, tl_vstr, &
1974#endif
1975#ifdef SOLVE3D
1976# ifdef ADJUST_STFLUX
1977 & ad_tflux, tl_tflux, &
1978# endif
1979 & ad_t, tl_t, &
1980 & ad_u, tl_u, &
1981 & ad_v, tl_v, &
1982#else
1983 & ad_ubar, tl_ubar, &
1984 & ad_vbar, tl_vbar, &
1985#endif
1986 & ad_zeta, tl_zeta)
1987 END IF
1988 END DO
1989
1990#ifdef TEST_ORTHOGONALIZATION
1991!
1992!-----------------------------------------------------------------------
1993! Test orthogonal properties of the new gradient.
1994!-----------------------------------------------------------------------
1995!
1996 DO rec=innLoop,1,-1
1997!
1998! Determine adjoint file to process.
1999!
2000 IF (ndefADJ(ng).gt.0) THEN
2001 lstr=LEN_TRIM(ADJbase(ng))
2002 WRITE (ncname,10) ADJbase(ng)(1:lstr-3), rec
2003 ELSE
2004 ncname=ADJname(ng)
2005 END IF
2006!
2007! Read in each previous gradient state solutions, G(0) to G(k), and
2008! compute its associated dot angaint orthogonalized G(k+1). Again,
2009! each gradient solution is loaded into TANGENT LINEAR STATE ARRAYS
2010! at index Lwrk.
2011!
2012 CALL read_state (ng, model, Istr, Iend, Jstr, Jend, &
2013 & LBi, UBi, LBj, UBj, &
2014 & Lwrk, rec, &
2015 & ndefADJ(ng), ncADJid(ng), ncname, &
2016#ifdef MASKING
2017 & rmask, umask, vmask, &
2018#endif
2019#ifdef ADJUST_WSTRESS
2020 & tl_ustr, tl_vstr, &
2021#endif
2022#ifdef SOLVE3D
2023# ifdef ADJUST_STFLUX
2024 & tl_stlfx, &
2025# endif
2026 & tl_t, tl_u, tl_v, &
2027#else
2028 & tl_ubar, tl_vbar, &
2029#endif
2030 & tl_zeta)
2031!
2032 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
2033 & LBi, UBi, LBj, UBj, &
2034 & NstateVar(ng), dot(0:), &
2035#ifdef MASKING
2036 & rmask, umask, vmask, &
2037#endif
2038#ifdef ADJUST_WSTRESS
2039 & ad_ustr(:,:,Lnew), tl_ustr(:,:,Lwrk), &
2040 & ad_vstr(:,:,Lnew), tl_vstr(:,:,Lwrk), &
2041#endif
2042#ifdef SOLVE3D
2043# ifdef ADJUST_STFLUX
2044 & ad_tflux(:,:,Lnew,:), tl_tflux(:,:,Lwrk,:), &
2045# endif
2046 & ad_t(:,:,:,Lnew,:), tl_t(:,:,:,Lwrk,:), &
2047 & ad_u(:,:,:,Lnew), tl_u(:,:,:,Lwrk), &
2048 & ad_v(:,:,:,Lnew), tl_v(:,:,:,Lwrk), &
2049#else
2050 & ad_ubar(:,:,Lnew), tl_ubar(:,:,Lwrk), &
2051 & ad_vbar(:,:,Lnew), tl_vbar(:,:,Lwrk), &
2052#endif
2053 & ad_zeta(:,:,Lnew), tl_zeta(:,:,Lwrk))
2054 dot_new(rec)=dot(0)
2055 END DO
2056!
2057! Report dot products. If everything is working correctly, at the
2058! end of the orthogonalization dot_new(rec) << dot_old(rec).
2059!
2060 IF (Master) THEN
2061 WRITE (stdout,20) outer, inner
2062 DO rec=innLoop,1,-1
2063 WRITE (stdout,30) DotProd(rec), rec-1
2064 END DO
2065 WRITE (stdout,*) ' '
2066 DO rec=innLoop,1,-1
2067 WRITE (stdout,40) innLoop, rec-1, dot_new(rec), &
2068 & rec-1, rec-1, dot_old(rec)
2069 END DO
2070 20 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
2071 & 'Gramm-Schmidt Orthogonalization:',/)
2072 30 FORMAT (12x,'Orthogonalization Factor = ',1p,e19.12,3x, &
2073 & '(Iter=',i3.3,')')
2074 40 FORMAT (2x,'Ortho Test: ', &
2075 & '<G(',i3.3,'),G(',i3.3,')> = ',1p,e15.8,1x, &
2076 & '<G(',i3.3,'),G(',i3.3,')> = ',1p,e15.8)
2077 END IF
2078#endif
2079
2080 RETURN
2081 END SUBROUTINE orthogonalize
2082
2083!
2084!***********************************************************************
2085 SUBROUTINE new_direction (ng, model, Istr, Iend, Jstr, Jend, &
2086 & LBi, UBi, LBj, UBj, &
2087 & Lwrk, Lnew, betaK, &
2088#ifdef MASKING
2089 & rmask, umask, vmask, &
2090#endif
2091#ifdef ADJUST_WSTRESS
2092 & ad_ustr, ad_vstr, &
2093#endif
2094#ifdef SOLVE3D
2095# ifdef ADJUST_STFLUX
2096 & ad_tflux, &
2097# endif
2098 & ad_t, ad_u, ad_v, &
2099#else
2100 & ad_ubar, ad_vbar, &
2101#endif
2102 & ad_zeta, &
2103#ifdef ADJUST_WSTRESS
2104 & d_sustr, d_svstr, &
2105#endif
2106#ifdef SOLVE3D
2107# ifdef ADJUST_STFLUX
2108 & d_stflx, &
2109# endif
2110 & d_t, d_u, d_v, &
2111#else
2112 & d_ubar, d_vbar, &
2113#endif
2114 & d_zeta)
2115!***********************************************************************
2116!
2117 USE mod_param
2118 USE mod_parallel
2119!
2120! Imported variable declarations.
2121!
2122 integer, intent(in) :: ng, model, Iend, Istr, Jend, Jstr
2123 integer, intent(in) :: LBi, UBi, LBj, UBj
2124 integer, intent(in) :: Lwrk, Lnew
2125
2126 real(r8), intent(in) :: betaK
2127!
2128#ifdef ASSUMED_SHAPE
2129# ifdef MASKING
2130 real(r8), intent(in) :: rmask(LBi:,LBj:)
2131 real(r8), intent(in) :: umask(LBi:,LBj:)
2132 real(r8), intent(in) :: vmask(LBi:,LBj:)
2133# endif
2134# ifdef ADJUST_WSTRESS
2135 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:)
2136 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:)
2137# endif
2138# ifdef SOLVE3D
2139# ifdef ADJUST_STFLUX
2140 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:)
2141# endif
2142 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
2143 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
2144 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
2145# else
2146 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
2147 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
2148# endif
2149 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
2150# ifdef ADJUST_WSTRESS
2151 real(r8), intent(inout) :: d_sustr(LBi:,LBj:)
2152 real(r8), intent(inout) :: d_svstr(LBi:,LBj:)
2153# endif
2154# ifdef SOLVE3D
2155# ifdef ADJUST_STFLUX
2156 real(r8), intent(inout) :: d_stflx(LBi:,LBj:,:)
2157# endif
2158 real(r8), intent(inout) :: d_t(LBi:,LBj:,:,:)
2159 real(r8), intent(inout) :: d_u(LBi:,LBj:,:)
2160 real(r8), intent(inout) :: d_v(LBi:,LBj:,:)
2161# else
2162 real(r8), intent(inout) :: d_ubar(LBi:,LBj:)
2163 real(r8), intent(inout) :: d_vbar(LBi:,LBj:)
2164# endif
2165 real(r8), intent(inout) :: d_zeta(LBi:,LBj:)
2166#else
2167# ifdef MASKING
2168 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
2169 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
2170 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
2171# endif
2172# ifdef ADJUST_WSTRESS
2173 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,2)
2174 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,2)
2175# endif
2176# ifdef SOLVE3D
2177# ifdef ADJUST_STFLUX
2178 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
2179# endif
2180 real(r8), intent(inout) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2181 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
2182 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
2183# else
2184 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,3)
2185 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,3)
2186# endif
2187 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,3)
2188# ifdef ADJUST_WSTRESS
2189 real(r8), intent(inout) :: d_sustr(LBi:UBi,LBj:UBj)
2190 real(r8), intent(inout) :: d_svstr(LBi:UBi,LBj:UBj)
2191# endif
2192# ifdef SOLVE3D
2193# ifdef ADJUST_STFLUX
2194 real(r8), intent(inout) :: d_stflx(LBi:UBi,LBj:UBj,NT(ng))
2195# endif
2196 real(r8), intent(inout) :: d_t(LBi:UBi,LBj:UBj,N(ng),NT(ng))
2197 real(r8), intent(inout) :: d_u(LBi:UBi,LBj:UBj,N(ng))
2198 real(r8), intent(inout) :: d_v(LBi:UBi,LBj:UBj,N(ng))
2199# else
2200 real(r8), intent(inout) :: d_ubar(LBi:UBi,LBj:UBj)
2201 real(r8), intent(inout) :: d_vbar(LBi:UBi,LBj:UBj)
2202# endif
2203 real(r8), intent(inout) :: d_zeta(LBi:UBi,LBj:UBj)
2204#endif
2205!
2206! Local variable declarations.
2207!
2208 integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
2209 integer :: i, j
2210#ifdef SOLVE3D
2211 integer :: itrc, k
2212#endif
2213
2214#include "set_bounds.h"
2215!
2216!-----------------------------------------------------------------------
2217! Compute new conjugate descent direction, d(k+1). Notice that the old
2218! descent direction is overwritten. Also the initial value is just
2219! d(0)=-G(0) since betaK=0 when inner=0.
2220!-----------------------------------------------------------------------
2221!
2222! Free-surface.
2223!
2224 DO j=JstrR,JendR
2225 DO i=IstrR,IendR
2226 d_zeta(i,j)=-ad_zeta(i,j,Lnew)+betaK*d_zeta(i,j)
2227#ifdef MASKING
2228 d_zeta(i,j)=d_zeta(i,j)*rmask(i,j)
2229#endif
2230 END DO
2231 END DO
2232
2233#ifndef SOLVE3D
2234!
2235! 2D momentum.
2236!
2237 DO j=JstrR,JendR
2238 DO i=Istr,IendR
2239 d_ubar(i,j)=-ad_ubar(i,j,Lnew)+betaK*d_ubar(i,j)
2240# ifdef MASKING
2241 d_ubar(i,j)=d_ubar(i,j)*umask(i,j)
2242# endif
2243 END DO
2244 END DO
2245 DO j=Jstr,JendR
2246 DO i=IstrR,IendR
2247 d_vbar(i,j)=-ad_vbar(i,j,Lnew)+betaK*d_vbar(i,j)
2248# ifdef MASKING
2249 d_vbar(i,j)=d_vbar(i,j)*vmask(i,j)
2250# endif
2251 END DO
2252 END DO
2253#endif
2254
2255#ifdef ADJUST_WSTRESS
2256!
2257! Surface momentum stress.
2258!
2259 DO j=JstrR,JendR
2260 DO i=Istr,IendR
2261 d_sustr(i,j)=-ad_ustr(i,j,Lnew)+betaK*d_sustr(i,j)
2262# ifdef MASKING
2263 d_sustr(i,j)=d_sustr(i,j)*umask(i,j)
2264# endif
2265 END DO
2266 END DO
2267 DO j=Jstr,JendR
2268 DO i=IstrR,IendR
2269 d_svstr(i,j)=-ad_vstr(i,j,Lnew)+betaK*d_svstr(i,j)
2270# ifdef MASKING
2271 d_svstr(i,j)=d_svstr(i,j)*vmask(i,j)
2272# endif
2273 END DO
2274 END DO
2275#endif
2276
2277#ifdef SOLVE3D
2278!
2279! 3D momentum.
2280!
2281 DO k=1,N(ng)
2282 DO j=JstrR,JendR
2283 DO i=Istr,IendR
2284 d_u(i,j,k)=-ad_u(i,j,k,Lnew)+betaK*d_u(i,j,k)
2285# ifdef MASKING
2286 d_u(i,j,k)=d_u(i,j,k)*umask(i,j)
2287# endif
2288 END DO
2289 END DO
2290 DO j=Jstr,JendR
2291 DO i=IstrR,IendR
2292 d_v(i,j,k)=-ad_v(i,j,k,Lnew)+betaK*d_v(i,j,k)
2293# ifdef MASKING
2294 d_v(i,j,k)=d_v(i,j,k)*vmask(i,j)
2295# endif
2296 END DO
2297 END DO
2298 END DO
2299!
2300! Tracers.
2301!
2302 DO itrc=1,NT(ng)
2303 DO k=1,N(ng)
2304 DO j=JstrR,JendR
2305 DO i=IstrR,IendR
2306 d_t(i,j,k,itrc)=-ad_t(i,j,k,Lnew,itrc)+ &
2307 & betaK*d_t(i,j,k,itrc)
2308# ifdef MASKING
2309 d_t(i,j,k,itrc)=d_t(i,j,k,itrc)*rmask(i,j)
2310# endif
2311 END DO
2312 END DO
2313 END DO
2314 END DO
2315
2316# ifdef ADJUST_STFLUX
2317!
2318! Surface tracers flux.
2319!
2320 DO itrc=1,NT(ng)
2321 DO j=JstrR,JendR
2322 DO i=IstrR,IendR
2323 d_stflx(i,j,itrc)=-ad_tflux(i,j,Lnew,itrc)+ &
2324 & betaK*d_stflx(i,j,itrc)
2325# ifdef MASKING
2326 d_stflx(i,j,itrc)=d_stflx(i,j,itrc)*rmask(i,j)
2327# endif
2328 END DO
2329 END DO
2330 END DO
2331# endif
2332#endif
2333
2334 RETURN
2335 END SUBROUTINE new_direction
2336
2337!
2338!**********************************************************************
2339 SUBROUTINE precond (ng, model, Istr, Iend, Jstr, Jend, &
2340 & LBi, UBi, LBj, UBj, &
2341 & NstateVars, Linp, Lwrk, Lscale, &
2342 & nConvRitz, Ritz, &
2343#ifdef MASKING
2344 & rmask, umask, vmask, &
2345#endif
2346#ifdef ADJUST_WSTRESS
2347 & ad_ustr, nl_ustr, tl_ustr, &
2348 & ad_vstr, nl_vstr, tl_vstr, &
2349#endif
2350#ifdef SOLVE3D
2351# ifdef ADJUST_STFLUX
2352 & ad_tflux, nl_tflux, tl_tflux, &
2353# endif
2354 & ad_t, nl_t, tl_t, &
2355 & ad_u, nl_u, tl_u, &
2356 & ad_v, nl_v, tl_v, &
2357#else
2358 & ad_ubar, nl_ubar, tl_ubar, &
2359 & ad_vbar, nl_vbar, tl_vbar, &
2360#endif
2361 & ad_zeta, nl_zeta, tl_zeta)
2362!***********************************************************************
2363!
2364 USE mod_param
2365 USE mod_parallel
2366 USE mod_ncparam
2367 USE mod_netcdf
2368 USE mod_iounits
2369#ifdef DISTRIBUTE
2370!
2371 USE distribute_mod, ONLY : mp_reduce
2372#endif
2373!
2374 implicit none
2375!
2376! Imported variable declarations.
2377!
2378 integer, intent(in) :: ng, model, Iend, Istr, Jend, Jstr
2379 integer, intent(in) :: LBi, UBi, LBj, UBj
2380 integer, intent(in) :: NstateVars, Linp, Lwrk, Lscale
2381 integer, intent(in) :: nConvRitz
2382!
2383 real(r8), intent(in) :: Ritz(:)
2384!
2385#ifdef ASSUMED_SHAPE
2386# ifdef MASKING
2387 real(r8), intent(in) :: rmask(LBi:,LBj:)
2388 real(r8), intent(in) :: umask(LBi:,LBj:)
2389 real(r8), intent(in) :: vmask(LBi:,LBj:)
2390# endif
2391#ifdef ADJUST_WSTRESS
2392 real(r8), intent(in) :: ad_ustr(LBi:,LBj:,:)
2393 real(r8), intent(in) :: ad_vstr(LBi:,LBj:,:)
2394#endif
2395#ifdef SOLVE3D
2396# ifdef ADJUST_STFLUX
2397 real(r8), intent(in) :: ad_tflux(LBi:,LBj:,:,:)
2398# endif
2399 real(r8), intent(in) :: ad_t(LBi:,LBj:,:,:,:)
2400 real(r8), intent(in) :: ad_u(LBi:,LBj:,:,:)
2401 real(r8), intent(in) :: ad_v(LBi:,LBj:,:,:)
2402#else
2403 real(r8), intent(in) :: ad_ubar(LBi:,LBj:,:)
2404 real(r8), intent(in) :: ad_vbar(LBi:,LBj:,:)
2405#endif
2406 real(r8), intent(in) :: ad_zeta(LBi:,LBj:,:)
2407# ifdef ADJUST_WSTRESS
2408 real(r8), intent(inout) :: nl_ustr(LBi:,LBj:,:)
2409 real(r8), intent(inout) :: nl_vstr(LBi:,LBj:,:)
2410# endif
2411# ifdef SOLVE3D
2412# ifdef ADJUST_STFLUX
2413 real(r8), intent(inout) :: nl_tflux(LBi:,LBj:,:,:)
2414# endif
2415 real(r8), intent(inout) :: nl_t(LBi:,LBj:,:,:,:)
2416 real(r8), intent(inout) :: nl_u(LBi:,LBj:,:,:)
2417 real(r8), intent(inout) :: nl_v(LBi:,LBj:,:,:)
2418# else
2419 real(r8), intent(inout) :: nl_ubar(LBi:,LBj:,:)
2420 real(r8), intent(inout) :: nl_vbar(LBi:,LBj:,:)
2421# endif
2422 real(r8), intent(inout) :: nl_zeta(LBi:,LBj:,:)
2423# ifdef ADJUST_WSTRESS
2424 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:)
2425 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:)
2426# endif
2427# ifdef SOLVE3D
2428# ifdef ADJUST_STFLUX
2429 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:)
2430# endif
2431 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
2432 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
2433 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
2434# else
2435 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
2436 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
2437# endif
2438 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
2439
2440#else
2441
2442# ifdef MASKING
2443 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
2444 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
2445 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
2446# endif
2447#ifdef ADJUST_WSTRESS
2448 real(r8), intent(in) :: ad_sustr(LBi:UBi,LBj:UBj,2)
2449 real(r8), intent(in) :: ad_svstr(LBi:UBi,LBj:UBj,2)
2450#endif
2451#ifdef SOLVE3D
2452# ifdef ADJUST_STFLUX
2453 real(r8), intent(in) :: ad_stflx(LBi:UBi,LBj:UBj,2,NT(ng))
2454# endif
2455 real(r8), intent(in) :: ad_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2456 real(r8), intent(in) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
2457 real(r8), intent(in) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
2458#else
2459 real(r8), intent(in) :: ad_ubar(LBi:UBi,LBj:UBj,3)
2460 real(r8), intent(in) :: ad_vbar(LBi:UBi,LBj:UBj,3)
2461#endif
2462 real(r8), intent(in) :: ad_zeta(LBi:UBi,LBj:UBj,3)
2463# ifdef ADJUST_WSTRESS
2464 real(r8), intent(inout) :: nl_ustr(LBi:UBi,LBj:UBj,2)
2465 real(r8), intent(inout) :: nl_vstr(LBi:UBi,LBj:UBj,2)
2466# endif
2467# ifdef SOLVE3D
2468# ifdef ADJUST_STFLUX
2469 real(r8), intent(inout) :: nl_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
2470# endif
2471 real(r8), intent(inout) :: nl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2472 real(r8), intent(inout) :: nl_u(LBi:UBi,LBj:UBj,N(ng),2)
2473 real(r8), intent(inout) :: nl_v(LBi:UBi,LBj:UBj,N(ng),2)
2474# else
2475 real(r8), intent(inout) :: nl_ubar(LBi:UBi,LBj:UBj,3)
2476 real(r8), intent(inout) :: nl_vbar(LBi:UBi,LBj:UBj,3)
2477# endif
2478 real(r8), intent(inout) :: nl_zeta(LBi:UBi,LBj:UBj,3)
2479# ifdef ADJUST_WSTRESS
2480 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,2)
2481 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,2)
2482# endif
2483# ifdef SOLVE3D
2484# ifdef ADJUST_STFLUX
2485 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
2486# endif
2487 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
2488 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
2489 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
2490# else
2491 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,3)
2492 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,3)
2493# endif
2494 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,3)
2495#endif
2496!
2497! Local variable declarations.
2498!
2499 integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
2500 integer :: NSUB, i, j, L1, L2, nvec
2501#ifdef SOLVE3D
2502 integer :: itrc, k
2503#endif
2504 real(r8) :: cff, fac, fac1, fac2
2505 real(r8), dimension(0:NstateVars) :: Dotprod
2506#ifdef DISTRIBUTE
2507 character (len=3) :: op_handle
2508#endif
2509
2510#include "set_bounds.h"
2511!
2512!-----------------------------------------------------------------------
2513! Apply the preconditioner. The approximated Hessian matrix is computed
2514! from the eigenvectors computed by the Lanczos algorithm which are
2515! stored in HSSname NetCDF file.
2516!-----------------------------------------------------------------------
2517!
2518! Copy ad_var(Linp) into tl_var(Lwrk)
2519!
2520 CALL state_copy (ng, Istr, Iend, Jstr, Jend, &
2521 & LBi, UBi, LBj, UBj, &
2522 & Linp, Lwrk, &
2523#ifdef ADJUST_WSTRESS
2524 & tl_ustr, ad_ustr, &
2525 & tl_vstr, ad_vstr, &
2526#endif
2527#ifdef SOLVE3D
2528# ifdef ADJUST_STFLUX
2529 & tl_tflux, ad_tflux, &
2530# endif
2531 & tl_t, ad_t, &
2532 & tl_u, ad_u, &
2533 & tl_v, ad_v, &
2534#else
2535 & tl_ubar, ad_ubar, &
2536 & tl_vbar, ad_vbar, &
2537#endif
2538 & tl_zeta, ad_zeta)
2539!
2540! Read the converged Hessian eigenvectors into NLM state array,
2541! index L1.
2542!
2543 DO nvec=1,nConvRitz
2544 L1=1
2545 L2=2
2546 CALL read_state (ng, model, Istr, Iend, Jstr, Jend, &
2547 & LBi, UBi, LBj, UBj, &
2548 & L1, nvec, &
2549 & 0, ncHSSid(ng), HSSname(ng), &
2550#ifdef MASKING
2551 & rmask, umask, vmask, &
2552#endif
2553#ifdef ADJUST_WSTRESS
2554 & nl_ustr, nl_vstr, &
2555#endif
2556#ifdef SOLVE3D
2557# ifdef ADJUST_STFLUX
2558 & nl_tflux, &
2559# endif
2560 & nl_t, nl_u, nl_v, &
2561#else
2562 & nl_ubar, nl_vbar, &
2563#endif
2564 & nl_zeta)
2565!
2566! Compute dot product between gradient and Hessian eigenvector.
2567!
2568 CALL state_dotprod (ng, model, Istr, Iend, Jstr, Jend, &
2569 & LBi, UBi, LBj, UBj, &
2570 & NstateVars, Dotprod(0:), &
2571#ifdef MASKING
2572 & rmask, umask, vmask, &
2573#endif
2574#ifdef ADJUST_WSTRESS
2575 & ad_ustr(:,:,Linp), nl_ustr(:,:,L1), &
2576 & ad_vstr(:,:,Linp), nl_vstr(:,:,L1), &
2577#endif
2578#ifdef SOLVE3D
2579# ifdef ADJUST_STFLUX
2580 & ad_tflux(:,:,Linp,:), nl_tflux(:,:,L1,:), &
2581# endif
2582 & ad_t(:,:,:,Linp,:), nl_t(:,:,:,L1,:), &
2583 & ad_u(:,:,:,Linp), nl_u(:,:,:,L1), &
2584 & ad_v(:,:,:,Linp), nl_v(:,:,:,L1), &
2585#else
2586 & ad_ubar(:,:,Linp), nl_ubar(:,:,L1), &
2587 & ad_vbar(:,:,Linp), nl_vbar(:,:,L1), &
2588#endif
2589 & ad_zeta(:,:,Linp), nl_zeta(:,:,L1))
2590!
2591! Lscale determines the form of the preconditioner:
2592!
2593! 1= Hessian
2594! -1= Inverse Hessian
2595! 2= Hessian square root
2596! -2= Inverse Hessian square root
2597!
2598! tl_var(Lwrk) = fac1 * tl_var(Lwrk) + fac2 * nl_var(L1)
2599!
2600 fac1=1.0_r8
2601
2602 IF (Lscale.eq.1) THEN
2603 fac2=(Ritz(nvec)-1.0_r8)*Dotprod(0)
2604 ELSE IF (Lscale.eq.-1) THEN
2605 fac2=(1.0_r8/Ritz(nvec)-1.0_r8)*Dotprod(0)
2606 ELSE IF (Lscale.eq.1) THEN
2607 fac2=(SQRT(Ritz(nvec))-1.0_r8)*Dotprod(0)
2608 ELSE IF (Lscale.eq.-1) THEN
2609 fac2=(1.0_r8/SQRT(Ritz(nvec))-1.0_r8)*Dotprod(0)
2610 END IF
2611
2612 CALL state_addition (ng, Istr, Iend, Jstr, Jend, &
2613 & LBi, UBi, LBj, UBj, &
2614 & Lwrk, L1, Lwrk, fac1, fac2, &
2615#ifdef MASKING
2616 & rmask, umask, vmask, &
2617#endif
2618#ifdef ADJUST_WSTRESS
2619 & tl_ustr, nl_ustr, &
2620 & tl_vstr, nl_vstr, &
2621#endif
2622#ifdef SOLVE3D
2623# ifdef ADJUST_STFLUX
2624 & tl_tflux, nl_tflux, &
2625# endif
2626 & tl_t, nl_t, &
2627 & tl_u, nl_u, &
2628 & tl_v, nl_v, &
2629#else
2630 & tl_ubar, nl_ubar, &
2631 & tl_vbar, nl_vbar, &
2632#endif
2633 & tl_zeta, nl_zeta)
2634 END DO
2635
2636 RETURN
2637 END SUBROUTINE precond
2638!
2639!***********************************************************************
2640 SUBROUTINE read_state (ng, model, Istr, Iend, Jstr, Jend, &
2641 & LBi, UBi, LBj, UBj, &
2642 & Lwrk, rec, &
2643 & ndef, ncfileid, ncname, &
2644#ifdef MASKING
2645 & rmask, umask, vmask, &
2646#endif
2647#ifdef ADJUST_WSTRESS
2648 & s_ustr, s_vstr, &
2649#endif
2650#ifdef SOLVE3D
2651# ifdef ADJUST_STFLUX
2652 & s_tflux, &
2653# endif
2654 & s_t, s_u, s_v, &
2655#else
2656 & s_ubar, s_vbar, &
2657#endif
2658 & s_zeta)
2659!***********************************************************************
2660!
2661 USE mod_param
2662 USE mod_parallel
2663 USE mod_iounits
2664 USE mod_ncparam
2665 USE mod_netcdf
2666 USE mod_scalars
2667!
2668! Imported variable declarations.
2669!
2670 integer, intent(in) :: ng, model, Iend, Istr, Jend, Jstr
2671 integer, intent(in) :: LBi, UBi, LBj, UBj
2672 integer, intent(in) :: Lwrk, rec, ndef, ncfileid
2673
2674 character (len=*), intent(in) :: ncname
2675!
2676#ifdef ASSUMED_SHAPE
2677# ifdef MASKING
2678 real(r8), intent(in) :: rmask(LBi:,LBj:)
2679 real(r8), intent(in) :: umask(LBi:,LBj:)
2680 real(r8), intent(in) :: vmask(LBi:,LBj:)
2681# endif
2682# ifdef ADJUST_WSTRESS
2683 real(r8), intent(inout) :: s_ustr(LBi:,LBj:,:)
2684 real(r8), intent(inout) :: s_vstr(LBi:,LBj:,:)
2685# endif
2686# ifdef SOLVE3D
2687# ifdef ADJUST_STFLUX
2688 real(r8), intent(inout) :: s_tflux(LBi:,LBj:,:,:)
2689# endif
2690 real(r8), intent(inout) :: s_t(LBi:,LBj:,:,:,:)
2691 real(r8), intent(inout) :: s_u(LBi:,LBj:,:,:)
2692 real(r8), intent(inout) :: s_v(LBi:,LBj:,:,:)
2693# else
2694 real(r8), intent(inout) :: s_ubar(LBi:,LBj:,:)
2695 real(r8), intent(inout) :: s_vbar(LBi:,LBj:,:)
2696# endif
2697 real(r8), intent(inout) :: s_zeta(LBi:,LBj:,:)
2698#else
2699# ifdef MASKING
2700 real(r8), intent(in) :: rmask(LBi:UBi,LBj:UBj)
2701 real(r8), intent(in) :: umask(LBi:UBi,LBj:UBj)
2702 real(r8), intent(in) :: vmask(LBi:UBi,LBj:UBj)
2703# endif
2704# ifdef ADJUST_WSTRESS
2705 real(r8), intent(inout) :: s_ubar(LBi:UBi,LBj:UBj,2)
2706 real(r8), intent(inout) :: s_vbar(LBi:UBi,LBj:UBj,2)
2707# endif
2708# ifdef SOLVE3D
2709# ifdef ADJUST_STFLUX
2710 real(r8), intent(inout) :: s_tflux(LBi:UBi,LBj:UBj,2,NT(ng))
2711# endif
2712 real(r8), intent(inout) :: s_t(LBi:UBi,LBj:UBj,N(ng),2,NT(ng))
2713 real(r8), intent(inout) :: s_u(LBi:UBi,LBj:UBj,N(ng),2)
2714 real(r8), intent(inout) :: s_v(LBi:UBi,LBj:UBj,N(ng),2)
2715# else
2716 real(r8), intent(inout) :: s_ubar(LBi:UBi,LBj:UBj,3)
2717 real(r8), intent(inout) :: s_vbar(LBi:UBi,LBj:UBj,3)
2718# endif
2719 real(r8), intent(inout) :: s_zeta(LBi:UBi,LBj:UBj,3)
2720#endif
2721!
2722! Local variable declarations.
2723!
2724 integer :: IstrR, IendR, JstrR, JendR, IstrU, JstrV
2725 integer :: i, j
2726#ifdef SOLVE3D
2727 integer :: itrc, k
2728#endif
2729 integer :: gtype, ncid, status
2730 integer, dimension(NV) :: vid
2731 integer, dimension(4) :: Vsize
2732
2733 integer :: nf_fread2d
2734#ifdef SOLVE3D
2735 integer :: nf_fread3d
2736#endif
2737
2738 real(r8) :: Fmin, Fmax, scale
2739
2740#include "set_bounds.h"
2741!
2742!-----------------------------------------------------------------------
2743! Read in requested model state record. Load data into state array
2744! index Lwrk.
2745!-----------------------------------------------------------------------
2746!
2747! Determine file and variables ids.
2748!
2749 IF (ndef.gt.0) THEN
2750 IF (InpThread) THEN
2751 status=nf_open(TRIM(ncname), nf_nowrite, ncid)
2752 IF (status.ne.nf_noerr) THEN
2753 WRITE (stdout,10) TRIM(ncname)
2754 exit_flag=2
2755 ioerror=status
2756 RETURN
2757 END IF
2758 END IF
2759 ELSE
2760 ncid=ncfileid
2761 END IF
2762#ifndef SOLVE3D
2763 status=nf_inq_varid(ncid, TRIM(Vname(1,idUbar)), vid(idUbar))
2764 status=nf_inq_varid(ncid, TRIM(Vname(1,idVbar)), vid(idVbar))
2765#endif
2766 status=nf_inq_varid(ncid, TRIM(Vname(1,idFsur)), vid(idFsur))
2767#ifdef ADJUST_WSTRESS
2768 status=nf_inq_varid(ncid, TRIM(Vname(1,idUsms)), vid(idUsms))
2769 status=nf_inq_varid(ncid, TRIM(Vname(1,idVsms)), vid(idVsms))
2770#endif
2771#ifdef SOLVE3D
2772 status=nf_inq_varid(ncid, TRIM(Vname(1,idUvel)), vid(idUvel))
2773 status=nf_inq_varid(ncid, TRIM(Vname(1,idVvel)), vid(idVvel))
2774 DO itrc=1,NT(ng)
2775 status=nf_inq_varid(ncid, TRIM(Vname(1,idTvar(itrc))), &
2776 & vid(idTvar(itrc)))
2777# ifdef ADJUST_STFLUX
2778 status=nf_inq_varid(ncid, TRIM(Vname(1,idTsur(itrc))), &
2779 & vid(idTsur(itrc)))
2780# endif
2781 END DO
2782#endif
2783 DO i=1,4
2784 Vsize(i)=0
2785 END DO
2786 scale=1.0_r8
2787!
2788! Read in free-surface.
2789!
2790 gtype=r2dvar
2791 status=nf_fread2d(ng, iTLM, ncid, vid(idFsur), rec, gtype, &
2792 & Vsize, LBi, UBi, LBj, UBj, &
2793 & scale, Fmin, Fmax, &
2794#ifdef MASKING
2795 & rmask(LBi,LBj), &
2796#endif
2797 & s_zeta(LBi,LBj,Lwrk))
2798 IF (status.ne.nf_noerr) THEN
2799 IF (Master) THEN
2800 WRITE (stdout,20) TRIM(Vname(1,idFsur)), rec, TRIM(ncname)
2801 END IF
2802 exit_flag=3
2803 ioerror=status
2804 RETURN
2805 END IF
2806
2807#ifndef SOLVE3D
2808!
2809! Read in 2D momentum.
2810!
2811 gtype=u2dvar
2812 status=nf_fread2d(ng, iTLM, ncid, vid(idUbar), rec, gtype, &
2813 & Vsize, LBi, UBi, LBj, UBj, &
2814 & scale, Fmin, Fmax, &
2815# ifdef MASKING
2816 & umask(LBi,LBj), &
2817# endif
2818 & s_ubar(LBi,LBj,Lwrk))
2819 IF (status.ne.nf_noerr) THEN
2820 IF (Master) THEN
2821 WRITE (stdout,20) TRIM(Vname(1,idUbar)), rec, TRIM(ncname)
2822 END IF
2823 exit_flag=3
2824 ioerror=status
2825 RETURN
2826 END IF
2827
2828 gtype=v2dvar
2829 status=nf_fread2d(ng, iTLM, ncid, vid(idVbar), rec, gtype, &
2830 & Vsize, LBi, UBi, LBj, UBj, &
2831 & scale, Fmin, Fmax, &
2832# ifdef MASKING
2833 & vmask(LBi,LBj), &
2834# endif
2835 & s_vbar(LBi,LBj,Lwrk))
2836 IF (status.ne.nf_noerr) THEN
2837 IF (Master) THEN
2838 WRITE (stdout,20) TRIM(Vname(1,idVbar)), rec, TRIM(ncname)
2839 END IF
2840 exit_flag=3
2841 ioerror=status
2842 RETURN
2843 END IF
2844#endif
2845
2846#ifdef ADJUST_WSTRESS
2847!
2848! Read surface momentum stress.
2849!
2850 gtype=u2dvar
2851 status=nf_fread2d(ng, iTLM, ncid, vid(idUsms), rec, gtype, &
2852 & Vsize, LBi, UBi, LBj, UBj, &
2853 & scale, Fmin, Fmax, &
2854# ifdef MASKING
2855 & umask(LBi,LBj), &
2856# endif
2857 & s_ustr(LBi,LBj,Lwrk))
2858 IF (status.ne.nf_noerr) THEN
2859 IF (Master) THEN
2860 WRITE (stdout,20) TRIM(Vname(1,idUsms)), rec, TRIM(ncname)
2861 END IF
2862 exit_flag=3
2863 ioerror=status
2864 RETURN
2865 END IF
2866
2867 gtype=v2dvar
2868 status=nf_fread2d(ng, iTLM, ncid, vid(idVsms), rec, gtype, &
2869 & Vsize, LBi, UBi, LBj, UBj, &
2870 & scale, Fmin, Fmax, &
2871# ifdef MASKING
2872 & vmask(LBi,LBj), &
2873# endif
2874 & s_vstr(LBi,LBj,Lwrk))
2875 IF (status.ne.nf_noerr) THEN
2876 IF (Master) THEN
2877 WRITE (stdout,20) TRIM(Vname(1,idVsms)), rec, TRIM(ncname)
2878 END IF
2879 exit_flag=3
2880 ioerror=status
2881 RETURN
2882 END IF
2883#endif
2884
2885#ifdef SOLVE3D
2886!
2887! Read in 3D momentum.
2888!
2889 gtype=u3dvar
2890 status=nf_fread3d(ng, iTLM, ncid, vid(idUvel), rec, gtype, &
2891 & Vsize, LBi, UBi, LBj, UBj, 1, N(ng), &
2892 & scale, Fmin, Fmax, &
2893# ifdef MASKING
2894 & umask(LBi,LBj), &
2895# endif
2896 & s_u(LBi,LBj,1,Lwrk))
2897 IF (status.ne.nf_noerr) THEN
2898 IF (Master) THEN
2899 WRITE (stdout,20) TRIM(Vname(1,idUvel)), rec, TRIM(ncname)
2900 END IF
2901 exit_flag=3
2902 ioerror=status
2903 RETURN
2904 END IF
2905
2906 gtype=v3dvar
2907 status=nf_fread3d(ng, iTLM, ncid, vid(idVvel), rec, gtype, &
2908 & Vsize, LBi, UBi, LBj, UBj, 1, N(ng), &
2909 & scale, Fmin, Fmax, &
2910# ifdef MASKING
2911 & vmask(LBi,LBj), &
2912# endif
2913 & s_v(LBi,LBj,1,Lwrk))
2914 IF (status.ne.nf_noerr) THEN
2915 IF (Master) THEN
2916 WRITE (stdout,20) TRIM(Vname(1,idVvel)), rec, TRIM(ncname)
2917 END IF
2918 exit_flag=3
2919 ioerror=status
2920 RETURN
2921 END IF
2922!
2923! Read in tracers.
2924!
2925 gtype=r3dvar
2926 DO itrc=1,NT(ng)
2927 status=nf_fread3d(ng, iTLM, ncid, vid(idTvar(itrc)), rec, &
2928 & gtype, Vsize, LBi, UBi, LBj, UBj, 1, N(ng), &
2929 & scale, Fmin, Fmax, &
2930# ifdef MASKING
2931 & rmask(LBi,LBj), &
2932# endif
2933 & s_t(LBi,LBj,1,Lwrk,itrc))
2934 IF (status.ne.nf_noerr) THEN
2935 IF (Master) THEN
2936 WRITE (stdout,20) TRIM(Vname(1,idTvar(itrc))), rec, &
2937 & TRIM(ncname)
2938 END IF
2939 exit_flag=3
2940 ioerror=status
2941 RETURN
2942 END IF
2943 END DO
2944
2945# ifdef ADJUST_STFLUX
2946!
2947! Read in surface tracers flux.
2948!
2949 gtype=r2dvar
2950 DO itrc=1,NT(ng)
2951 status=nf_fread2d(ng, iTLM, ncid, vid(idTsur(itrc)), rec, &
2952 & gtype, Vsize, LBi, UBi, LBj, UBj, &
2953 & scale, Fmin, Fmax, &
2954# ifdef MASKING
2955 & rmask(LBi,LBj), &
2956# endif
2957 & s_tflux(LBi,LBj,Lwrk,itrc))
2958 IF (status.ne.nf_noerr) THEN
2959 IF (Master) THEN
2960 WRITE (stdout,20) TRIM(Vname(1,idTsur(itrc))), rec, &
2961 & TRIM(ncname)
2962 END IF
2963 exit_flag=3
2964 ioerror=status
2965 RETURN
2966 END IF
2967 END DO
2968# endif
2969#endif
2970!
2971! If multiple files, close current file.
2972!
2973 IF (ndef.gt.0) THEN
2974 status=nf_close(ncid)
2975 END IF
2976!
2977 10 FORMAT (' READ_STATE - unable to open NetCDF file: ',a)
2978 20 FORMAT (' READ_STATE - error while reading variable: ',a,2x, &
2979 & 'at time record = ',i3,/,16x,'in NetCDF file: ',a)
2980
2981 RETURN
2982 END SUBROUTINE read_state
2983
2984 SUBROUTINE cg_write (ng, innLoop, outLoop)
2985!
2986!=======================================================================
2987! !
2988! This routine writes conjugate gradient vectors into 4DVAR NetCDF !
2989! for restart purposes. !
2990! !
2991!=======================================================================
2992!
2993 USE mod_param
2994 USE mod_parallel
2995 USE mod_fourdvar
2996 Use mod_iounits
2997 USE mod_ncparam
2998 USE mod_netcdf
2999 USE mod_scalars
3000!
3001 implicit none
3002!
3003! Imported variable declarations
3004!
3005 integer, intent(in) :: ng, innLoop, outLoop
3006!
3007! Local variable declarations.
3008!
3009 logical, save :: First = .TRUE.
3010
3011 integer :: i, status
3012 integer :: start(2), total(2)
3013 integer, save :: varid(7)
3014!
3015!-----------------------------------------------------------------------
3016! Write out conjugate gradient vectors.
3017!-----------------------------------------------------------------------
3018!
3019 IF (OutThread) THEN
3020 IF (First) THEN
3021 First=.FALSE.
3022 DO i=1,7
3023 varid(i)=0
3024 END DO
3025 END IF
3026!
3027! Write out outer and inner iteration.
3028!
3029 IF (varid(1).eq.0) THEN
3030 status=nf_inq_varid(ncMODid(ng),'outer',varid(1))
3031 END IF
3032 status=nf_put_var1_int(ncMODid(ng),varid(1),1,outer)
3033 IF (status.ne.nf_noerr) THEN
3034 WRITE (stdout,10) 'outer', TRIM(MODname(ng))
3035 exit_flag=3
3036 ioerror=status
3037 RETURN
3038 END IF
3039
3040 IF (varid(2).eq.0) THEN
3041 status=nf_inq_varid(ncMODid(ng),'inner',varid(2))
3042 END IF
3043 status=nf_put_var1_int(ncMODid(ng),varid(2),1,inner)
3044 IF (status.ne.nf_noerr) THEN
3045 WRITE (stdout,10) 'inner', TRIM(MODname(ng))
3046 exit_flag=3
3047 ioerror=status
3048 RETURN
3049 END IF
3050!
3051! Write out number of converged Ritz eigenvalues.
3052!
3053 IF ((innLoop.eq.0).and.(outloop.eq.1)) THEN
3054 IF (varid(3).eq.0) THEN
3055 status=nf_inq_varid(ncMODid(ng),'nConvRitz',varid(3))
3056 END IF
3057 status=nf_put_var1_int(ncMODid(ng),varid(3),1,nConvRitz)
3058 IF (status.ne.nf_noerr) THEN
3059 WRITE (stdout,10) 'nConvRitz', TRIM(MODname(ng))
3060 exit_flag=3
3061 ioerror=status
3062 RETURN
3063 END IF
3064 END IF
3065!
3066! Write out converged Ritz eigenvalues.
3067!
3068 IF ((innLoop.eq.0).and.(outloop.eq.1)) THEN
3069 IF (varid(4).eq.0) THEN
3070 status=nf_inq_varid(ncMODid(ng),'Ritz',varid(4))
3071 END IF
3072 start(1)=1
3073 total(1)=nConvRitz
3074 status=nf_put_vara_TYPE(ncMODid(ng), varid(4), start, &
3075 & total, Ritz)
3076 IF (status.ne.nf_noerr) THEN
3077 WRITE (stdout,10) 'Ritz', TRIM(MODname(ng))
3078 exit_flag=3
3079 ioerror=status
3080 RETURN
3081 END IF
3082 END IF
3083!
3084! Write out conjugate gradient norms.
3085!
3086 IF (varid(5).eq.0) THEN
3087 status=nf_inq_varid(ncMODid(ng),'cg_alpha',varid(5))
3088 END IF
3089 start(1)=1
3090 total(1)=Ninner+1
3091 start(2)=1
3092 total(2)=Nouter
3093 status=nf_put_vara_TYPE(ncMODid(ng), varid(5), start, &
3094 & total, cg_alpha(0,1))
3095 IF (status.ne.nf_noerr) THEN
3096 WRITE (stdout,10) 'cg_alpha', TRIM(MODname(ng))
3097 exit_flag=3
3098 ioerror=status
3099 RETURN
3100 END IF
3101!
3102 IF (varid(6).eq.0) THEN
3103 status=nf_inq_varid(ncMODid(ng),'cg_beta',varid(6))
3104 END IF
3105 start(1)=1
3106 total(1)=Ninner+1
3107 start(2)=1
3108 total(2)=Nouter
3109 status=nf_put_vara_TYPE(ncMODid(ng), varid(6), start, &
3110 & total, cg_beta(0,1))
3111 IF (status.ne.nf_noerr) THEN
3112 WRITE (stdout,10) 'cg_beta', TRIM(MODname(ng))
3113 exit_flag=3
3114 ioerror=status
3115 RETURN
3116 END IF
3117!
3118 IF (varid(7).eq.0) THEN
3119 status=nf_inq_varid(ncMODid(ng),'cg_tau',varid(7))
3120 END IF
3121 start(1)=1
3122 total(1)=Ninner+1
3123 start(2)=1
3124 total(2)=Nouter
3125 status=nf_put_vara_TYPE(ncMODid(ng), varid(7), start, total, &
3126 & cg_tau(0,1))
3127 IF (status.ne.nf_noerr) THEN
3128 WRITE (stdout,10) 'cg_tau', TRIM(MODname(ng))
3129 exit_flag=3
3130 ioerror=status
3131 RETURN
3132 END IF
3133 END IF
3134!
3135!-----------------------------------------------------------------------
3136! Synchronize observations NetCDF file to disk.
3137!-----------------------------------------------------------------------
3138!
3139 IF (OutThread) THEN
3140 status=nf_sync(ncMODid(ng))
3141 IF (status.ne.nf_noerr) THEN
3142 WRITE (stdout,20)
3143 exit_flag=3
3144 ioerror=status
3145 RETURN
3146 END IF
3147 END IF
3148
3149 10 FORMAT (/,' CG_WRITE - error while writing variable: ',a,/, &
3150 & 12x,'into NetCDF file: ',a)
3151 20 FORMAT (/,' CG_WRITE - unable to synchronize 4DVAR', &
3152 & 1x,'NetCDF file to disk.')
3153
3154 END SUBROUTINE cg_write