ROMS
Loading...
Searching...
No Matches
i4dvar.F
Go to the documentation of this file.
1#include "cppdefs.h"
2 MODULE i4dvar_mod
3
4#ifdef I4DVAR
5!
6!git $Id$
7!================================================== Hernan G. Arango ===
8! Copyright (c) 2002-2025 The ROMS Group Andrew M. Moore !
9! Licensed under a MIT/X style license !
10! See License_ROMS.md !
11!=======================================================================
12! !
13! This module splits the I4D-Var data assimilation algorithm into !
14! its logical components routines. !
15! !
16! background_initialize: !
17! !
18! Initializes the nonlinear model kernel. It is separated from !
19! the 'background' phase to allow ESM coupling. !
20! !
21! background: !
22! !
23! Timesteps the nonlinear model to compute the basic state Xb(t) !
24! used to linearize the tangent linear and adjoint models. It !
25! interpolates the nonlinear model trajectory to the observations !
26! locations. !
27! !
28! increment: !
29! !
30! Minimizes of the 4D-Var cost function over Ninner inner loops !
31! iterations to compute the data assimilation increment, dXa. !
32! !
33! analysis: !
34! !
35! Computes the nonlinear model new initial conditions by adding !
36! the 4D-Var data assimilation increment to the background !
37! initial state, Xa = Xb(t=0) + dXa. Then, writes it into NLM !
38! initial conditions NetCDF INI(ng)%name (record Lini). !
39! !
40! posterior_analysis_initialize: !
41! !
42! Initializes the nonlinear model with the 4D-Var analysis, Xa. !
43! It is separated from 'posterior_analysis' phase to allow ESM !
44! coupling. !
45! !
46! posterior_analysis: !
47! !
48! Timesteps the nonlinear model to compute the posterior analysis !
49! state. It interpolates solution at observation locations !
50! !
51! prior_error: !
52! !
53! Processes the prior background error covariance and its !
54! normalization coefficients. !
55! !
56! !
57! References: !
58! !
59! Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver, !
60! and J. Zavala-Garay, 2011: The Regional Ocean Modeling System !
61! (ROMS) 4-dimensional variational data assimilations systems, !
62! Part I - System overview and formulation, Prog. Oceanogr., 91, !
63! 34-49, doi:10.1016/j.pocean.2011.05.004. !
64! !
65! Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani, !
66! B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson, !
67! 2011: The Regional Ocean Modeling System (ROMS) 4-dimensional !
68! variational data assimilations systems, Part II - Performance !
69! and application to the California Current System, Prog. !
70! Oceanogr., 91, 50-73, doi:10.1016/j.pocean.2011.05.003. !
71! !
72!=======================================================================
73!
74 USE mod_kinds
75 USE mod_param
76 USE mod_parallel
77 USE mod_fourdvar
78 USE mod_iounits
79 USE mod_ncparam
80 USE mod_netcdf
81# if defined PIO_LIB && defined DISTRIBUTE
83# endif
84 USE mod_scalars
85 USE mod_stepping
86!
87# ifdef ADJUST_BOUNDARY
89# endif
90# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
92# endif
94 USE mod_ocean, ONLY : initialize_ocean
95!
96# ifdef BALANCE_OPERATOR
97 USE ad_balance_mod, ONLY : ad_balance
98# endif
101 USE ad_wrt_his_mod, ONLY : ad_wrt_his
102# ifdef STD_MODEL
103 USE background_std_mod, ONLY : background_std
104# endif
105 USE back_cost_mod, ONLY : back_cost
106 USE cgradient_mod, ONLY : cgradient
107# ifdef SPLIT_4DVAR
109# endif
110 USE close_io_mod, ONLY : close_file
111 USE cost_grad_mod, ONLY : cost_grad
112 USE def_hessian_mod, ONLY : def_hessian
113# ifdef SPLIT_4DVAR
114 USE def_ini_mod, ONLY : def_ini
115# endif
116 USE def_mod_mod, ONLY : def_mod
117 USE def_norm_mod, ONLY : def_norm
118# ifdef STD_MODEL
119 USE def_std_mod, ONLY : def_std
120# endif
121 USE get_state_mod, ONLY : get_state
122 USE ini_adjust_mod, ONLY : ini_adjust
124# if defined MASKING && defined SPLIT_4DVAR
125 USE set_masks_mod, ONLY : set_masks
126# endif
127 USE strings_mod, ONLY : founderror
128 USE sum_grad_mod, ONLY : sum_grad
129# ifdef BALANCE_OPERATOR
130 USE tl_balance_mod, ONLY : tl_balance
131# endif
133# ifdef SPLIT_4DVAR
134 USE tl_def_his_mod, ONLY : tl_def_his
135 USE tl_def_ini_mod, ONLY : tl_def_ini
136# endif
138 USE tl_wrt_ini_mod, ONLY : tl_wrt_ini
139# ifdef EVOLVED_LCZ
140 USE wrt_evolved_mod, ONLY : wrt_evolved
141# endif
142# if defined ADJUST_BOUNDARY || defined ADJUST_STFLUX || \
143 defined adjust_wstress
144 USE wrt_ini_mod, ONLY : wrt_frc
145# endif
146 USE wrt_ini_mod, ONLY : wrt_ini
147# ifdef STD_MODEL
148 USE wrt_std_mod, ONLY : wrt_std
149# endif
150# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
152# endif
153!
154 implicit none
155!
156 PUBLIC :: background_initialize
157 PUBLIC :: background
158 PUBLIC :: increment
159 PUBLIC :: analysis
161 PUBLIC :: posterior_analysis
162 PUBLIC :: prior_error
163!
164! Set module internal parameters.
165!
166 logical :: lgetnrm ! process covariance normalization
167 logical :: lgetstd ! process covariance standard deviation
168# ifdef SPLIT_4DVAR
169 logical :: ldone ! 4D-Var cycle finish switch
170# endif
171!
172 integer :: ltlm1 = 1 ! trial x-space TLM IC record in ITL
173 integer :: ltlm2 = 2 ! previous v-space TLM IC record in ITL
174 integer :: ltlm3 = 3 ! trial v-space TLM IC record in ITL
175 integer :: ladj1 = 1 ! initial cost gradient
176 integer :: ladj2 = 2 ! new cost gradient (not normalized)
177 integer :: lini = 1 ! NLM initial conditions record in INI
178 integer :: lbck = 2 ! background record in INI
179 integer :: rec1 = 1
180 integer :: rec2 = 2
181 integer :: rec3 = 3
182 integer :: rec4 = 4
183# ifdef SPLIT_4DVAR
184 integer :: rec5 = 5
185# endif
186!
187 CONTAINS
188!
189 SUBROUTINE background_initialize (my_outer)
190!
191!=======================================================================
192! !
193! This routine initializes ROMS nonlinear model trajectory Xb_n-1(0) !
194! for each (n-1) outer loops. It is separated from the 'background' !
195! 4D-Var phase in ESM coupling applications that use generic methods !
196! for 'initialize', 'run', and 'finalize'. !
197! !
198! On Input: !
199! !
200! my_outer Outer-loop counter (integer) !
201! !
202!=======================================================================
203!
204! Imported variable declarations
205!
206 integer, intent(in) :: my_outer
207!
208! Local variable declarations.
209!
210 integer :: lstr, ng, tile
211# ifdef STD_MODEL
212 integer :: lstd
213# endif
214 integer :: fcount, inprec, tindex
215# ifdef PROFILE
216 integer :: thread
217# endif
218!
219 character (len=*), parameter :: myfile = &
220 & __FILE__//", background_initialize"
221!
222 sourcefile=myfile
223!
224!-----------------------------------------------------------------------
225! Initialize nonlinear model kernel.
226!-----------------------------------------------------------------------
227
228# ifdef PROFILE
229!
230! Start profile clock.
231!
232 DO ng=1,ngrids
233 DO thread=thread_range
234 CALL wclock_on (ng, inlm, 86, __line__, myfile)
235 END DO
236 END DO
237# endif
238!
239! Clear nonlinear mixing arrays.
240!
241 DO ng=1,ngrids
242 DO tile=first_tile(ng),last_tile(ng),+1
243 CALL initialize_mixing (ng, tile, inlm)
244 END DO
245 END DO
246
247# ifdef SPLIT_4DVAR
248!
249!-----------------------------------------------------------------------
250! If split 4D-Var algorithm, set several variables that computed or
251! assigned in other 4D-Var phase executable.
252!-----------------------------------------------------------------------
253!
254! If outer>1, set ERstr=Nrun to compute the open boundary conditions
255! (OBC_time) and surface forcing (SF_time) adjustment times, if needed.
256!
257 IF (my_outer.gt.1) THEN
258 nrun=1+(my_outer-1)*(ninner+1)
259 erstr=nrun
260 END IF
261
262# if defined ADJUST_BOUNDARY || defined ADJUST_STFLUX || \
263 defined adjust_wstress
264!
265! If outer>1, reset reset surface forcing and lateral boundary output
266! time indices to LADJ2 since they are changed in the inner loops.
267!
268 IF (my_outer.gt.1) THEN
269 DO ng=1,ngrids
270 lfout(ng)=ladj2
271# ifdef ADJUST_BOUNDARY
272 lbout(ng)=ladj2
273# endif
274 END DO
275 END IF
276# endif
277!
278! If outer>1, open Nonlinear model initial conditions NetCDF file
279! (INI struc) and inquire about its variable IDs. Reset the value of
280! INI(ng)%Rindex, which is needed in "initial" to the 4D-Var analysis
281! record.
282!
283 IF (my_outer.gt.1) THEN
284 DO ng=1,ngrids
285 ldefini(ng)=.false.
286 CALL def_ini (ng)
287 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
288 ini(ng)%Rindex=lini
289 END DO
290 END IF
291!
292! If outer>1, open 4D-Var NetCDF file (DAV struc) and inquire about its
293! variables.
294!
295 IF (my_outer.gt.1) THEN
296 DO ng=1,ngrids
297 ldefmod(ng)=.false.
298 CALL def_mod (ng)
299 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
300 havenlmod(ng)=.false. ! do not read NLmodVal
301 END DO
302 END IF
303!
304! If outer>1, read in the global observation screening and quality
305! control scale ObsScaleGlobal ("obs_scale") that it is written in the
306! first outer loop. Such data is in memory in the unsplit algorithm.
307! Recall that in I4D-Var, observation variables are read and load into
308! the arrays elements 1:Nobs(ng) for each survey time. That is, only
309! the values needed are read.
310!
311! HGA: What to do in 4D-Var with nested grids?
312!
313 IF (my_outer.gt.1) THEN
314 ng=1
315 SELECT CASE (dav(ng)%IOtype)
316 CASE (io_nf90)
317 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
318 & vname(1,idobss), obsscaleglobal, &
319 & ncid = dav(ng)%ncid)
320
321# if defined PIO_LIB && defined DISTRIBUTE
322 CASE (io_pio)
323 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
324 & vname(1,idobss), obsscaleglobal, &
325 & piofile = dav(ng)%pioFile)
326# endif
327 END SELECT
328 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
329 END IF
330!
331! If outer>1, read in initial (outer=1) nonlinear cost function values
332! from DAV file and load it into "CostNorm", its values are used for
333! reporting normalized values.
334!
335 IF (my_outer.gt.1) THEN
336 DO ng=1,ngrids
337 SELECT CASE (dav(ng)%IOtype)
338 CASE (io_nf90)
339 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
340 & 'NLcost_function', &
341 & fourdvar(ng)%CostNorm(0:), &
342 & ncid = dav(ng)%ncid, &
343 & start = (/1,1/), &
344 & total = (/nobsvar(ng)+1,1/))
345
346# if defined PIO_LIB && defined DISTRIBUTE
347 CASE (io_pio)
348 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
349 & 'NLcost_function', &
350 & fourdvar(ng)%CostNorm(0:), &
351 & piofile = dav(ng)%pioFile, &
352 & start = (/1,1/), &
353 & total = (/nobsvar(ng)+1,1/))
354# endif
355 END SELECT
356 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
357 END DO
358 END IF
359
360# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
361 defined adjust_boundary
362!
363! If outer>1, read in surface forcing and or lateral boundary
364! increments (X-space) from ITL file Rec5.
365!
366 IF (my_outer.gt.1) THEN
367 DO ng=1,ngrids
368 ldefitl(ng)=.false.
369 CALL tl_def_ini (ng)
370 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
371!
372 inprec=rec5
373 CALL get_state (ng, itlm, 1, itl(ng), inprec, ltlm1)
374 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
375 END DO
376 END IF
377# endif
378
379# ifdef FORWARD_FLUXES
380!
381! If outer>1, set structure to process nonlinear surface fluxes from
382! the initial (outer=1) NLM quicksave trajectory. Notice that 2 is
383! substracted since the counter starts at 'outer0'.
384!
385 IF (my_outer.gt.1) THEN
386 DO ng=1,ngrids
387 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-2
388 lstr=len_trim(qck(ng)%name)
389 qck(ng)%base=qck(ng)%name(1:lstr-3)
390 END DO
391!
392 CALL edit_multifile ('QCK2BLK')
393 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
394 END IF
395# endif
396# endif
397!
398!-----------------------------------------------------------------------
399! Initialize nonlinear model.
400!-----------------------------------------------------------------------
401!
402! Activate switch to write out initial misfit between model and
403! observations.
404!
405 DO ng=1,ngrids
406 IF (my_outer.eq.1) THEN
407 wrtmisfit(ng)=.true.
408 ELSE
409 wrtmisfit(ng)=.false.
410 END IF
411 END DO
412!
413! Set nonlinear output history file name. Create a basic state file
414! for each outher loop.
415!
416 DO ng=1,ngrids
417 idefhis(ng)=-1
418 ldefhis(ng)=.true.
419 lwrthis(ng)=.true.
420 lreadfwd(ng)=.false.
421 WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer-1
422 lstr=len_trim(his(ng)%name)
423 his(ng)%base=his(ng)%name(1:lstr-3)
424 END DO
425!
426! Set the nonlinear model to output the quicksave history file as a
427! function of the outer loop. It may be used as the basic state
428! trajectory for the surface fluxes (wind stress, shortwave, heat
429! flux, and E-P) because they can be saved frequently to resolve
430! the daily cycle while avoiding large files.
431!
432 DO ng=1,ngrids
433 ldefqck(ng)=.true.
434 lwrtqck(ng)=.true.
435# ifdef FORWARD_FLUXES
436 lreadblk(ng)=.false.
437# endif
438 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-1
439 lstr=len_trim(qck(ng)%name)
440 qck(ng)%base=qck(ng)%name(1:lstr-3)
441 END DO
442!
443! Initialize nonlinear model. If outer=1, the model is initialized
444! with the background or reference state. Otherwise, the model is
445! initialized with the estimated initial conditions from previous
446! iteration, X(0) = X(0) + deltaX(0).
447!
448 DO ng=1,ngrids
449 wrtnlmod(ng)=.true.
450 wrttlmod(ng)=.false.
451 rst(ng)%Rindex=0
452 fcount=rst(ng)%load
453 rst(ng)%Nrec(fcount)=0
454 END DO
455!
456 CALL initial
457 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
458!
459! If first pass, save nonlinear initial conditions (currently in time
460! index 1, background) into next record (Lbck) of INI(ng)%name NetCDF
461! file. The record "Lbck" becomes the background state record and the
462! record "Lini" becomes current nonlinear initial conditions. Both
463! records are used in the algorithm below.
464!
465 IF (my_outer.eq.1) THEN
466 tindex=1
467 DO ng=1,ngrids
468 ini(ng)%Rindex=1
469 fcount=ini(ng)%load
470 ini(ng)%Nrec(fcount)=1
471# ifdef DISTRIBUTE
472 CALL wrt_ini (ng, myrank, tindex)
473# else
474 CALL wrt_ini (ng, -1, tindex)
475# endif
476 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
477 END DO
478 END IF
479
480# ifdef STD_MODEL
481!
482! Estimate initial conditions standard deviation fields for the
483! background error covariance using the method of Mogensen et al.
484! (2012), which assumes that the background errors are proportional
485! to the vertical derivatives of the state. The error have simmilar
486! shape as the prior profile but iwth a displacement over the water
487! column.
488!
489 IF (my_outer.eq.1) THEN
490 lstd=1 ! time index in error arrays
491 DO ng=1,ngrids
492 DO tile=first_tile(ng),last_tile(ng),+1
493 CALL background_std (ng, tile, tindex, lstd)
494 END DO
495 END DO
496!
497! Write out estimated background standard deviation fields.
498!
499 DO ng=1,ngrids
500 IF (lwrtstd(ng)) THEN
501 std(5,ng)%Rindex=0
502 fcount=std(5,ng)%load
503 std(5,ng)%Nrec(fcount)=1
504# ifdef DISTRIBUTE
505 CALL wrt_std (ng, myrank, lstd)
506# else
507 CALL wrt_std (ng, -1, lstd)
508# endif
509 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
510 END IF
511 END DO
512 END IF
513# endif
514
515# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
516!
517! Compute the reference zeta and biconjugate gradient arrays
518! required for the balance of free surface.
519!
520 IF (balance(isfsur)) THEN
521 DO ng=1,ngrids
522 DO tile=first_tile(ng),last_tile(ng),+1
523 CALL balance_ref (ng, tile, lini)
524 CALL biconj (ng, tile, inlm, lini)
525 END DO
526 wrtzetaref(ng)=.true.
527 END DO
528 END IF
529# endif
530!
531! If first pass, define output 4DVAR NetCDF file containing all
532! processed data at observation locations.
533!
534 IF (my_outer.eq.1) THEN
535 DO ng=1,ngrids
536 ldefmod(ng)=.true.
537 CALL def_mod (ng)
538 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
539 END DO
540 END IF
541!
542 10 FORMAT (a,'_outer',i0,'.nc')
543!
544 RETURN
545 END SUBROUTINE background_initialize
546!
547 SUBROUTINE background (my_outer, RunInterval)
548!
549!=======================================================================
550! !
551! This routine computes the backgound state trajectory, Xb_n-1(t), !
552! used to linearize the tangent linear and adjoint models in the !
553! inner loops. It interpolates the background at the observations !
554! locations, and computes the accept/reject quality control flag, !
555! ObsScale. !
556! !
557! On Input: !
558! !
559! my_outer Outer-loop counter (integer) !
560! RunInterval NLM kernel time stepping window (seconds) !
561! !
562!=======================================================================
563!
564! Imported variable declarations
565!
566 integer, intent(in) :: my_outer
567!
568 real(dp), intent(in) :: runinterval
569!
570! Local variable declarations.
571!
572 logical :: donestepping
573!
574 integer :: i, lstr, ng
575# ifdef PROFILE
576 integer :: thread
577# endif
578# if defined MODEL_COUPLING && !defined MCT_LIB
579 integer :: nstrstep, nendstep, extra
580!
581 real(dp) :: endtime, nexttime
582# endif
583!
584 character (len=*), parameter :: myfile = &
585 & __FILE__//", background"
586!
587 sourcefile=myfile
588
589# if !(defined MODEL_COUPLING && defined ESMF_LIB)
590!
591!-----------------------------------------------------------------------
592! Set nonlinear model initial conditions.
593!-----------------------------------------------------------------------
594!
595 CALL background_initialize (my_outer)
596 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
597# endif
598!
599!-----------------------------------------------------------------------
600! Run nonlinear model. Save nonlinear tracjectory needed by the
601! adjoint and tangent linear models. Interpolate nonlinear model
602! to observation locations (compute and save H x). It processes
603! and writes the observations accept/reject flag (ObsScale) once
604! to allow background quality control, if any.
605# if defined MODEL_COUPLING && !defined MCT_LIB
606! Since the ROMS kernel has a delayed output and line diagnostics by
607! one timestep, subtact an extra value to the report of starting and
608! ending timestep for clarity. Usually, the model coupling interval
609! is of the same size as ROMS timestep.
610# endif
611!-----------------------------------------------------------------------
612
613# ifdef PROFILE
614!
615! Start profile clock.
616!
617 DO ng=1,ngrids
618 DO thread=thread_range
619 CALL wclock_on (ng, inlm, 86, __line__, myfile)
620 END DO
621 END DO
622# endif
623!
624! Initialize various parameters and switches.
625!
626 myruninterval=runinterval
627!
628 DO ng=1,ngrids
629# ifdef AVERAGES
630 idefavg(ng)=-1
631 ldefavg(ng)=.true.
632 lwrtavg(ng)=.true.
633 WRITE (avg(ng)%name,10) trim(avg(ng)%head), my_outer
634 lstr=len_trim(avg(ng)%name)
635 avg(ng)%base=avg(ng)%name(1:lstr-3)
636# endif
637# ifdef DIAGNOSTICS
638 idefdia(ng)=-1
639 ldefdia(ng)=.true.
640 lwrtdia(ng)=.true.
641 WRITE (dia(ng)%name,10) trim(dia(ng)%head), my_outer
642 lstr=len_trim(dia(ng)%name)
643 dia(ng)%base=dia(ng)%name(1:lstr-3)
644# endif
645 wrtobsscale(ng)=.true.
646# if defined MODEL_COUPLING && !defined MCT_LIB
647!
648 nexttime=time(ng)+runinterval
649 endtime=initime(ng)+(ntimes(ng)-1)*dt(ng)
650 IF ((nexttime.eq.endtime).and.(ng.eq.1)) THEN
651 extra=0 ! last time interval
652 ELSE
653 extra=1
654 END IF
655 step_counter(ng)=0
656 nstrstep=iic(ng)
657 nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
658 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0, &
659 & nstrstep, nendstep
660# else
661 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0, &
662 & ntstart(ng), ntend(ng)
663# endif
664 END DO
665!
666# ifdef SOLVE3D
667 CALL main3d (runinterval)
668# else
669 CALL main2d (runinterval)
670# endif
671 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
672!
673!-----------------------------------------------------------------------
674! If completed stepping, write out into NetCDF files.
675# if defined MODEL_COUPLING && !defined MCT_LIB
676! In coupled applications, RunInterval is much less than ntimes*dt,
677! so we need to wait until the last coupling interval is finished.
678! Otherwise, the control switches will be turned off prematurely.
679# endif
680!-----------------------------------------------------------------------
681!
682# if defined MODEL_COUPLING && !defined MCT_LIB
683 IF (nendstep.eq.ntend(1)) THEN
684 donestepping=.true.
685 ELSE
686 donestepping=.false.
687 END IF
688# else
689 donestepping=.true.
690# endif
691!
692 IF (donestepping) THEN
693
694# if defined ADJUST_BOUNDARY || defined ADJUST_STFLUX || \
695 defined adjust_wstress
696!
697! Write out initial and background surface forcing into initial
698! INI(ng)%name NetCDF file for latter use.
699!
700 DO ng=1,ngrids
701# ifdef DISTRIBUTE
702 CALL wrt_frc (ng, myrank, lfout(ng), lini)
703# else
704 CALL wrt_frc (ng, -1, lfout(ng), lini)
705# endif
706 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
707!
708 IF (my_outer.eq.1) THEN
709# ifdef DISTRIBUTE
710 CALL wrt_frc (ng, myrank, lfout(ng), lbck)
711# else
712 CALL wrt_frc (ng, -1, lfout(ng), lbck)
713# endif
714 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
715 END IF
716 END DO
717# endif
718!
719! Write out nonlinear model misfit cost function into DAV(ng)%name
720! NetCDF file.
721!
722 sourcefile=myfile
723 DO ng=1,ngrids
724 SELECT CASE (dav(ng)%IOtype)
725 CASE (io_nf90)
726 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
727 & 'NLcost_function', &
728 & fourdvar(ng)%NLobsCost(0:), &
729 & (/1,my_outer/), &
730 & (/nobsvar(ng)+1,1/), &
731 & ncid = dav(ng)%ncid)
732
733# if defined PIO_LIB && defined DISTRIBUTE
734 CASE (io_pio)
735 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
736 & 'NLcost_function', &
737 & fourdvar(ng)%NLobsCost(0:), &
738 & (/1,my_outer/), &
739 & (/nobsvar(ng)+1,1/), &
740 & piofile = dav(ng)%pioFile)
741# endif
742 END SELECT
743 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
744 END DO
745 END IF
746
747# ifdef PROFILE
748!
749! Stop profile clock
750!
751 DO ng=1,ngrids
752 DO thread=thread_range
753 CALL wclock_off (ng, inlm, 86, __line__, myfile)
754 END DO
755 END DO
756# endif
757!
758 10 FORMAT (a,'_outer',i0,'.nc')
759 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
760 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
761 ', TimeSteps: ',i0,' - ',i0,')',/)
762!
763 RETURN
764 END SUBROUTINE background
765!
766 SUBROUTINE increment (my_outer, RunInterval)
767!
768!=======================================================================
769! !
770! This routine computes the 4D-Var data assimilation state increment, !
771! dXa, by iterating the inner loops and minimizing the cost function. !
772! !
773! On Input: !
774! !
775! my_outer Outer-loop counter (integer) !
776! RunInterval TLM/ADM kernels time stepping window (seconds) !
777! !
778!=======================================================================
779!
780! Imported variable declarations
781!
782 logical :: lweak = .false.
783!
784 integer, intent(in) :: my_outer
785!
786 real(dp), intent(in) :: runinterval
787!
788! Local variable declarations.
789!
790 integer :: i, ifile, lstr, my_inner, ng, tile
791 integer :: fcount, inprec, lcon, lsav
792# ifdef PROFILE
793 integer :: thread
794# endif
795!
796 real(r8) :: rate
797# ifdef SPLIT_4DVAR
798 real(dp) :: stime
799# endif
800!
801 character (len=10) :: suffix
802
803 character (len=*), parameter :: myfile = &
804 & __FILE__//", increment"
805!
806 sourcefile=myfile
807!
808!=======================================================================
809! Compute 4D-Var increment.
810!=======================================================================
811
812# ifdef PROFILE
813!
814! Start profile clock.
815!
816 DO ng=1,ngrids
817 DO thread=thread_range
818 CALL wclock_on (ng, itlm, 87, __line__, myfile)
819 END DO
820 END DO
821# endif
822
823# ifdef SPLIT_4DVAR
824!
825!-----------------------------------------------------------------------
826! If split 4D-Var algorithm, set several variables that are computed
827! or assigned in other 4D-Var phase executable.
828!-----------------------------------------------------------------------
829!
830! Reset Nrun counter to a value greater than one.
831!
832 IF (my_outer.gt.1) THEN
833 nrun=1+(my_outer-1)*(ninner+1)
834 END IF
835!
836! Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
837! Activate "haveNLmod" to read its values in calls to "obs_read". Its
838! values were written in the "background" phase.
839!
840 DO ng=1,ngrids
841 ldefmod(ng)=.false.
842 CALL def_mod (ng)
843 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
844 havenlmod(ng)=.true.
845 END DO
846!
847! In the split 4D-Var algorithm, we need to read the global observation
848! screening and quality control scale ObsScaleGlobal ("obs_scale") that
849! it is written in the "background" phase. Such data is in memory in
850! the unsplit algorithm. Recall that in I4D-Var, observation variables
851! are read and load into the arrays elements 1:Nobs(ng) for each survey
852! time. That is, only the values needed are read.
853!
854! HGA: What to do in 4D-Var with nested grids?
855!
856 ng=1
857 SELECT CASE (dav(ng)%IOtype)
858 CASE (io_nf90)
859 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
860 & vname(1,idobss), obsscaleglobal, &
861 & ncid = dav(ng)%ncid)
862
863# if defined PIO_LIB && defined DISTRIBUTE
864 CASE (io_pio)
865 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
866 & vname(1,idobss), obsscaleglobal, &
867 & piofile = dav(ng)%pioFile)
868# endif
869 END SELECT
870 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
871!
872! Read in first value of NLM cost function computed at outer=1, and
873! load it into "CostNorm" that is used in the reporting of normalized
874! values.
875!
876 sourcefile=myfile
877 DO ng=1,ngrids
878 SELECT CASE (dav(ng)%IOtype)
879 CASE (io_nf90)
880 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
881 & 'NLcost_function', &
882 & fourdvar(ng)%CostNorm(0:), &
883 & ncid = dav(ng)%ncid, &
884 & start = (/1,1/), &
885 & total = (/nobsvar(ng)+1,1/))
886
887# if defined PIO_LIB && defined DISTRIBUTE
888 CASE (io_pio)
889 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
890 & 'NLcost_function', &
891 & fourdvar(ng)%CostNorm(0:), &
892 & piofile = dav(ng)%pioFile, &
893 & start = (/1,1/), &
894 & total = (/nobsvar(ng)+1,1/))
895# endif
896 END SELECT
897 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
898 END DO
899!
900! Read in current outer loop NLM cost function.
901!
902 DO ng=1,ngrids
903 SELECT CASE (dav(ng)%IOtype)
904 CASE (io_nf90)
905 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
906 & 'NLcost_function', &
907 & fourdvar(ng)%NLobsCost(0:), &
908 & ncid = dav(ng)%ncid, &
909 & start = (/1,my_outer/), &
910 & total = (/nobsvar(ng)+1,1/))
911
912# if defined PIO_LIB && defined DISTRIBUTE
913 CASE (io_pio)
914 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
915 & 'NLcost_function', &
916 & fourdvar(ng)%NLobsCost(0:), &
917 & piofile = dav(ng)%pioFile, &
918 & start = (/1,my_outer/), &
919 & total = (/nobsvar(ng)+1,1/))
920# endif
921 END SELECT
922 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
923 END DO
924!
925! If outer>1, read in previous background and TLM cost function to
926! compute FOURDVAR(ng)%CostFun(0). Its value is in memory in the
927! unsplit driver.
928!
929 IF (my_outer.gt.1) THEN
930 DO ng=1,ngrids
931 SELECT CASE (dav(ng)%IOtype)
932 CASE (io_nf90)
933 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
934 & 'TLcost_function', &
935 & fourdvar(ng)%CostFunOld(0:), &
936 & ncid = dav(ng)%ncid, &
937 & start = (/nrun-1/), &
938 & total = (/1/))
939 IF (founderror(exit_flag, noerror, __line__, &
940 & myfile)) RETURN
941!
942 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
943 & 'back_function', &
944 & fourdvar(ng)%BackCost(0:), &
945 & ncid = dav(ng)%ncid, &
946 & start = (/nrun-1/), &
947 & total = (/1/))
948 IF (founderror(exit_flag, noerror, __line__, &
949 & myfile)) RETURN
950
951# if defined PIO_LIB && defined DISTRIBUTE
952 CASE (io_pio)
953 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
954 & 'TLcost_function', &
955 & fourdvar(ng)%CostFunOld(0:), &
956 & piofile = dav(ng)%pioFile, &
957 & start = (/nrun-1/), &
958 & total = (/1/))
959 IF (founderror(exit_flag, noerror, __line__, &
960 & myfile)) RETURN
961!
962 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
963 & 'back_function', &
964 & fourdvar(ng)%BackCost(0:), &
965 & piofile = dav(ng)%pioFile, &
966 & start = (/nrun-1/), &
967 & total = (/1/))
968 IF (founderror(exit_flag, noerror, __line__, &
969 & myfile)) RETURN
970# endif
971 END SELECT
972!
973 fourdvar(ng)%CostFun(0)=fourdvar(ng)%CostFunOld(0)+ &
974 & fourdvar(ng)%BackCost(0)
975 END DO
976 END IF
977!
978! If outer>1, read several conjugate gradient variables from 4D-VAR
979! NetCDF (DAV struc) file. In the unsplit case, such values are
980! available in memory.
981!
982 IF (my_outer.gt.1) THEN
983 DO ng=1,ngrids
984 CALL cg_read_cgradient (ng, itlm, my_outer)
985 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
986 END DO
987 END IF
988!
989! If outer>1, open tangent linear initial conditions NetCDF file
990! (ITL struc) and inquire about its variables IDs.
991!
992 IF (my_outer.gt.1) THEN
993 DO ng=1,ngrids
994 ldefitl(ng)=.false.
995 CALL tl_def_ini (ng)
996 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
997 END DO
998 END IF
999!
1000! If outer>1, open tangent linear history NetCDF file (TLM struc) and
1001! inquire about its variables IDs.
1002!
1003 IF (my_outer.gt.1) THEN
1004 DO ng=1,ngrids
1005 ldeftlm(ng)=.false.
1006 CALL tl_def_his (ng, ldeftlm(ng))
1007 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1008 END DO
1009 END IF
1010!
1011! Set nonlinear output history file to be used as the basis state
1012! trajectory. The 4D-Var increment phase is computed by a different
1013! executable and needs to know some of the HIS structure information.
1014!
1015 DO ng=1,ngrids
1016 WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer-1
1017 lstr=len_trim(his(ng)%name)
1018 his(ng)%base=his(ng)%name(1:lstr-3)
1019 IF (his(ng)%Nfiles.gt.1) THEN
1020 DO ifile=1,his(ng)%Nfiles
1021 WRITE (suffix,"('_',i4.4,'.nc')") ifile
1022 his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
1023 END DO
1024 his(ng)%name=trim(his(ng)%files(1))
1025 ELSE
1026 his(ng)%files(1)=trim(his(ng)%name)
1027 END IF
1028 END DO
1029!
1030! Set the nonlinear model quicksave-history file as the basic state for
1031! the surface fluxes computed in "bulk_flux", which may be available at
1032! more frequent intervals while avoiding large files. Since the 4D-Var
1033! increment phase is calculated by a different executable and needs to
1034! know some of the QCK structure information.
1035!
1036 DO ng=1,ngrids
1037 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-1
1038 lstr=len_trim(qck(ng)%name)
1039 qck(ng)%base=qck(ng)%name(1:lstr-3)
1040 IF (qck(ng)%Nfiles.gt.1) THEN
1041 DO ifile=1,qck(ng)%Nfiles
1042 WRITE (suffix,"('_',i4.4,'.nc')") ifile
1043 qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
1044 END DO
1045 qck(ng)%name=trim(qck(ng)%files(1))
1046 ELSE
1047 qck(ng)%files(1)=trim(qck(ng)%name)
1048 END IF
1049 END DO
1050!
1051! Read in 4D-Var starting time (sec) from nonlinear trajectory.
1052! Initialize "tday" which are needed to write the correct time in
1053! the ITL NetCDF file. It is alse need for boundary and surface
1054! forcing adjustments, if any.
1055!
1056 inprec=1
1057 DO ng=1,ngrids
1058 SELECT CASE (his(ng)%IOtype)
1059 CASE (io_nf90)
1060 CALL netcdf_get_fvar (ng, itlm, his(ng)%name, &
1061 & vname(1,idtime), stime, &
1062 & start = (/inprec/), &
1063 & total = (/1/))
1064
1065# if defined PIO_LIB && defined DISTRIBUTE
1066 CASE (io_pio)
1067 CALL pio_netcdf_get_fvar (ng, itlm, his(ng)%name, &
1068 & vname(1,idtime), stime, &
1069 & start = (/inprec/), &
1070 & total = (/1/))
1071# endif
1072 END SELECT
1073 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1074 initime(ng)=stime
1075
1076# ifdef ADJUST_BOUNDARY
1077!
1078! Set time (sec) for the open boundary adjustment.
1079!
1080 obc_time(1,ng)=stime
1081 DO i=2,nbrec(ng)
1082 obc_time(i,ng)=obc_time(i-1,ng)+nobc(ng)*dt(ng)
1083 END DO
1084# endif
1085
1086# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
1087!
1088! Set time (sec) for the surface forcing adjustment.
1089!
1090 sf_time(1,ng)=stime
1091 DO i=2,nfrec(ng)
1092 sf_time(i,ng)=sf_time(i-1,ng)+nsff(ng)*dt(ng)
1093 END DO
1094# endif
1095 END DO
1096# endif
1097!
1098!-----------------------------------------------------------------------
1099! Initialize.
1100!-----------------------------------------------------------------------
1101!
1102! Set various switches.
1103!
1104 DO ng=1,ngrids
1105 lwrtcost(ng)=.true.
1106# ifdef AVERAGES
1107 ldefavg(ng)=.false.
1108 lwrtavg(ng)=.false.
1109# endif
1110# ifdef DIAGNOSTICS
1111 ldefdia(ng)=.false.
1112 lwrtdia(ng)=.false.
1113# endif
1114 wrtnlmod(ng)=.false.
1115 wrtobsscale(ng)=.false.
1116 wrttlmod(ng)=.true.
1117 END DO
1118!
1119! Set structure for the nonlinear forward trajectory to be processed
1120! by the tangent linear and adjoint models. Also, set switches to
1121! process the FWD structure in routine "check_multifile". Notice that
1122! it is possible to split solution into multiple NetCDF files to reduce
1123! their size.
1124!
1125 CALL edit_multifile ('HIS2FWD')
1126 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1127 DO ng=1,ngrids
1128 lreadfwd(ng)=.true.
1129 END DO
1130
1131# ifdef FORWARD_FLUXES
1132!
1133! Set the BLK structure to contain the nonlinear model surface fluxes
1134! needed by the tangent linear and adjoint models. Also, set switches
1135! to process that structure in routine "check_multifile". Notice that
1136! it is possible to split the solution into multiple NetCDF files to
1137! reduce their size.
1138!
1139! The switch LreadFRC is deactivated because all the atmospheric
1140! forcing, including shortwave radiation, is read from the NLM
1141! surface fluxes or is assigned during ESM coupling. Such fluxes
1142! are available from the QCK structure. There is no need for reading
1143! and processing from the FRC structure input forcing-files.
1144!
1145 CALL edit_multifile ('QCK2BLK')
1146 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1147 DO ng=1,ngrids
1148 lreadblk(ng)=.true.
1149 lreadfrc(ng)=.false.
1150 lreadqck(ng)=.false.
1151 END DO
1152# endif
1153!
1154! Clear the nonlinear state arrays to make sure that the RHS terms
1155! rzeta, rubar, rvbar, ru, and rv are zero. Otherwise, those arrays
1156! will have the last computed values when running the nonlinear model
1157! if not processing the forward trajectory RHS terms. It needs to be
1158! done to get identical solutions with the split schemes.
1159!
1160 DO ng=1,ngrids
1161 DO tile=first_tile(ng),last_tile(ng),+1
1162 CALL initialize_ocean (ng, tile, inlm)
1163 END DO
1164 END DO
1165!
1166! The minimization algorithm requires to save all the gradient
1167! solutions for each inner loop iteration. They are used for
1168! orthogonalization in the conjugate gradient algorithm. Thus,
1169! we need to reset adjoint file record indices.
1170!
1171 DO ng=1,ngrids
1172 adm(ng)%Rindex=0
1173 fcount=adm(ng)%load
1174 adm(ng)%Nrec(fcount)=0
1175 END DO
1176!
1177! An adjoint NetCDF is created for each outer loop.
1178!
1179 DO ng=1,ngrids
1180 idefadj(ng)=-1
1181 ldefadj(ng)=.true.
1182 WRITE (adm(ng)%name,10) trim(adm(ng)%head), my_outer
1183 lstr=len_trim(adm(ng)%name)
1184 adm(ng)%base=adm(ng)%name(1:lstr-3)
1185 END DO
1186!
1187! Define output Hessian NetCDF file containing the eigenvectors
1188! approximation to the Hessian matrix computed from the Lanczos
1189! algorithm. Notice that the file name is a function of the
1190! outer loop. That is, a file is created for each outer loop.
1191!
1192 DO ng=1,ngrids
1193 ldefhss(ng)=.true.
1194 WRITE (hss(ng)%name,10) trim(hss(ng)%head), my_outer
1195 lstr=len_trim(hss(ng)%name)
1196 hss(ng)%base=hss(ng)%name(1:lstr-3)
1197 CALL def_hessian (ng)
1198 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1199 END DO
1200!
1201! Notice that inner loop iteration start from zero. This is needed to
1202! compute the minimization initial increment deltaX(0), its associated
1203! gradient G(0), and descent direction d(0) used in the conjugate
1204! gradient algorithm.
1205!
1206 inner_loop : DO my_inner=0,ninner
1207 inner=my_inner
1208!
1209!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1210! Time-step tangent linear model: compute cost function.
1211!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1212!
1213! If first pass inner=0, initialize tangent linear state (increments,
1214! deltaX) from rest. Otherwise, use trial initial conditions estimated
1215! by the conjugate gradient algorithm in previous inner loop. The TLM
1216! initial conditions are read from ITL(ng)%name, record 1.
1217!
1218 DO ng=1,ngrids
1219 itl(ng)%Rindex=1
1220 CALL tl_initial (ng)
1221 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1222 END DO
1223!
1224! On first pass, initialize records 2, 3 and 4 of the ITL file to zero.
1225!
1226 IF ((my_inner.eq.0).and.(my_outer.eq.1)) THEN
1227 DO ng=1,ngrids
1228# ifdef DISTRIBUTE
1229 CALL tl_wrt_ini (ng, myrank, ltlm1, rec2)
1230# else
1231 CALL tl_wrt_ini (ng, -1, ltlm1, rec2)
1232# endif
1233 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1234# ifdef DISTRIBUTE
1235 CALL tl_wrt_ini (ng, myrank, ltlm1, rec3)
1236# else
1237 CALL tl_wrt_ini (ng, -1, ltlm1, rec3)
1238# endif
1239 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1240# ifdef DISTRIBUTE
1241 CALL tl_wrt_ini (ng, myrank, ltlm1, rec4)
1242# else
1243 CALL tl_wrt_ini (ng, -1, ltlm1, rec4)
1244# endif
1245 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1246 END DO
1247 END IF
1248
1249# ifdef MULTIPLE_TLM
1250!
1251! If multiple TLM history NetCDF files, activate writing and determine
1252! output filename. The multiple file option is used to perturb initial
1253! state and create ensembles. The TLM final trajectory is written for
1254! each inner loop on separated NetCDF files.
1255!
1256 DO ng=1,ngrids
1257 ideftlm(ng)=-1
1258 ldeftlm(ng)=.true.
1259 lwrttlm(ng)=.true.
1260 WRITE (tlm(ng)%name,20) trim(tlm(ng)%head), nrun
1261 lstr=len_trim(tlm(ng)%name)
1262 tlm(ng)%base=tlm(ng)%name(1:lstr-3)
1263 END DO
1264# endif
1265!
1266! Run tangent linear model. Compute misfit observation cost function,
1267! Jo.
1268!
1269 DO ng=1,ngrids
1270 IF (master) THEN
1271 WRITE (stdout,30) 'TL', ng, my_outer, my_inner, &
1272 & ntstart(ng), ntend(ng)
1273 END IF
1274 END DO
1275!
1276# ifdef SOLVE3D
1277 CALL tl_main3d (runinterval)
1278# else
1279 CALL tl_main2d (runinterval)
1280# endif
1281 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1282
1283# ifdef EVOLVED_LCZ
1284!
1285! Write evolved tangent Lanczos vector into hessian netcdf file for use
1286! later.
1287!
1288! NOTE: When using this option, it is important to set LhessianEV and
1289! Lprecond to FALSE in s4dvar.in, otherwise the evolved Lanczos vectors
1290! with be overwritten by the Hessian eigenvectors. The fix to this
1291! is to define a new NetCDF file that contains the evolved Lanczos
1292! vectors.
1293!
1294 IF (my_inner.ne.0) THEN
1295 DO ng=1,ngrids
1296# ifdef DISTRIBUTE
1297 CALL wrt_evolved (ng, myrank, kstp(ng), nrhs(ng))
1298# else
1299 CALL wrt_evolved (ng, -1, kstp(ng), nrhs(ng))
1300# endif
1301 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1302 END DO
1303 END IF
1304# endif
1305
1306# ifdef MULTIPLE_TLM
1307!
1308! If multiple TLM history NetCDF files, close current NetCDF file.
1309!
1310 sourcefile=myfile
1311 DO ng=1,ngrids
1312 CALL close_file (ng, itlm, tlm(ng))
1313 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1314 END DO
1315# endif
1316!
1317!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1318! Time step adjoint model backwards: compute cost function gradient.
1319!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1320!
1321! Initialize the adjoint model always from rest.
1322!
1323 DO ng=1,ngrids
1324 CALL ad_initial (ng)
1325 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1326 END DO
1327!
1328! Time-step adjoint model backwards. The adjoint model is forced with
1329! the adjoint of the observation misfit (Jo) term.
1330!
1331 DO ng=1,ngrids
1332 IF (master) THEN
1333 WRITE (stdout,30) 'AD', ng, my_outer, my_inner, &
1334 & ntstart(ng), ntend(ng)
1335 END IF
1336 END DO
1337!
1338# ifdef SOLVE3D
1339 CALL ad_main3d (runinterval)
1340# else
1341 CALL ad_main2d (runinterval)
1342# endif
1343 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1344!
1345! Clear adjoint arrays.
1346!
1347 DO ng=1,ngrids
1348 DO tile=first_tile(ng),last_tile(ng),+1
1349 CALL initialize_ocean (ng, tile, iadm)
1350# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
1351 CALL initialize_forces (ng, tile, iadm)
1352# endif
1353# ifdef ADJUST_BOUNDARY
1354 CALL initialize_boundary (ng, tile, iadm)
1355# endif
1356 END DO
1357 END DO
1358!
1359!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1360! Conjugate gradient algorithm.
1361!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1362!
1363! Read TLM v-space initial conditions, record 3 in ITL(ng)%name, and
1364! load it into time index LTLM1. This is needed to compute background
1365! cost function. Also read in new (x-space) gradient vector, GRADx(Jo),
1366! from adjoint history file ADM(ng)%name. Read in the sum of all the
1367! previous outer-loop increments which are always in record 4 of
1368! the ITL file.
1369!
1370 DO ng=1,ngrids
1371 IF (my_inner.eq.0) THEN
1372 inprec=rec1
1373 CALL get_state (ng, itlm, 8, itl(ng), inprec, ltlm1)
1374 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1375 ELSE
1376 inprec=rec3
1377 CALL get_state (ng, itlm, 8, itl(ng), inprec, ltlm1)
1378 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1379 END IF
1380
1381 inprec=rec4
1382 CALL get_state (ng, itlm, 8, itl(ng), inprec, ltlm2)
1383 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1384
1385 inprec=adm(ng)%Rindex
1386 CALL get_state (ng, iadm, 4, adm(ng), inprec, ladj2)
1387 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1388
1389# ifdef BALANCE_OPERATOR
1390 inprec=lini
1391 CALL get_state (ng, inlm, 2, ini(ng), inprec, lini)
1392 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1393 nrhs(ng)=lini
1394# endif
1395 END DO
1396!
1397! Convert observation cost function gradient, GRADx(Jo), from model
1398! space (x-space) to minimization space (v-space):
1399!
1400! GRADv(Jo) = B^(T/2) GRADx(Jo), operator: S G L^(T/2) W^(-1/2)
1401!
1402! First, multiply the adjoint solution, GRADx(Jo), by the background-
1403! error standard deviations, S. Second, convolve result with the
1404! adjoint diffusion operator, G L^(T/2) W^(-1/2). Then, backgound
1405! cost function contribution (BackCost) and cost function gradient
1406! (v-space) by adding background and observation contributions:
1407!
1408! GRADv(J) = GRADv(Jb) + GRADv(Jo) = deltaV + GRADv(Jo)
1409!
1410 DO ng=1,ngrids
1411 DO tile=first_tile(ng),last_tile(ng),+1
1412# ifdef BALANCE_OPERATOR
1413 CALL ad_balance (ng, tile, lini, ladj2)
1414# endif
1415 CALL ad_variability (ng, tile, ladj2, lweak)
1416 CALL ad_convolution (ng, tile, ladj2, lweak, 2)
1417 CALL cost_grad (ng, tile, ltlm1, ltlm2, ladj2)
1418 END DO
1419 END DO
1420!
1421! Compute current total cost function.
1422!
1423 DO ng=1,ngrids
1424 IF (nrun.eq.1) THEN
1425 DO i=0,nobsvar(ng)
1426 fourdvar(ng)%CostFunOld(i)=fourdvar(ng)%CostNorm(i)
1427 fourdvar(ng)%CostFun(i)=fourdvar(ng)%CostNorm(i)
1428 END DO
1429 ELSE
1430 DO i=0,nobsvar(ng)
1431 fourdvar(ng)%CostFunOld(i)=fourdvar(ng)%CostFun(i)
1432 END DO
1433 END IF
1434 END DO
1435!
1436! Prepare for background cost function (Jb) calculation:
1437!
1438! Read the convolved gradient from inner=0 (which is permanently
1439! saved in record 1 of the adjoint file) ALWAYS into record 1.
1440!
1441 IF (my_inner.gt.0) THEN
1442 DO ng=1,ngrids
1443 inprec=ladj1
1444 CALL get_state (ng, iadm, 3, adm(ng), inprec, ladj1)
1445 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1446 END DO
1447 END IF
1448!
1449! Compute background cost function (Jb) for inner=0:
1450!
1451! If first pass of inner loop, read in the sum of previous v-space
1452! gradients from record 4 of ITL file using the TLM model variables
1453! as temporary storage. Also add background cost function to Cost0.
1454!
1455 IF (my_inner.eq.0) THEN
1456 DO ng=1,ngrids
1457 inprec=rec4
1458 CALL get_state (ng, itlm, 2, itl(ng), inprec, ltlm2)
1459 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1460!
1461 DO tile=first_tile(ng),last_tile(ng),+1
1462 CALL back_cost (ng, tile, ltlm2)
1463 END DO
1464!
1465 fourdvar(ng)%Cost0(my_outer)=fourdvar(ng)%Cost0(my_outer)+ &
1466 & fourdvar(ng)%BackCost(0)
1467 END DO
1468 END IF
1469!
1470! Compute current total cost function.
1471!
1472 DO ng=1,ngrids
1473 IF (nrun.eq.1) THEN
1474 DO i=0,nobsvar(ng)
1475 fourdvar(ng)%CostNorm(i)=fourdvar(ng)%CostNorm(i)+ &
1476 & fourdvar(ng)%BackCost(i)
1477 fourdvar(ng)%CostFunOld(i)=fourdvar(ng)%CostNorm(i)
1478 fourdvar(ng)%CostFun(i)=fourdvar(ng)%CostNorm(i)
1479 END DO
1480 ELSE
1481 DO i=0,nobsvar(ng)
1482 fourdvar(ng)%CostFunOld(i)=fourdvar(ng)%CostFun(i)
1483 END DO
1484 END IF
1485 END DO
1486!
1487! Determine the descent direction in which the quadractic total cost
1488! function decreases. Then, determine the TLM initial conditions,
1489! deltaV(LTLM1), and its gradient, GRADv{J(Lnew)} at the new
1490! direction. Also, Compute TLM v-space trial initial conditions for
1491! next inner loop, deltaV(LTLM2). The new gradient minimize the
1492! quadratic cost function spanned by current and previous inner loop
1493! iterations. This is achieved by orthogonalizing (Gramm-Schmidt
1494! algorithm) against all previous inner loop gradients.
1495!
1496 DO ng=1,ngrids
1497 DO tile=first_tile(ng),last_tile(ng),+1
1498 CALL cgradient (ng, tile, itlm, my_inner, my_outer)
1499 END DO
1500 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1501 END DO
1502!
1503! Report background (Jb) and observations (Jo) cost function values
1504! normalized by their first minimization value. It also reports the
1505! percentage change on total cost function value with respect to
1506! previous iteration. Compute the optimality of the minimization to
1507! check the statistical hypotheses between the background and
1508! observations errors: the cost function value at the minimum, Jmin,
1509! is ideally equal to half the number of observations assimilated
1510! (Optimality=1=2*Jmin/Nobs), for a linear system.
1511!
1512 IF (master) THEN
1513 DO ng=1,ngrids
1514 IF (nrun.gt.1) THEN
1515 rate=100.0_r8*abs(fourdvar(ng)%CostFun(0)- &
1516 & fourdvar(ng)%CostFunOld(0))/ &
1517 & fourdvar(ng)%CostFunOld(0)
1518 ELSE
1519 rate=0.0_r8
1520 END IF
1521 optimality(ng)=2.0_r8*fourdvar(ng)%CostFun(0)/ &
1522 & (fourdvar(ng)%ObsCount(0)- &
1523 & fourdvar(ng)%ObsReject(0))
1524 WRITE (stdout,40) my_outer, my_inner, &
1525 & fourdvar(ng)%BackCost(0)/ &
1526 & fourdvar(ng)%CostNorm(0), &
1527 & fourdvar(ng)%ObsCost(0)/ &
1528 & fourdvar(ng)%CostNorm(0), &
1529 & rate
1530 IF (my_inner.eq.0) THEN
1531 DO i=0,nobsvar(ng)
1532 IF (fourdvar(ng)%NLobsCost(i).ne.0.0_r8) THEN
1533 IF (i.eq.0) THEN
1534 WRITE (stdout,50) my_outer, my_inner, &
1535 & fourdvar(ng)%NLobsCost(i)/ &
1536 & fourdvar(ng)%CostNorm(i)
1537 ELSE
1538 WRITE (stdout,60) my_outer, my_inner, &
1539 & fourdvar(ng)%NLobsCost(i)/ &
1540 & fourdvar(ng)%CostNorm(i), &
1541 & trim(obsname(i))
1542 END IF
1543 END IF
1544 fourdvar(ng)%NLobsCost(i)=0.0
1545 END DO
1546 END IF
1547 WRITE (stdout,70) my_outer, my_inner, optimality(ng)
1548 END DO
1549 END IF
1550!
1551! Save total v-space cost function gradient, GRADv{J(Lnew)}, into
1552! ADM(ng)%name history NetCDF file. Noticed that the lastest adjoint
1553! solution record is over-written in the NetCDF file for future use.
1554! The switch "LwrtState2d" is activated to write out state arrays
1555! instead ad_*_sol arrays (HGA).
1556!
1557 DO ng=1,ngrids
1558# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
1559 lfout(ng)=ladj2
1560# endif
1561# ifdef ADJUST_BOUNDARY
1562 lbout(ng)=ladj2
1563# endif
1564 kstp(ng)=ladj2
1565# ifdef SOLVE3D
1566 nstp(ng)=ladj2
1567# endif
1568 adm(ng)%Rindex=adm(ng)%Rindex-1
1569 lwrtstate2d(ng)=.true.
1570# ifdef DISTRIBUTE
1571 CALL ad_wrt_his (ng, myrank)
1572# else
1573 CALL ad_wrt_his (ng, -1)
1574# endif
1575 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1576 lwrtstate2d(ng)=.false.
1577 END DO
1578!
1579! Write out trial v-space TLM initial conditions, currently in time
1580! index LTM2, into record 3 of ITL(ng)%name NetCDF file.
1581!
1582 DO ng=1,ngrids
1583# ifdef DISTRIBUTE
1584 CALL tl_wrt_ini (ng, myrank, ltlm2, rec3)
1585# else
1586 CALL tl_wrt_ini (ng, -1, ltlm2, rec3)
1587# endif
1588 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1589 END DO
1590!
1591! Read current outer loop nonlinear model initial conditions and
1592! background state vectors.
1593!
1594 DO ng=1,ngrids
1595 inprec=lini
1596 CALL get_state (ng, inlm, 2, ini(ng), inprec, lini)
1597 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1598
1599 inprec=lbck
1600 CALL get_state (ng, inlm, 9, ini(ng), inprec, lbck)
1601 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1602 END DO
1603!
1604! Convert increment vector, deltaV, from minimization space (v-space)
1605! to model space (x-space):
1606!
1607! deltaX = B^(1/2) deltaV
1608! or
1609! deltaX = W^(1/2) L^(1/2) G S
1610!
1611! First, convolve estimated increment vector (v-space) by with the
1612! tangent linear diffusion operator, W^(1/2) L^(1/2) G. Second,
1613! multiply result by the background-error standard deviation, S.
1614!
1615 lcon=ltlm2
1616 DO ng=1,ngrids
1617 DO tile=first_tile(ng),last_tile(ng),+1
1618 CALL tl_convolution (ng, tile, lcon, lweak, 2)
1619 CALL tl_variability (ng, tile, lcon, lweak)
1620# ifdef BALANCE_OPERATOR
1621 CALL tl_balance (ng, tile, lini, lcon)
1622# endif
1623 END DO
1624 END DO
1625!
1626! Write out trial x-space (convolved) TLM initial conditions, currently
1627! in time index Lcon, into record 1 of ITL(ng)%name NetCDF file.
1628!
1629 DO ng=1,ngrids
1630# ifdef DISTRIBUTE
1631 CALL tl_wrt_ini (ng, myrank, lcon, rec1)
1632# else
1633 CALL tl_wrt_ini (ng, -1, lcon, rec1)
1634# endif
1635 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1636 END DO
1637!
1638!-----------------------------------------------------------------------
1639! Update counters.
1640!-----------------------------------------------------------------------
1641!
1642 DO ng=1,ngrids
1643 lsav=lnew(ng)
1644 lnew(ng)=lold(ng)
1645 lold(ng)=lsav
1646 nrun=nrun+1
1647 END DO
1648!
1649 END DO inner_loop
1650!
1651! Turn of switch to write cost functions in the DAV NetCDF file. They
1652! are written in calls to "tl_wrt_ini" only inside the inner loops.
1653!
1654 DO ng=1,ngrids
1655 lwrtcost(ng)=.false.
1656 END DO
1657!
1658! Close adjoint NetCDF file.
1659!
1660 sourcefile=myfile
1661 DO ng=1,ngrids
1662 CALL close_file (ng, iadm, adm(ng))
1663 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1664 END DO
1665!
1666! Close Hessian NetCDF file.
1667!
1668 DO ng=1,ngrids
1669 CALL close_file (ng, iadm, hss(ng))
1670 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1671 END DO
1672
1673# ifdef PROFILE
1674!
1675! Stop profile clock
1676!
1677 DO ng=1,ngrids
1678 DO thread=thread_range
1679 CALL wclock_off (ng, itlm, 87, __line__, myfile)
1680 END DO
1681 END DO
1682# endif
1683!
1684 10 FORMAT (a,'_outer',i0,'.nc')
1685 20 FORMAT (a,'_member',i3.3,'.nc')
1686 30 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
1687 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
1688 ', TimeSteps: ',i0,' - ',i0,')',/)
1689 40 FORMAT (/,' (',i3.3,',',i3.3,'): TLM Cost Jb, J = ', &
1690 & 1p,e17.10,0p,1x,1p,e17.10,0p,t68,1p,e11.4,' %')
1691 50 FORMAT (/,'>(',i3.3,',',i3.3,'): NLM Cost J = ', &
1692 & 18x,1p,e17.10,0p)
1693 60 FORMAT (' (',i3.3,',',i3.3,'): NLM Cost J = ', &
1694 & 18x,1p,e17.10,0p,t69,a)
1695 70 FORMAT (/,1x,'(',i3.3,',',i3.3,'): Optimality (2*J/Nobs) = ', &
1696 & 1p,e17.10,/)
1697!
1698 RETURN
1699 END SUBROUTINE increment
1700!
1701 SUBROUTINE analysis (my_outer, RunInterval)
1702!
1703!=======================================================================
1704! !
1705! This routine computes 4D-Var data assimilation analysis, Xa. The !
1706! nonlinear model initial conditions are computed by adding the !
1707! 4D-Var increments to the current background: Xa = Xb + dXa. !
1708! !
1709! On Input: !
1710! !
1711! my_outer Outer-loop counter (integer) !
1712! RunInterval NLM kernel time stepping window (seconds) !
1713! !
1714!=======================================================================
1715!
1716! Imported variable declarations
1717!
1718 logical :: lweak = .false.
1719!
1720 integer, intent(in) :: my_outer
1721!
1722 real(dp), intent(in) :: runinterval
1723!
1724! Local variable declarations.
1725!
1726 integer :: ng, tile
1727 integer :: fcount, lcon, nrec
1728 integer :: inprec, outrec
1729# ifdef PROFILE
1730 integer :: thread
1731# endif
1732!
1733 character (len=*), parameter :: myfile = &
1734 & __FILE__//", analysis"
1735!
1736 sourcefile=myfile
1737!
1738!=======================================================================
1739! Compute new nonlinear initial conditions by adding minimization
1740! increment to previous outer loop initial conditions:
1741!
1742! Xi(outer+1) = Xi(outer) + deltaX(Lcon)
1743!
1744!=======================================================================
1745
1746# ifdef PROFILE
1747!
1748! Start profile clock.
1749!
1750 DO ng=1,ngrids
1751 DO thread=thread_range
1752 CALL wclock_on (ng, inlm, 88, __line__, myfile)
1753 END DO
1754 END DO
1755# endif
1756!
1757! Clear nonlinear state variables.
1758!
1759 DO ng=1,ngrids
1760 DO tile=first_tile(ng),last_tile(ng),+1
1761 CALL initialize_ocean (ng, tile, inlm)
1762 END DO
1763 END DO
1764
1765# ifdef SPLIT_4DVAR
1766!
1767!-----------------------------------------------------------------------
1768! If split 4D-Var algorithm, set several variables that computed or
1769! assigned in other 4D-Var phase executable.
1770!-----------------------------------------------------------------------
1771!
1772! Open nonlinear model initial conditions NetCDF file (INI struc) and
1773! inquire about its variable IDs.
1774!
1775 DO ng=1,ngrids
1776 ldefini(ng)=.false.
1777 CALL def_ini (ng)
1778 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1779 END DO
1780!
1781! Open tangent linear model initial conditions NetCDF file (ITL struc)
1782! and inquire about its variable IDs.
1783!
1784 DO ng=1,ngrids
1785 ldefitl(ng)=.false.
1786 CALL tl_def_ini (ng)
1787 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1788 END DO
1789!
1790! Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
1791! It is used in "tl_wrt_ini" to write out observation and background
1792! cost functions.
1793!
1794 DO ng=1,ngrids
1795 ldefmod(ng)=.false.
1796 CALL def_mod (ng)
1797 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1798 nrec=dav(ng)%Rindex
1799 END DO
1800!
1801! Set Nrun>1, set Nrun couter. The counter is only allow to increase
1802! in the "increment" phase.
1803!
1804 nrun=1+(my_outer-1)*(ninner+1)
1805 erstr=nrun
1806
1807# ifdef MASKING
1808!
1809! Set internal mask arrays to process I/O NetCDF files.
1810!
1811 DO ng=1,ngrids
1812 DO tile=first_tile(ng),last_tile(ng),+1
1813 CALL set_masks (ng, tile, inlm)
1814 END DO
1815 END DO
1816# endif
1817# endif
1818!
1819!-----------------------------------------------------------------------
1820! Compute nonlinear model initial from the 4D-Var analysis.
1821!-----------------------------------------------------------------------
1822!
1823! In order to use the correct fields, the model time indices are set
1824! to Lini.
1825!
1826! The appropriate TLM correction for the NLM model resides in record 1
1827! of the ITL file.
1828!
1829 DO ng=1,ngrids
1830 kstp(ng)=lini
1831# ifdef SOLVE3D
1832 nstp(ng)=lini
1833# endif
1834 inprec=lini
1835 CALL get_state (ng, inlm, 1, ini(ng), inprec, lini)
1836 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1837!
1838 inprec=ltlm1
1839 CALL get_state (ng, itlm, 1, itl(ng), inprec, ltlm1)
1840 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1841!
1842 DO tile=first_tile(ng),last_tile(ng),+1
1843 CALL ini_adjust (ng, tile, ltlm1, lini)
1844!! CALL ini_fields (ng, tile, iNLM)
1845 END DO
1846 END DO
1847!
1848! Write out new nonlinear model initial conditions into record Lini
1849! of INI(ng)%name.
1850!
1851 DO ng=1,ngrids
1852 fcount=ini(ng)%load
1853 ini(ng)%Nrec(fcount)=1
1854 outrec=lini
1855# ifdef DISTRIBUTE
1856 CALL wrt_ini (ng, myrank, lini, outrec)
1857# else
1858 CALL wrt_ini (ng, -1, lini, outrec)
1859# endif
1860 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1861 END DO
1862!
1863! Gather the v-space increments from the final inner-loop and
1864! save in record 4 of the ITL file. The current v-space increment
1865! is in record 3 and the sum so far is in record 4.
1866!
1867 DO ng=1,ngrids
1868 inprec=rec3
1869 CALL get_state (ng, itlm, 8, itl(ng), inprec, ltlm1)
1870 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1871!
1872 inprec=rec4
1873 CALL get_state (ng, itlm, 8, itl(ng), inprec, ltlm2)
1874 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1875!
1876 DO tile=first_tile(ng),last_tile(ng),+1
1877 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1878 END DO
1879 END DO
1880!
1881! Write the current sum into record 4 of the ITL file.
1882!
1883 DO ng=1,ngrids
1884# ifdef DISTRIBUTE
1885 CALL tl_wrt_ini (ng, myrank, ltlm2, rec4)
1886# else
1887 CALL tl_wrt_ini (ng, -1, ltlm2, rec4)
1888# endif
1889 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1890 END DO
1891
1892# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
1893 defined adjust_boundary
1894!
1895! Set index containing the surface forcing increments used the run
1896! the nonlinear model in the outer loop and read the forcing
1897! increments. For bulk fluxes, we read Rec1 because the stress
1898! fluxes change by virtue of the changing initial conditions.
1899! When not using bulk fluxes, we read Rec4 because the background
1900! stress and flux is prescribed by input files which are not
1901! overwritten so we need to correct the background using the
1902! sum of the increments from all previous outer-loops.
1903! If using Rec4 we need to convert from v-space to x-space
1904! by applying the convolution.
1905! Note that Lfinp=Lbinp so the the forcing and boundary
1906! adjustments are both processsed correctly.
1907# ifdef BALANCE_OPERATOR
1908! Currently, We do not need the call to tl_balance below, but we
1909! might later if we impose a balance constraint on the wind stress
1910! corrections.
1911# endif
1912!
1913! AMM: CHECK WHAT HAPPENS WITH SECONDARY PRECONDITIONING.
1914!
1915 DO ng=1,ngrids
1916 lfinp(ng)=ltlm1
1917# if defined BULK_FLUXES && !defined FORWARD_FLUXES
1918 inprec=rec1
1919 CALL get_state (ng, itlm, 1, itl(ng), inprec, lfinp(ng))
1920# endif
1921# if defined FORWARD_FLUXES || !defined BULK_FLUXES
1922 inprec=rec4
1923 CALL get_state (ng, itlm, 1, itl(ng), rec4, lfinp(ng))
1924 lcon=lfinp(ng)
1925!
1926 DO tile=first_tile(ng),last_tile(ng),+1
1927 CALL tl_convolution (ng, tile, lcon, lweak, 2)
1928 CALL tl_variability (ng, tile, lcon, lweak)
1929# ifdef BALANCE_OPERATOR
1930!! CALL tl_balance (ng, tile, Lini, Lcon)
1931# endif
1932 END DO
1933# endif
1934 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1935 END DO
1936
1937# ifdef SPLIT_4DVAR
1938!
1939! If split 4D-Var, write out surface forcing and lateral boundary
1940! condition increments (X-space) needed for the nonlinear model into
1941! Rec5.
1942!
1943 DO ng=1,ngrids
1944# ifdef DISTRIBUTE
1945 CALL tl_wrt_ini (ng, myrank, lfinp(ng), rec5)
1946# else
1947 CALL tl_wrt_ini (ng, -1, lfinp(ng), rec5)
1948# endif
1949 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1950 END DO
1951# endif
1952# endif
1953!
1954! Clear tangent linear state variables.
1955!
1956 DO ng=1,ngrids
1957 DO tile=first_tile(ng),last_tile(ng),+1
1958 CALL initialize_ocean (ng, tile, itlm)
1959 END DO
1960 END DO
1961!
1962! Close current forward NetCDF file.
1963!
1964 sourcefile=myfile
1965 DO ng=1,ngrids
1966 CALL close_file (ng, inlm, fwd(ng))
1967 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1968 END DO
1969
1970# ifdef PROFILE
1971!
1972! Stop profile clock
1973!
1974 DO ng=1,ngrids
1975 DO thread=thread_range
1976 CALL wclock_off (ng, inlm, 88, __line__, myfile)
1977 END DO
1978 END DO
1979# endif
1980!
1981 RETURN
1982 END SUBROUTINE analysis
1983!
1985!
1986!=======================================================================
1987! !
1988! This routine initializes the nonlinear kernel with the 4D-Var new !
1989! state estimate Xa = Xb + dXa. It is separated from the 'analysis' !
1990! 4D-Var phase in ESM coupling applications that use generic methods !
1991! for 'initialize', 'run', and 'finalize'. !
1992! !
1993! On Input: !
1994! !
1995! my_outer Outer-loop counter (integer) !
1996! !
1997!=======================================================================
1998!
1999! Local variable declarations.
2000!
2001 integer :: i, ifile, lstr, ng, tile
2002 integer :: fcount, inprec
2003# ifdef PROFILE
2004 integer :: thread
2005# endif
2006!
2007 character (len=*), parameter :: myfile = &
2008 & __FILE__//", posterior_analysis_initialize"
2009!
2010 sourcefile=myfile
2011!
2012!-----------------------------------------------------------------------
2013! Initialize nonlinear model kernel with new 4D-Var state estimate.
2014!-----------------------------------------------------------------------
2015
2016# ifdef PROFILE
2017!
2018! Start profile clock.
2019!
2020 DO ng=1,ngrids
2021 DO thread=thread_range
2022 CALL wclock_on (ng, inlm, 88, __line__, myfile)
2023 END DO
2024 END DO
2025# endif
2026!
2027! Clear nonlinear mixing arrays.
2028!
2029 DO ng=1,ngrids
2030 DO tile=first_tile(ng),last_tile(ng),+1
2031 CALL initialize_mixing (ng, tile, inlm)
2032 END DO
2033 END DO
2034
2035# ifdef SPLIT_4DVAR
2036!
2037!-----------------------------------------------------------------------
2038! If split 4D-Var algorithm, set several variables that computed or
2039! assigned in other 4D-Var phase executable.
2040!-----------------------------------------------------------------------
2041!
2044!
2045! Set Nrun>1, to read in surface forcing and open boundary conditions
2046! increments in "initial", if appropriate.
2047!
2048 nrun=1+nouter*(ninner+1)
2049!
2050! Set ERstr=Nrun, to set the open boundary condition (OBC_time) and
2051! surface forcing (SF_time) adjustment times, if needed.
2052!
2053 erstr=nrun
2054!
2055! Open Nonlinear model initial conditions NetCDF file (INI struc) and
2056! inquire about its variable IDs.
2057!
2058 DO ng=1,ngrids
2059 ldefini(ng)=.false.
2060 CALL def_ini (ng)
2061 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2062 END DO
2063!
2064! Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
2065!
2066 DO ng=1,ngrids
2067 ldefmod(ng)=.false.
2068 CALL def_mod (ng)
2069 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2070 havenlmod(ng)=.false. ! do not read NLmodVal
2071 END DO
2072!
2073! In the split 4D-Var algorithm, we need to read the global observation
2074! screening and quality control scale ObsScaleGlobal ("obs_scale") that
2075! it is written in the "background" phase. Such data is in memory in
2076! the unsplit algorithm. Recall that in I4D-Var, observation variables
2077! are read and load into the arrays elements 1:Nobs(ng) for each survey
2078! time. That is, only the values needed are read.
2079!
2080! HGA: What to do in 4D-Var with nested grids?
2081!
2082 ng=1
2083 SELECT CASE (dav(ng)%IOtype)
2084 CASE (io_nf90)
2085 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
2086 & vname(1,idobss), obsscaleglobal, &
2087 & ncid = dav(ng)%ncid)
2088
2089# if defined PIO_LIB && defined DISTRIBUTE
2090 CASE (io_pio)
2091 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
2092 & vname(1,idobss), obsscaleglobal, &
2093 & piofile = dav(ng)%pioFile)
2094# endif
2095 END SELECT
2096 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2097!
2098! Read initial nonlinear cost function values from DAV file and load it
2099! into CostNorm, its values are used for reporting normalized values.
2100!
2101 DO ng=1,ngrids
2102 SELECT CASE (dav(ng)%IOtype)
2103 CASE (io_nf90)
2104 CALL netcdf_get_fvar (ng, inlm, dav(ng)%name, &
2105 & 'NLcost_function', &
2106 & fourdvar(ng)%CostNorm(0:), &
2107 & ncid = dav(ng)%ncid, &
2108 & start = (/1,1/), &
2109 & total = (/nobsvar(ng)+1,1/))
2110
2111# if defined PIO_LIB && defined DISTRIBUTE
2112 CASE (io_pio)
2113 CALL pio_netcdf_get_fvar (ng, inlm, dav(ng)%name, &
2114 & 'NLcost_function', &
2115 & fourdvar(ng)%CostNorm(0:), &
2116 & piofile = dav(ng)%pioFile, &
2117 & start = (/1,1/), &
2118 & total = (/nobsvar(ng)+1,1/))
2119# endif
2120 END SELECT
2121 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2122 END DO
2123
2124# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
2125 defined adjust_boundary
2126!
2127! If split 4D-Var, read in surface forcing and or lateral boundary
2128! increments (X-space) from ITL file Rec5.
2129!
2130 DO ng=1,ngrids
2131 ldefitl(ng)=.false.
2132 CALL tl_def_ini (ng)
2133 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2134!
2135 inprec=rec5
2136 CALL get_state (ng, itlm, 1, itl(ng), inprec, ltlm1)
2137 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2138 END DO
2139!
2140! If split 4D-Var, reset reset surface forcing and lateral boundary
2141! output time indices to LADJ2 since they are changed in the inner
2142! loops.
2143!
2144 DO ng=1,ngrids
2145 lfout(ng)=ladj2
2146# ifdef ADJUST_BOUNDARY
2147 lbout(ng)=ladj2
2148# endif
2149 END DO
2150# endif
2151
2152# ifdef FORWARD_FLUXES
2153!
2154! If not first NLM run, set BLK structure to the previous quicksave
2155! trajectory. There is a logic in "get_data" that reads from BLK.
2156! (HGA: This probably legacy code that it is no longer needed)
2157!
2158 DO ng=1,ngrids
2159 WRITE (qck(ng)%name,10) trim(qck(ng)%head), nouter-1
2160 lstr=len_trim(qck(ng)%name)
2161 qck(ng)%base=qck(ng)%name(1:lstr-3)
2162 END DO
2163!
2164 CALL edit_multifile ('QCK2BLK')
2165 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2166# endif
2167# endif
2168!
2169!-----------------------------------------------------------------------
2170! Initialize nonlinear model.
2171!-----------------------------------------------------------------------
2172!
2173! Set nonlinear output history file name. Create a basic state file
2174! for each outher loop. Notice that the LreadBLK and LreadFWD switches
2175! are turned off to suppress processing of the structures when
2176! "check_multifile" during nonlinear model execution.
2177!
2178 DO ng=1,ngrids
2179 idefhis(ng)=-1
2180 ldefhis(ng)=.true.
2181 lwrthis(ng)=.true.
2182# ifdef FORWARD_FLUXES
2183 lreadblk(ng)=.false.
2184# endif
2185 lreadfwd(ng)=.false.
2186 his(ng)%Rindex=0
2187 fcount=his(ng)%load
2188 his(ng)%Nrec(fcount)=0
2189 WRITE (his(ng)%name,10) trim(fwd(ng)%head), nouter
2190 lstr=len_trim(his(ng)%name)
2191 his(ng)%base=his(ng)%name(1:lstr-3)
2192 END DO
2193!
2194! Set the nonlinear model output quicksave-history file for the
2195! posterior analysis trajectory. Notice that the LreadBLK switch
2196! is turned off to suppress processing of the structure in
2197! "check_multifile" during nonlinear model execution.
2198!
2199 DO ng=1,ngrids
2200 idefqck(ng)=-1
2201 ldefqck(ng)=.true.
2202 lwrtqck(ng)=.true.
2203# ifdef FORWARD_FLUXES
2204 lreadblk(ng)=.false.
2205# endif
2206 WRITE (qck(ng)%name,10) trim(qck(ng)%head), nouter
2207 lstr=len_trim(qck(ng)%name)
2208 qck(ng)%base=qck(ng)%name(1:lstr-3)
2209 END DO
2210!
2211! Initialize nonlinear model with estimated initial conditions. Reset
2212! the value of INI(ng)%Rindex, which is needed in "initial" to the
2213! 4D-Var analysis record.
2214!
2215 DO ng=1,ngrids
2216 wrtnlmod(ng)=.true.
2217 wrttlmod(ng)=.false.
2218 wrtmisfit(ng)=.false.
2219 rst(ng)%Rindex=0
2220 fcount=rst(ng)%load
2221 rst(ng)%Nrec(fcount)=0
2222 ini(ng)%Rindex=lini
2223 END DO
2224!
2225 CALL initial
2226 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2227!
2228! Clear NLobsCost. Activate switch to write out final misfit between
2229! model and observations.
2230!
2231 DO ng=1,ngrids
2232 DO i=0,nobsvar(ng)
2233 fourdvar(ng)%NLobsCost(i)=0.0_r8
2234 END DO
2235 wrtmisfit(ng)=.true.
2236 END DO
2237
2238# ifdef PROFILE
2239!
2240! Stop profile clock
2241!
2242 DO ng=1,ngrids
2243 DO thread=thread_range
2244 CALL wclock_off (ng, inlm, 88, __line__, myfile)
2245 END DO
2246 END DO
2247# endif
2248!
2249 10 FORMAT (a,'_outer',i0,'.nc')
2250!
2251 RETURN
2252 END SUBROUTINE posterior_analysis_initialize
2253!
2254 SUBROUTINE posterior_analysis (RunInterval)
2255!
2256!=======================================================================
2257! !
2258! This routine initialize the NLM with estimated 4D-Var state and !
2259! interpolates solution at observation locations for posterior !
2260! analysis. !
2261! !
2262! On Input: !
2263! !
2264! RunInterval NLM kernel time stepping window (seconds) !
2265! !
2266!=======================================================================
2267!
2268! Imported variable declarations
2269!
2270 real(dp), intent(in) :: runinterval
2271!
2272! Local variable declarations.
2273!
2274 logical :: donestepping
2275!
2276 integer :: i, lstr, ng
2277# ifdef PROFILE
2278 integer :: thread
2279# endif
2280# if defined MODEL_COUPLING && !defined MCT_LIB
2281 integer :: nstrstep, nendstep, extra
2282!
2283 real(dp) :: endtime, nexttime
2284# endif
2285!
2286 character (len=*), parameter :: myfile = &
2287 & __FILE__//", posterior_analysis"
2288!
2289 sourcefile=myfile
2290
2291# if !(defined MODEL_COUPLING && defined ESMF_LIB)
2292!
2293!-----------------------------------------------------------------------
2294! Set nonlinear model initial conditions.
2295!-----------------------------------------------------------------------
2296!
2298 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2299# endif
2300!
2301!-----------------------------------------------------------------------
2302! Run nonlinear model. Interpolate nonlinear model to observation
2303! locations.
2304# if defined MODEL_COUPLING && !defined MCT_LIB
2305! Since the ROMS kernel has a delayed output and line diagnostics by
2306! one timestep, subtact an extra value to the report of starting and
2307! ending timestep for clarity. Usually, the model coupling interval
2308! is of the same size as ROMS timestep.
2309# endif
2310!-----------------------------------------------------------------------
2311
2312# ifdef PROFILE
2313!
2314! Start profile clock.
2315!
2316 DO ng=1,ngrids
2317 DO thread=thread_range
2318 CALL wclock_on (ng, inlm, 88, __line__, myfile)
2319 END DO
2320 END DO
2321# endif
2322!
2323! Initialize various parameters and switches.
2324!
2325 myruninterval=runinterval
2326!
2327 DO ng=1,ngrids
2328# ifdef AVERAGES
2329 idefavg(ng)=-1
2330 ldefavg(ng)=.true.
2331 lwrtavg(ng)=.true.
2332 WRITE (avg(ng)%name,10) trim(avg(ng)%head), outer
2333 lstr=len_trim(avg(ng)%name)
2334 avg(ng)%base=avg(ng)%name(1:lstr-3)
2335# endif
2336# ifdef DIAGNOSTICS
2337 idefdia(ng)=-1
2338 ldefdia(ng)=.true.
2339 lwrtdia(ng)=.true.
2340 WRITE (dia(ng)%name,10) trim(dia(ng)%head), outer
2341 lstr=len_trim(dia(ng)%name)
2342 dia(ng)%base=dia(ng)%name(1:lstr-3)
2343# endif
2344# if defined MODEL_COUPLING && !defined MCT_LIB
2345!
2346 nexttime=time(ng)+runinterval
2347 endtime=initime(ng)+(ntimes(ng)-1)*dt(ng)
2348 IF ((nexttime.eq.endtime).and.(ng.eq.1)) THEN
2349 extra=0 ! last time interval
2350 ELSE
2351 extra=1
2352 END IF
2353 step_counter(ng)=0
2354 nstrstep=iic(ng)
2355 nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
2356 IF (master) WRITE (stdout,20) 'NL', ng, nstrstep, nendstep
2357# else
2358 IF (master) WRITE (stdout,20) 'NL', ng, ntstart(ng), ntend(ng)
2359# endif
2360 END DO
2361!
2362# ifdef SOLVE3D
2363 CALL main3d (runinterval)
2364# else
2365 CALL main2d (runinterval)
2366# endif
2367 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2368!
2369!-----------------------------------------------------------------------
2370! If completed stepping, write out and report cost function.
2371# if defined MODEL_COUPLING && !defined MCT_LIB
2372! In coupled applications, RunInterval is much less than ntimes*dt,
2373! so we need to wait until the last coupling interval is finished.
2374! Otherwise, the control switches will be turned off prematurely.
2375# endif
2376!-----------------------------------------------------------------------
2377!
2378# if defined MODEL_COUPLING && !defined MCT_LIB
2379 IF (nendstep.eq.ntend(1)) THEN
2380 donestepping=.true.
2381 ELSE
2382 donestepping=.false.
2383 END IF
2384# else
2385 donestepping=.true.
2386# endif
2387!
2388 IF (donestepping) THEN
2389!
2390! Write out nonlinear model final misfit cost function into
2391! DAV(ng)%name NetCDF file. Notice that it is written in the
2392! Nouter+1 record.
2393!
2394 sourcefile=myfile
2395 DO ng=1,ngrids
2396 SELECT CASE (dav(ng)%IOtype)
2397 CASE (io_nf90)
2398 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
2399 & 'NLcost_function', &
2400 & fourdvar(ng)%NLobsCost(0:), &
2401 & (/1,nouter+1/), &
2402 & (/nobsvar(ng)+1,1/), &
2403 & ncid = dav(ng)%ncid)
2404
2405# if defined PIO_LIB && defined DISTRIBUTE
2406 CASE (io_pio)
2407 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
2408 & 'NLcost_function', &
2409 & fourdvar(ng)%NLobsCost(0:), &
2410 & (/1,nouter+1/), &
2411 & (/nobsvar(ng)+1,1/), &
2412 & piofile = dav(ng)%pioFile)
2413# endif
2414 END SELECT
2415 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2416 END DO
2417!
2418! Report the final value of the nonlinear model misfit cost function.
2419!
2420 IF (master) THEN
2421 DO ng=1,ngrids
2422 DO i=0,nobsvar(ng)
2423 IF (fourdvar(ng)%NLobsCost(i).ne.0.0_r8) THEN
2424 IF (i.eq.0) THEN
2425 WRITE (stdout,30) outer, inner, &
2426 & fourdvar(ng)%NLobsCost(i)/ &
2427 & fourdvar(ng)%CostNorm(i)
2428 ELSE
2429 WRITE (stdout,40) outer, inner, &
2430 & fourdvar(ng)%NLobsCost(i)/ &
2431 & fourdvar(ng)%CostNorm(i), &
2432 & trim(obsname(i))
2433 END IF
2434 END IF
2435 END DO
2436 END DO
2437 END IF
2438!
2439! Done. Set history file ID to closed state since we manipulated
2440! its indices with the forward file ID which was closed above.
2441!
2442 DO ng=1,ngrids
2443 his(ng)%ncid=-1
2444 END DO
2445 END IF
2446
2447# ifdef PROFILE
2448!
2449! Stop profile clock
2450!
2451 DO ng=1,ngrids
2452 DO thread=thread_range
2453 CALL wclock_off (ng, inlm, 88, __line__, myfile)
2454 END DO
2455 END DO
2456# endif
2457!
2458 10 FORMAT (a,'_outer',i0,'.nc')
2459 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
2460 & ' (Grid: ',i0,', TimeSteps: ',i0,' - ',i0,')',/)
2461 30 FORMAT (/,'>(',i3.3,',',i3.3,'): NLM Cost J = ', &
2462 & 18x,1p,e17.10,0p)
2463 40 FORMAT (' (',i3.3,',',i3.3,'): NLM Cost J = ', &
2464 & 18x,1p,e17.10,0p,t69,a)
2465!
2466 RETURN
2467 END SUBROUTINE posterior_analysis
2468!
2469 SUBROUTINE prior_error (ng)
2470!
2471!=======================================================================
2472! !
2473! This routine processes background prior error covariance standard !
2474! deviations and normalization coefficients. !
2475! !
2476! On Input: !
2477! !
2478! ng Nested grid number !
2479! !
2480!=======================================================================
2481!
2482! Imported variable declarations
2483!
2484 integer, intent(in) :: ng
2485!
2486! Local variable declarations.
2487!
2488 integer :: tile
2489 integer :: nrmrec, stdrec, tindex
2490!
2491 character (len=*), parameter :: myfile = &
2492 & __FILE__//", prior_error"
2493!
2494 sourcefile=myfile
2495!
2496!-----------------------------------------------------------------------
2497! Set application grid, metrics, and associated variables and
2498! parameters.
2499!-----------------------------------------------------------------------
2500!
2501! The ROMS application grid configuration is done once. It is usually
2502! done in the "initial" kernel routine. However, since we are calling
2503! the "normalization" routine here, we need several grid variables and
2504! parameter. Also, if reading only water points, we need to know the
2505! land/sea mask arrays to unpack.
2506!
2507 IF (setgridconfig(ng)) THEN
2508 CALL set_grid (ng, inlm)
2509 END IF
2510!
2511!-----------------------------------------------------------------------
2512! Read in standard deviation factors for error covariance.
2513!-----------------------------------------------------------------------
2514!
2515 IF (lgetstd) THEN
2516!
2517! Initial conditions standard deviation. They are loaded in Tindex=1
2518! of the e_var(...,Tindex) state variables.
2519!
2520 IF (lreadstd(ng)) THEN
2521 stdrec=1
2522 tindex=1
2523# ifdef STD_MODEL
2524 CALL get_state (ng, 10, 10, std(5,ng), stdrec, tindex)
2525# else
2526 CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
2527# endif
2528 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2529 END IF
2530
2531# ifdef ADJUST_BOUNDARY
2532!
2533! Open boundary conditions standard deviation.
2534!
2535 stdrec=1
2536 tindex=1
2537 CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
2538 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2539# endif
2540
2541# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
2542!
2543! Surface forcing standard deviation.
2544!
2545 stdrec=1
2546 tindex=1
2547 CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
2548 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2549# endif
2550
2551# ifdef STD_MODEL
2552!
2553! If computing the standard deviation (STD) directly from the
2554! background (prior) field as an alternative to climatological
2555! read from the input NetCDF file, define the output STD NetCDF
2556! file.
2557!
2558 IF (ldefstd(ng)) THEN
2559 CALL def_std (ng)
2560 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2561 END IF
2562# endif
2563 END IF
2564!
2565!-----------------------------------------------------------------------
2566! Error covariance normalization coefficients.
2567!-----------------------------------------------------------------------
2568!
2569 IF (lgetnrm) THEN
2570!
2571! Compute or read in background-error covariance normalization factors.
2572! If computing, write out factors to NetCDF. This is an expensive
2573! computation that needs to be computed only once for a particular
2574! application grid.
2575!
2576 IF (any(lwrtnrm(:,ng))) THEN
2577 CALL def_norm (ng, inlm, 1)
2578 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2579
2580# ifdef ADJUST_BOUNDARY
2581 CALL def_norm (ng, inlm, 3)
2582 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2583# endif
2584
2585# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
2586 CALL def_norm (ng, inlm, 4)
2587 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2588# endif
2589!
2590 DO tile=first_tile(ng),last_tile(ng),+1
2591 CALL normalization (ng, tile, 2)
2592 END DO
2593 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2594 ldefnrm(1:4,ng)=.false.
2595 lwrtnrm(1:4,ng)=.false.
2596
2597 ELSE
2598
2599 nrmrec=1
2600 CALL get_state (ng, 14, 14, nrm(1,ng), nrmrec, 1)
2601 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2602
2603# ifdef ADJUST_BOUNDARY
2604 CALL get_state (ng, 16, 16, nrm(3,ng), nrmrec, 1)
2605 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2606# endif
2607
2608# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
2609 CALL get_state (ng, 17, 17, nrm(4,ng), nrmrec, 1)
2610 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2611# endif
2612
2613 END IF
2614 END IF
2615!
2616 RETURN
2617 END SUBROUTINE prior_error
2618#endif
2619 END MODULE i4dvar_mod
subroutine ad_initial(ng)
Definition ad_initial.F:4
subroutine ad_main2d
Definition ad_main2d.F:586
subroutine ad_main3d(runinterval)
Definition ad_main3d.F:4
subroutine edit_multifile(task)
subroutine initial
Definition initial.F:3
subroutine main2d
Definition main2d.F:746
subroutine main3d(runinterval)
Definition main3d.F:4
subroutine, public ad_balance(ng, tile, lbck, linp)
Definition ad_balance.F:59
subroutine, public ad_convolution(ng, tile, linp, lweak, ifac)
subroutine, public ad_variability(ng, tile, linp, lweak)
subroutine, public ad_wrt_his(ng, tile)
Definition ad_wrt_his.F:66
subroutine back_cost(ng, tile, lsum)
Definition back_cost.F:71
subroutine, public cgradient(ng, tile, model, innloop, outloop)
Definition cgradient.F:205
subroutine, public cg_read_cgradient(ng, model, outloop)
Definition cgradient.F:7536
subroutine, public close_file(ng, model, s, ncname, lupdate)
Definition close_io.F:43
subroutine, public cost_grad(ng, tile, linp1, linp2, lout)
Definition cost_grad.F:43
subroutine, public def_hessian(ng)
Definition def_hessian.F:57
subroutine, public def_ini(ng)
Definition def_ini.F:41
subroutine, public def_mod(ng)
Definition def_mod.F:49
subroutine, public def_norm(ng, model, ifile)
Definition def_norm.F:53
subroutine, public get_state(ng, model, msg, s, inirec, tindex)
Definition get_state.F:90
logical lgetnrm
Definition i4dvar.F:166
integer rec5
Definition i4dvar.F:184
subroutine, public prior_error(ng)
Definition i4dvar.F:2470
integer lini
Definition i4dvar.F:177
integer rec1
Definition i4dvar.F:179
integer ladj2
Definition i4dvar.F:176
integer ladj1
Definition i4dvar.F:175
integer rec2
Definition i4dvar.F:180
logical ldone
Definition i4dvar.F:169
subroutine, public increment(my_outer, runinterval)
Definition i4dvar.F:767
subroutine, public background_initialize(my_outer)
Definition i4dvar.F:190
logical lgetstd
Definition i4dvar.F:167
integer lbck
Definition i4dvar.F:178
integer rec4
Definition i4dvar.F:182
subroutine, public background(my_outer, runinterval)
Definition i4dvar.F:548
subroutine, public posterior_analysis(runinterval)
Definition i4dvar.F:2255
integer ltlm2
Definition i4dvar.F:173
integer rec3
Definition i4dvar.F:181
integer ltlm1
Definition i4dvar.F:172
subroutine, public analysis(my_outer, runinterval)
Definition i4dvar.F:1702
subroutine, public posterior_analysis_initialize
Definition i4dvar.F:1985
integer ltlm3
Definition i4dvar.F:174
subroutine, public ini_adjust(ng, tile, linp, lout)
Definition ini_adjust.F:69
subroutine, public initialize_boundary(ng, tile, model)
subroutine, public initialize_forces(ng, tile, model)
type(t_fourdvar), dimension(:), allocatable fourdvar
real(r8), dimension(:), allocatable optimality
logical, dimension(:), allocatable wrtmisfit
integer, dimension(:), allocatable nobsvar
logical, dimension(:), allocatable havenlmod
logical, dimension(:), allocatable wrttlmod
logical, dimension(:), allocatable wrtnlmod
logical, dimension(:), allocatable wrtobsscale
character(len=40), dimension(:), allocatable obsname
logical, dimension(:), allocatable wrtzetaref
type(t_io), dimension(:,:), allocatable std
type(t_io), dimension(:), allocatable his
type(t_io), dimension(:,:), allocatable nrm
type(t_io), dimension(:), allocatable adm
type(t_io), dimension(:), allocatable hss
type(t_io), dimension(:), allocatable tlm
type(t_io), dimension(:), allocatable itl
type(t_io), dimension(:), allocatable dav
type(t_io), dimension(:), allocatable qck
type(t_io), dimension(:), allocatable fwd
type(t_io), dimension(:), allocatable rst
type(t_io), dimension(:), allocatable ini
type(t_io), dimension(:), allocatable avg
integer stdout
character(len=256) sourcefile
type(t_io), dimension(:), allocatable dia
integer, parameter r8
Definition mod_kinds.F:28
integer, parameter dp
Definition mod_kinds.F:25
subroutine, public initialize_mixing(ng, tile, model)
integer idobss
integer, parameter io_nf90
Definition mod_ncparam.F:95
integer, dimension(:), allocatable idefadj
integer, dimension(:), allocatable idefqck
integer, dimension(:), allocatable idefavg
integer, parameter io_pio
Definition mod_ncparam.F:96
integer, dimension(:), allocatable idefdia
integer, dimension(:), allocatable ideftlm
integer isfsur
character(len=maxlen), dimension(6, 0:nv) vname
integer idtime
integer, dimension(:), allocatable idefhis
subroutine, public initialize_ocean(ng, tile, model)
Definition mod_ocean.F:1526
integer, dimension(:), allocatable first_tile
logical master
integer, dimension(:), allocatable last_tile
integer, parameter inlm
Definition mod_param.F:662
integer, parameter iadm
Definition mod_param.F:665
integer ngrids
Definition mod_param.F:113
integer, parameter itlm
Definition mod_param.F:663
integer ninner
logical, dimension(:,:), allocatable lwrtnrm
logical, dimension(:), allocatable lwrtqck
logical, dimension(:), allocatable lwrtdia
logical, dimension(:), allocatable lreadqck
integer nouter
integer, dimension(:), allocatable ntimes
logical, dimension(:), allocatable ldefitl
integer, dimension(:), allocatable iic
logical, dimension(:), allocatable lreadstd
real(dp), dimension(:), allocatable dt
logical, dimension(:), allocatable ldefhss
integer, dimension(:), allocatable nobc
logical, dimension(:), allocatable setgridconfig
logical, dimension(:), allocatable balance
logical, dimension(:), allocatable ldefdia
logical, dimension(:), allocatable lreadfrc
integer, dimension(:), allocatable nfrec
logical, dimension(:), allocatable ldefini
logical, dimension(:), allocatable ldefavg
logical, dimension(:), allocatable ldefqck
logical, dimension(:), allocatable lwrtavg
logical, dimension(:), allocatable ldefadj
logical, dimension(:), allocatable ldefhis
integer, dimension(:), allocatable ntend
logical, dimension(:), allocatable ldefmod
logical, dimension(:,:), allocatable ldefnrm
logical, dimension(:), allocatable lwrtcost
integer exit_flag
logical, dimension(:), allocatable lwrtstate2d
real(dp), dimension(:,:), allocatable obc_time
real(dp) myruninterval
integer erstr
logical, dimension(:), allocatable lwrthis
logical, dimension(:), allocatable lwrttlm
real(dp), dimension(:), allocatable time
logical, dimension(:), allocatable ldeftlm
integer, dimension(:), allocatable nsff
integer, dimension(:), allocatable ntstart
integer nrun
integer, dimension(:), allocatable step_counter
real(dp), dimension(:,:), allocatable sf_time
integer, dimension(:), allocatable nbrec
logical, dimension(:), allocatable lreadfwd
integer inner
real(dp), dimension(:), allocatable initime
integer noerror
integer outer
logical, dimension(:), allocatable lreadblk
integer, dimension(:), allocatable lold
integer, dimension(:), allocatable lbout
integer, dimension(:), allocatable kstp
integer, dimension(:), allocatable nrhs
integer, dimension(:), allocatable lfinp
integer, dimension(:), allocatable lnew
integer, dimension(:), allocatable lfout
integer, dimension(:), allocatable nstp
subroutine, public normalization(ng, tile, ifac)
subroutine, public set_masks(ng, tile, model)
Definition set_masks.F:44
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
subroutine, public sum_grad(ng, tile, linp, lout)
Definition sum_grad.F:26
subroutine, public tl_balance(ng, tile, lbck, linp)
Definition tl_balance.F:59
subroutine, public tl_convolution(ng, tile, linp, lweak, ifac)
subroutine, public tl_def_his(ng, ldef)
Definition tl_def_his.F:51
subroutine, public tl_def_ini(ng)
Definition tl_def_ini.F:43
subroutine, public tl_variability(ng, tile, linp, lweak)
subroutine, public tl_wrt_ini(ng, tile, tindex, outrec)
Definition tl_wrt_ini.F:74
subroutine, public wrt_evolved(ng, tile, kout, nout)
Definition wrt_evolved.F:55
subroutine, public wrt_ini(ng, tile, tindex, outrec)
Definition wrt_ini.F:89
subroutine, public wrt_frc(ng, tile, tindex, outrec)
Definition wrt_ini.F:1011
subroutine, public biconj(ng, tile, model, lbck)
subroutine, public balance_ref(ng, tile, lbck)
subroutine set_grid(ng, model)
Definition set_grid.F:3
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3
subroutine tl_initial(ng)
Definition tl_initial.F:4
subroutine tl_main2d
Definition tl_main2d.F:429
subroutine tl_main3d(runinterval)
Definition tl_main3d.F:4