ROMS
Loading...
Searching...
No Matches
tl_rbl4dvar_roms.h
Go to the documentation of this file.
1 MODULE roms_kernel_mod
2!
3!git $Id$
4!=================================================== Andrew M. Moore ===
5! Copyright (c) 2002-2025 The ROMS Group Hernan G. Arango !
6! Licensed under a MIT/X style license !
7! See License_ROMS.md !
8!=======================================================================
9! !
10! ROMS Strong/Weak Constraint 4-Dimensional Variational Data !
11! Assimilation and its Tangent Linear Driver: Restricted !
12! B-preconditioned Lanczos (RBL4D-Var). !
13! !
14! This driver is used for the dual formulation (observation space), !
15! strong or weak constraint 4D-Var where errors may be considered !
16! in both model and observations. !
17! !
18! This driver is used for strong/weak constraint 4D-Var where errors !
19! may be considered in both model and observations. !
20! !
21! The routines in this driver control the initialization, time- !
22! stepping, and finalization of ROMS model following ESMF !
23! conventions: !
24! !
25! ROMS_initialize !
26! ROMS_run !
27! ROMS_finalize !
28! !
29! References: !
30! !
31! Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver, !
32! and J. Zavala-Garay, 2011: The Regional Ocean Modeling System !
33! (ROMS) 4-dimensional variational data assimilations systems, !
34! Part I - System overview and formulation, Prog. Oceanogr., 91, !
35! 34-49, doi:10.1016/j.pocean.2011.05.004. !
36! !
37! Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani, !
38! B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson, !
39! 2011: The Regional Ocean Modeling System (ROMS) 4-dimensional !
40! variational data assimilations systems, Part II - Performance !
41! and application to the California Current System, Prog. !
42! Oceanogr., 91, 50-73, doi:10.1016/j.pocean.2011.05.003. !
43! !
44! Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani, !
45! B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson, !
46! 2011: The Regional Ocean Modeling System (ROMS) 4-dimensional !
47! variational data assimilations systems, Part III - Observation !
48! impact and observation sensitivity in the California Current !
49! System, Prog. Oceanogr., 91, 74-94, !
50! doi:10.1016/j.pocean.2011.05.005. !
51! !
52!=======================================================================
53!
54 USE mod_param
55 USE mod_parallel
56 USE mod_arrays
57 USE mod_fourdvar
58 USE mod_iounits
59 USE mod_ncparam
60 USE mod_netcdf
61#if defined PIO_LIB && defined DISTRIBUTE
63#endif
64 USE mod_scalars
65 USE mod_stepping
66!
67#ifdef ADJUST_BOUNDARY
69#endif
71 USE mod_ocean, ONLY : initialize_ocean
72!
73 USE ad_wrt_his_mod, ONLY : ad_wrt_his
75 USE congrad_mod, ONLY : congrad
78 USE def_mod_mod, ONLY : def_mod
79 USE def_norm_mod, ONLY : def_norm
80 USE get_state_mod, ONLY : get_state
81 USE inp_par_mod, ONLY : inp_par
82#ifdef MCT_LIB
83# ifdef ATM_COUPLING
84 USE mct_coupler_mod, ONLY : initialize_ocn2atm_coupling
85# endif
86# ifdef WAV_COUPLING
87 USE mct_coupler_mod, ONLY : initialize_ocn2wav_coupling
88# endif
89#endif
93 USE tl_def_ini_mod, ONLY : tl_def_ini
95 USE wrt_ini_mod, ONLY : wrt_ini
96 USE wrt_rst_mod, ONLY : wrt_rst
97#if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
99#endif
100!
101 implicit none
102!
103 PUBLIC :: roms_initialize
104 PUBLIC :: roms_run
105 PUBLIC :: roms_finalize
106!
107 CONTAINS
108!
109 SUBROUTINE roms_initialize (first, mpiCOMM)
110!
111!=======================================================================
112! !
113! This routine allocates and initializes ROMS state variables !
114! and internal and external parameters. !
115! !
116!=======================================================================
117!
118! Imported variable declarations.
119!
120 logical, intent(inout) :: first
121!
122 integer, intent(in), optional :: mpiCOMM
123!
124! Local variable declarations.
125!
126 logical :: allocate_vars = .true.
127!
128#ifdef DISTRIBUTE
129 integer :: MyError, MySize
130#endif
131 integer :: STDrec, Tindex
132 integer :: chunk_size, ng, thread
133#ifdef _OPENMP
134 integer :: my_threadnum
135#endif
136!
137 character (len=*), parameter :: MyFile = &
138 & __FILE__//", ROMS_initialize"
139
140#ifdef DISTRIBUTE
141!
142!-----------------------------------------------------------------------
143! Set distribute-memory (mpi) world communictor.
144!-----------------------------------------------------------------------
145!
146 IF (PRESENT(mpicomm)) THEN
147 ocn_comm_world=mpicomm
148 ELSE
149 ocn_comm_world=mpi_comm_world
150 END IF
151 CALL mpi_comm_rank (ocn_comm_world, myrank, myerror)
152 CALL mpi_comm_size (ocn_comm_world, mysize, myerror)
153#endif
154!
155!-----------------------------------------------------------------------
156! On first pass, initialize model parameters a variables for all
157! nested/composed grids. Notice that the logical switch "first"
158! is used to allow multiple calls to this routine during ensemble
159! configurations.
160!-----------------------------------------------------------------------
161!
162 IF (first) THEN
163 first=.false.
164!
165! Initialize parallel control switches. These scalars switches are
166! independent from standard input parameters.
167!
169!
170! Set the ROMS standard output unit to write verbose execution info.
171! Notice that the default standard out unit in Fortran is 6.
172!
173! In some applications like coupling or disjointed mpi-communications,
174! it is advantageous to write standard output to a specific filename
175! instead of the default Fortran standard output unit 6. If that is
176! the case, it opens such formatted file for writing.
177!
178 IF (set_stdoutunit) THEN
180 set_stdoutunit=.false.
181 END IF
182!
183! Read in model tunable parameters from standard input. Allocate and
184! initialize variables in several modules after the number of nested
185! grids and dimension parameters are known.
186!
187 CALL inp_par (inlm)
188 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
189!
190! Set domain decomposition tile partition range. This range is
191! computed only once since the "first_tile" and "last_tile" values
192! are private for each parallel thread/node.
193!
194#if defined _OPENMP
195 mythread=my_threadnum()
196#elif defined DISTRIBUTE
198#else
199 mythread=0
200#endif
201 DO ng=1,ngrids
202 chunk_size=(ntilex(ng)*ntilee(ng)+numthreads-1)/numthreads
203 first_tile(ng)=mythread*chunk_size
204 last_tile(ng)=first_tile(ng)+chunk_size-1
205 END DO
206!
207! Initialize internal wall clocks. Notice that the timings does not
208! includes processing standard input because several parameters are
209! needed to allocate clock variables.
210!
211 IF (master) THEN
212 WRITE (stdout,10)
213 10 FORMAT (/,' Process Information:',/)
214 END IF
215!
216 DO ng=1,ngrids
217 DO thread=thread_range
218 CALL wclock_on (ng, inlm, 0, __line__, myfile)
219 END DO
220 END DO
221!
222! Allocate and initialize modules variables.
223!
224 CALL roms_allocate_arrays (allocate_vars)
226 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
227
228 END IF
229
230#if defined MCT_LIB && (defined ATM_COUPLING || defined WAV_COUPLING)
231!
232!-----------------------------------------------------------------------
233! Initialize coupling streams between model(s).
234!-----------------------------------------------------------------------
235!
236 DO ng=1,ngrids
237# ifdef ATM_COUPLING
238 CALL initialize_ocn2atm_coupling (ng, myrank)
239# endif
240# ifdef WAV_COUPLING
241 CALL initialize_ocn2wav_coupling (ng, myrank)
242# endif
243 END DO
244#endif
245!
246!-----------------------------------------------------------------------
247! Read in standard deviation factors for error covariance.
248!-----------------------------------------------------------------------
249!
250! Initial conditions standard deviation. They are loaded in Tindex=1
251! of the e_var(...,Tindex) state variables.
252!
253 stdrec=1
254 tindex=1
255 DO ng=1,ngrids
256 CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
257 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
258 END DO
259!
260! Model error standard deviation. They are loaded in Tindex=2
261! of the e_var(...,Tindex) state variables.
262!
263 stdrec=1
264 tindex=2
265 IF (nsa.eq.2) THEN
266 DO ng=1,ngrids
267 CALL get_state (ng, 11, 11, std(2,ng), stdrec, tindex)
268 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
269 END DO
270 END IF
271
272#ifdef ADJUST_BOUNDARY
273!
274! Open boundary conditions standard deviation.
275!
276 stdrec=1
277 tindex=1
278 DO ng=1,ngrids
279 CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
280 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
281 END DO
282#endif
283#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
284!
285! Surface forcing standard deviation.
286!
287 stdrec=1
288 tindex=1
289 DO ng=1,ngrids
290 CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
291 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
292 END DO
293#endif
294!
295 RETURN
296 END SUBROUTINE roms_initialize
297!
298 SUBROUTINE roms_run (RunInterval)
299!
300!=======================================================================
301! !
302! This routine time-steps ROMS nonlinear, tangent linear and !
303! adjoint models. !
304! !
305!=======================================================================
306!
307! Imported variable declarations
308!
309 real(dp), intent(in) :: RunInterval ! seconds
310!
311! Local variable declarations.
312!
313 logical :: Lcgini, Linner, Lposterior
314!
315 integer :: my_inner, my_outer
316 integer :: Lbck, Lini, Rec1, Rec2
317 integer :: i, lstr, ng, status, tile
318 integer :: Fcount, Mstr, NRMrec
319
320 integer, dimension(Ngrids) :: indxSave
321 integer, dimension(Ngrids) :: Nrec
322!
323 character (len=11) :: driver
324 character (len=20) :: string
325
326 character (len=*), parameter :: MyFile = &
327 & __FILE__//", ROMS_run"
328!
329!=======================================================================
330! Run model for all nested grids, if any.
331!=======================================================================
332!
333! Initialize relevant parameters.
334!
335 DO ng=1,ngrids
336#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
337 lfinp(ng)=1 ! forcing index for input
338 lfout(ng)=1 ! forcing index for output history files
339#endif
340#ifdef ADJUST_BOUNDARY
341 lbinp(ng)=1 ! boundary index for input
342 lbout(ng)=1 ! boundary index for output history files
343#endif
344 lold(ng)=1 ! old minimization time index
345 lnew(ng)=2 ! new minimization time index
346 END DO
347 lini=1 ! NLM initial conditions record in INI
348 lbck=2 ! background record in INI
349 rec1=1
350 rec2=2
351 nrun=1
352 outer=0
353 inner=0
354 erstr=1
356 driver='tl_rbl4dvar'
357!
358!-----------------------------------------------------------------------
359! Configure weak constraint RBL4D-Var algorithm.
360!-----------------------------------------------------------------------
361!
362! Initialize the switch to gather weak constraint forcing.
363!
364 DO ng=1,ngrids
365 wrtforce(ng)=.false.
366 END DO
367!
368! Initialize and set nonlinear model initial conditions.
369!
370 DO ng=1,ngrids
371 wrtnlmod(ng)=.true.
372 wrtrpmod(ng)=.false.
373 wrttlmod(ng)=.false.
374 END DO
375!
376 CALL initial
377 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
378!
379! Save nonlinear initial conditions (currently in time index 1,
380! background) into record "Lbck" of INI(ng)%name NetCDF file. The
381! record "Lbck" becomes the background state record and the record
382! "Lini" becomes current nonlinear initial conditions.
383!
384 DO ng=1,ngrids
385 ini(ng)%Rindex=1
386 fcount=ini(ng)%load
387 ini(ng)%Nrec(fcount)=1
388#ifdef DISTRIBUTE
389 CALL wrt_ini (ng, myrank, 1)
390#else
391 CALL wrt_ini (ng, -1, 1)
392#endif
393 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
394 END DO
395!
396! Set nonlinear output history file as the initial basic state
397! trajectory.
398!
399 DO ng=1,ngrids
400 ldefhis(ng)=.true.
401 lwrthis(ng)=.true.
402 WRITE (his(ng)%name,10) trim(fwd(ng)%head), outer
403 lstr=len_trim(his(ng)%name)
404 his(ng)%base=his(ng)%name(1:lstr-3)
405 END DO
406!
407!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
408! Model-error covariance normalization and stardard deviation factors.
409!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
410!
411! Compute or read in the error covariance normalization factors.
412! If computing, write out factors to NetCDF. This is an expensive
413! computation that needs to be computed only once for a particular
414! application grid and decorrelation scales.
415!
416 DO ng=1,ngrids
417 IF (any(lwrtnrm(:,ng))) THEN
418 CALL def_norm (ng, inlm, 1)
419 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
420
421 IF (nsa.eq.2) THEN
422 CALL def_norm (ng, inlm, 2)
423 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
424 END IF
425
426#ifdef ADJUST_BOUNDARY
427 CALL def_norm (ng, inlm, 3)
428 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
429#endif
430
431#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
432 CALL def_norm (ng, inlm, 4)
433 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
434#endif
435
436 DO tile=first_tile(ng),last_tile(ng),+1
437 CALL normalization (ng, tile, 2)
438 END DO
439 ldefnrm(1:4,ng)=.false.
440 lwrtnrm(1:4,ng)=.false.
441 ELSE
442 nrmrec=1
443 CALL get_state (ng, 14, 14, nrm(1,ng), nrmrec, 1)
444 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
445
446 IF (nsa.eq.2) THEN
447 CALL get_state (ng, 15, 15, nrm(2,ng), nrmrec, 2)
448 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
449 END IF
450
451#ifdef ADJUST_BOUNDARY
452 CALL get_state (ng, 16, 16, nrm(3,ng), nrmrec, 1)
453 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
454#endif
455
456#if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
457 CALL get_state (ng, 17, 17, nrm(4,ng), nrmrec, 1)
458 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
459#endif
460 END IF
461 END DO
462!
463! Define tangent linear initial conditions file.
464!
465 DO ng=1,ngrids
466 ldefitl(ng)=.true.
467 CALL tl_def_ini (ng)
468 ldefitl(ng)=.false.
469 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
470 END DO
471!
472! Define impulse forcing NetCDF file.
473!
474 DO ng=1,ngrids
475 ldeftlf(ng)=.true.
476 CALL def_impulse (ng)
477 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
478 END DO
479!
480! Define output 4DVAR NetCDF file containing all processed data
481! at observation locations.
482!
483 DO ng=1,ngrids
484 ldefmod(ng)=.true.
485 CALL def_mod (ng)
486 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
487 END DO
488!
489!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
490! Run nonlinear model and compute background state trajectory, X_n-1(t)
491! and the background values at the observation points and times. It
492! processes and writes the observations accept/reject flag (ObsScale)
493! once to allow background quality control, if any.
494!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
495!
496 DO ng=1,ngrids
497 wrtobsscale(ng)=.true.
498 sporadicimpulse(ng)=.false.
499 frequentimpulse(ng)=.false.
500 IF (master) THEN
501 WRITE (stdout,20) 'NL', ng, ntstart(ng), ntend(ng)
502 END IF
503 END DO
504!
505#ifdef SOLVE3D
506 CALL main3d (runinterval)
507#else
508 CALL main2d (runinterval)
509#endif
510 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
511
512 DO ng=1,ngrids
513 wrtnlmod(ng)=.false.
514 wrtobsscale(ng)=.false.
515 END DO
516!
517! Set structure for the nonlinear forward trajectory to be processed
518! by the tangent linear and adjoint models. Also, set switches to
519! process the FWD structure in routine "check_multifile". Notice that
520! it is possible to split solution into multiple NetCDF files to reduce
521! their size.
522!
523 CALL edit_multifile ('HIS2FWD')
524 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
525 DO ng=1,ngrids
526 lreadfwd(ng)=.true.
527 END DO
528
529#ifdef FORWARD_FLUXES
530!
531! Set the BLK structure to contain the nonlinear model surface fluxes
532! needed by the tangent linear and adjoint models. Also, set switches
533! to process that structure in routine "check_multifile". Notice that
534! it is possible to split the solution into multiple NetCDF files to
535! reduce their size.
536!
537! The switch LreadFRC is deactivated because all the atmospheric
538! forcing, including shortwave radiation, is read from the NLM
539! surface fluxes or is assigned during ESM coupling. Such fluxes
540! are available from the QCK structure. There is no need for reading
541! and processing from the FRC structure input forcing-files.
542!
543 CALL edit_multifile ('QCK2BLK')
544 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
545 DO ng=1,ngrids
546 lreadblk(ng)=.true.
547 lreadfrc(ng)=.false.
548 lreadqck(ng)=.false.
549 END DO
550#endif
551!
552! Report data penalty function.
553!
554 DO ng=1,ngrids
555 IF (master) THEN
556 DO i=0,nobsvar(ng)
557 IF (i.eq.0) THEN
558 string='Total'
559 ELSE
560 string=obsname(i)
561 END IF
562 IF (fourdvar(ng)%NLPenalty(i).ne.0.0_r8) THEN
563 WRITE (stdout,30) outer, inner, 'NLM', &
564 & fourdvar(ng)%NLPenalty(i), &
565 & trim(string)
566 END IF
567 END DO
568 END IF
569!
570! Write out initial data penalty function to NetCDF file.
571!
572 sourcefile=myfile
573 SELECT CASE (dav(ng)%IOtype)
574 CASE (io_nf90)
575 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
576 & 'NL_iDataPenalty', &
577 & fourdvar(ng)%NLPenalty(0:), &
578 & (/1/), (/nobsvar(ng)+1/), &
579 & ncid = dav(ng)%ncid)
580
581#if defined PIO_LIB && defined DISTRIBUTE
582 CASE (io_pio)
583 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
584 & 'NL_iDataPenalty', &
585 & fourdvar(ng)%NLPenalty(0:), &
586 & (/1/), (/nobsvar(ng)+1/), &
587 & piofile = dav(ng)%pioFile)
588#endif
589 END SELECT
590 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
591!
592! Clean penalty array before next run of NL model.
593!
594 fourdvar(ng)%NLPenalty=0.0_r8
595 END DO
596!
597!-----------------------------------------------------------------------
598! Solve the system (following Courtier, 1997):
599!
600! (H M_n B (M_n)' H' + Cobs) * w_n = d_n
601!
602! d_n = yo - H * Xb_n
603!
604! where M_n is the tangent linear model matrix, Cobs is the
605! observation-error covariance, B is the background error covariance
606! and dx_n=B M' H' w_n is the analysis increment so that Xa=Xb+dx_n.
607! d_n is the misfit between observations (yo) and model (H * Xb_n),
608! and H is the linearized observation operator.
609!
610! Here, _n denotes a sequence of outer-loop estimates.
611!
612! The system does not need to be solved explicitly by inverting the
613! symmetric matrix, P_n:
614!
615! P_n = H M_n B (M_n)' H' + Cobs
616!
617! but by computing the action of P_n on any vector PSI, such that
618!
619! P_n * PSI = H M_n B (M_n)' H' * PSI + Cobs * PSI
620!
621! The (H M_n B (M_n)' H') matrix is not explicitly computed but
622! evaluated by one integration backward of the adjoint model and
623! one integration forward of the tangent linear model for any
624! forcing vector PSI.
625!
626! A preconditioned conjugate gradient algorithm is used to compute
627! an approximation PSI for w_n.
628!
629!-----------------------------------------------------------------------
630!
631 outer_loop1 : DO my_outer=1,nouter
632 outer=my_outer
633 inner=0
634!
635! Set basic state trajectory (X_n-1) file to previous outer loop file
636! (outer-1).
637!
638 DO ng=1,ngrids
639 WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
640 lstr=len_trim(his(ng)%name)
641 his(ng)%base=his(ng)%name(1:lstr-3)
642 END DO
643!
644! Clear tangent linear forcing arrays before entering inner-loop.
645! This is very important since these arrays are non-zero and must
646! be zero when running the tangent linear model.
647!
648 DO ng=1,ngrids
649 DO tile=first_tile(ng),last_tile(ng),+1
650 CALL initialize_forces (ng, tile, itlm)
651# ifdef ADJUST_BOUNDARY
652 CALL initialize_boundary (ng, tile, itlm)
653# endif
654 END DO
655 END DO
656
657#if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
658!
659! Compute the reference zeta and biconjugate gradient arrays
660! required for the balance of free surface.
661!
662 IF (balance(isfsur)) THEN
663 DO ng=1,ngrids
664 CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
665 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
666!
667 DO tile=first_tile(ng),last_tile(ng),+1
668 CALL balance_ref (ng, tile, lini)
669 CALL biconj (ng, tile, inlm, lini)
670 END DO
671 wrtzetaref(ng)=.true.
672 END DO
673 END IF
674#endif
675!
676 inner_loop1 : DO my_inner=0,ninner
677 inner=my_inner
678!
679! Initialize conjugate gradient algorithm depending on hot start or
680! outer loop index.
681!
682 IF (inner.eq.0) THEN
683 lcgini=.true.
684 DO ng=1,ngrids
685 CALL congrad (ng, irpm, outer, inner, ninner, lcgini)
686 END DO
687 END IF
688!
689! If initialization step, skip the inner-loop computations.
690!
691 linner=.false.
692 IF ((inner.ne.0).or.(nrun.ne.1)) THEN
693 IF (((inner.eq.0).and.lhotstart).or.(inner.ne.0)) THEN
694 linner=.true.
695 END IF
696 END IF
697!
698! Start inner loop computations.
699!
700 inner_compute1 : IF (linner) THEN
701!
702!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
703! Integrate adjoint model forced with any vector PSI at the observation
704! locations and generate adjoint trajectory, Lambda_n(t).
705!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
706!
707! Initialize the adjoint model from rest.
708!
709 DO ng=1,ngrids
710 CALL ad_initial (ng)
712 & __line__, myfile)) RETURN
713 wrtmisfit(ng)=.false.
714 END DO
715!
716! Set adjoint history NetCDF parameters. Define adjoint history
717! file only once to avoid opening too many files.
718!
719 DO ng=1,ngrids
720 wrtforce(ng)=.true.
721 IF (nrun.gt.1) ldefadj(ng)=.false.
722 fcount=adm(ng)%load
723 adm(ng)%Nrec(fcount)=0
724 adm(ng)%Rindex=0
725 END DO
726!
727! Time-step adjoint model backwards forced with current PSI vector.
728!
729 DO ng=1,ngrids
730 IF (master) THEN
731 WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
732 END IF
733 END DO
734!
735#ifdef SOLVE3D
736 CALL ad_main3d (runinterval)
737#else
738 CALL ad_main2d (runinterval)
739#endif
740 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
741!
742! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
743! record into the adjoint history file. Note that the weak-constraint
744! forcing is delayed by nADJ time-steps.
745!
746 DO ng=1,ngrids
747#ifdef DISTRIBUTE
748 CALL ad_wrt_his (ng, myrank)
749#else
750 CALL ad_wrt_his (ng, -1)
751#endif
753 & __line__, myfile)) RETURN
754 END DO
755!
756! Write out adjoint initial condition record into the adjoint
757! history file.
758!
759 DO ng=1,ngrids
760 wrtforce(ng)=.false.
761#ifdef DISTRIBUTE
762 CALL ad_wrt_his (ng, myrank)
763#else
764 CALL ad_wrt_his (ng, -1)
765#endif
767 & __line__, myfile)) RETURN
768 END DO
769!
770! Convolve adjoint trajectory with error covariances.
771!
772 lposterior=.false.
773 CALL error_covariance (itlm, driver, outer, inner, &
774 & lbck, lini, lold, lnew, &
775 & rec1, rec2, lposterior)
776 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
777!
778! Convert the current adjoint solution in ADM(ng)%name to impulse
779! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
780! To facilitate the forcing to the TLM and RPM, the forcing is
781! processed and written in increasing time coordinates (recall that
782! the adjoint solution in ADM(ng)%name is backwards in time).
783!
784 IF (master) THEN
785 WRITE (stdout,40) outer, inner
786 END IF
787 DO ng=1,ngrids
788 tlf(ng)%Rindex=0
789#ifdef DISTRIBUTE
790 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
791#else
792 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
793#endif
795 & __line__, myfile)) RETURN
796 END DO
797!
798!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
799! Integrate tangent linear model forced by the convolved adjoint
800! trajectory (impulse forcing) to compute R_n * PSI at observation
801! points.
802!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
803!
804! Initialize tangent linear model from initial impulse which is now
805! stored in file ITL(ng)%name.
806!
807 DO ng=1,ngrids
808 wrtnlmod(ng)=.false.
809 wrttlmod(ng)=.true.
810 END DO
811!
812! If weak constraint, the impulses are time-interpolated at each
813! time-steps.
814!
815 DO ng=1,ngrids
816 IF (frcrec(ng).gt.3) THEN
817 frequentimpulse(ng)=.true.
818 END IF
819 END DO
820!
821! Initialize tangent linear model from ITL(ng)%name, record Rec1.
822!
823 DO ng=1,ngrids
824 itl(ng)%Rindex=rec1
825 CALL tl_initial (ng)
827 & __line__, myfile)) RETURN
828 END DO
829!
830! Activate switch to write out initial misfit between model and
831! observations.
832!
833 IF ((outer.eq.1).and.(inner.eq.1)) THEN
834 DO ng=1,ngrids
835 wrtmisfit(ng)=.true.
836 END DO
837 END IF
838!
839! Set tangent linear history NetCDF parameters. Define tangent linear
840! history file at the beggining of each inner loop to avoid opening
841! too many NetCDF files.
842!
843 IF (inner.gt.1) ldeftlm(ng)=.false.
844 fcount=tlm(ng)%load
845 tlm(ng)%Nrec(fcount)=0
846 tlm(ng)%Rindex=0
847!
848! Run tangent linear model forward and force with convolved adjoint
849! trajectory impulses. Compute (H M B M' H')_n * PSI at observation
850! points which are used in the conjugate gradient algorithm.
851!
852 DO ng=1,ngrids
853 IF (master) THEN
854 WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
855 END IF
856 END DO
857!
858#ifdef SOLVE3D
859 CALL tl_main3d (runinterval)
860#else
861 CALL tl_main2d (runinterval)
862#endif
863 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
864
865 DO ng=1,ngrids
866 wrtnlmod(ng)=.false.
867 wrttlmod(ng)=.false.
868 END DO
869!
870!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
871! Use conjugate gradient algorithm to find a better approximation
872! PSI to coefficients Beta_n.
873!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
874!
875 nrun=nrun+1
876 DO ng=1,ngrids
877 lcgini=.false.
878 CALL congrad (ng, itlm, outer, inner, ninner, lcgini)
880 & __line__, myfile)) RETURN
881 END DO
882
883 END IF inner_compute1
884
885 END DO inner_loop1
886!
887! Close tangent linear NetCDF file.
888!
889 DO ng=1,ngrids
890 status=nf90_close(tlm(ng)%ncid)
891 tlm(ng)%ncid=-1
892 END DO
893!
894!-----------------------------------------------------------------------
895! Once the w_n, have been approximated with sufficient accuracy,
896! compute estimates of Lambda_n and Xhat_n by carrying out one
897! backward intergration of the adjoint model and one forward
898! itegration of the nonlinear model.
899!-----------------------------------------------------------------------
900!
901! Initialize the adjoint model always from rest.
902!
903 DO ng=1,ngrids
904 CALL ad_initial (ng)
905 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
906 END DO
907!
908! Set adjoint history NetCDF parameters. Define adjoint history
909! file one to avoid opening to many files.
910!
911 DO ng=1,ngrids
912 wrtforce(ng)=.true.
913 IF (nrun.gt.1) ldefadj(ng)=.false.
914 fcount=adm(ng)%load
915 adm(ng)%Nrec(fcount)=0
916 adm(ng)%Rindex=0
917 END DO
918!
919! Time-step adjoint model backwards forced with estimated coefficients,
920! Beta_n.
921!
922 DO ng=1,ngrids
923 IF (master) THEN
924 WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
925 END IF
926 END DO
927!
928#ifdef SOLVE3D
929 CALL ad_main3d (runinterval)
930#else
931 CALL ad_main2d (runinterval)
932#endif
933 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
934!
935! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
936! record into the adjoint history file. Note that the weak-constraint
937! forcing is delayed by nADJ time-steps.
938!
939 DO ng=1,ngrids
940#ifdef DISTRIBUTE
941 CALL ad_wrt_his (ng, myrank)
942#else
943 CALL ad_wrt_his (ng, -1)
944#endif
945 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
946 END DO
947!
948! Write out adjoint initial condition record into the adjoint
949! history file.
950!
951 DO ng=1,ngrids
952 wrtforce(ng)=.false.
953#ifdef DISTRIBUTE
954 CALL ad_wrt_his (ng, myrank)
955#else
956 CALL ad_wrt_his (ng, -1)
957#endif
958 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
959 END DO
960!
961! Convolve adjoint trajectory with error covariances.
962!
963 lposterior=.false.
964 CALL error_covariance (inlm, driver, outer, inner, &
965 & lbck, lini, lold, lnew, &
966 & rec1, rec2, lposterior)
967 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
968!
969! Convert the current adjoint solution in ADM(ng)%name to impulse
970! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
971! To facilitate the forcing to the TLM and RPM, the forcing is
972! processed and written in increasing time coordinates (recall that
973! the adjoint solution in ADM(ng)%name is backwards in time).
974!
975 IF (master) THEN
976 WRITE (stdout,40) outer, inner
977 END IF
978 DO ng=1,ngrids
979 tlf(ng)%Rindex=0
980#ifdef DISTRIBUTE
981 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
982#else
983 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
984#endif
985 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
986 END DO
987!
988!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
989! Run nonlinear model and compute a "new estimate" of the state
990! trajectory, X_n(t).
991!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
992!
993! Set new basic state trajectory for next outer loop.
994!
995 DO ng=1,ngrids
996 idefhis(ng)=-1
997 ldefhis(ng)=.true.
998 lwrthis(ng)=.true.
999 wrtnlmod(ng)=.true.
1000 wrttlmod(ng)=.false.
1001#ifdef FORWARD_FLUXES
1002 lreadblk(ng)=.false.
1003#endif
1004 lreadfwd(ng)=.false.
1005 WRITE (his(ng)%name,10) trim(fwd(ng)%head), outer
1006 lstr=len_trim(his(ng)%name)
1007 his(ng)%base=his(ng)%name(1:lstr-3)
1008 END DO
1009!
1010! If weak constraint, the impulses are time-interpolated at each
1011! time-steps.
1012!
1013 DO ng=1,ngrids
1014 IF (frcrec(ng).gt.3) THEN
1015 frequentimpulse(ng)=.true.
1016 END IF
1017 END DO
1018!
1019! Clear tangent arrays before running nonlinear model (important).
1020!
1021 DO ng=1,ngrids
1022 DO tile=first_tile(ng),last_tile(ng),+1
1023 CALL initialize_ocean (ng, tile, itlm)
1024 CALL initialize_forces (ng, tile, itlm)
1025# ifdef ADJUST_BOUNDARY
1026 CALL initialize_boundary (ng, tile, itlm)
1027# endif
1028 END DO
1029 END DO
1030!
1031! Initialize nonlinear model INI(ng)%name file, record outer+2.
1032! Notice that NetCDF record index counter is saved because this
1033! counter is used to write initial conditions.
1034!
1035 DO ng=1,ngrids
1036 indxsave(ng)=ini(ng)%Rindex
1037 ini(ng)%Rindex=outer+2
1038 END DO
1039!
1040 CALL initial
1041 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1042
1043 DO ng=1,ngrids
1044 ini(ng)%Rindex=indxsave(ng)
1045 END DO
1046!
1047! Activate switch to write out final misfit between model and
1048! observations.
1049!
1050 IF (outer.eq.nouter) THEN
1051 DO ng=1,ngrids
1052 wrtmisfit(ng)=.true.
1053 END DO
1054 END IF
1055!
1056! Run nonlinear forced by convolved adjoint trajectory impulses and
1057! compute new basic state trajectory X_n.
1058!
1059 DO ng=1,ngrids
1060 IF (master) THEN
1061 WRITE (stdout,20) 'NL', ng, ntstart(ng), ntend(ng)
1062 END IF
1063 END DO
1064!
1065#ifdef SOLVE3D
1066 CALL main3d (runinterval)
1067#else
1068 CALL main2d (runinterval)
1069#endif
1070 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1071
1072 DO ng=1,ngrids
1073 wrtnlmod(ng)=.false.
1074 wrttlmod(ng)=.false.
1075 END DO
1076
1077! Set structure for the nonlinear forward trajectory to be processed
1078! by the tangent linear and adjoint models. Also, set switches to
1079! process the FWD structure in routine "check_multifile". Notice that
1080! it is possible to split solution into multiple NetCDF files to reduce
1081! their size.
1082!
1083 CALL edit_multifile ('HIS2FWD')
1084 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1085 DO ng=1,ngrids
1086 lreadfwd(ng)=.true.
1087 END DO
1088
1089#ifdef FORWARD_FLUXES
1090!
1091! Set the BLK structure to contain the nonlinear model surface fluxes
1092! needed by the tangent linear and adjoint models. Also, set switches
1093! to process that structure in routine "check_multifile". Notice that
1094! it is possible to split the solution into multiple NetCDF files to
1095! reduce their size.
1096!
1097 CALL edit_multifile ('QCK2BLK')
1098 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1099 DO ng=1,ngrids
1100 lreadblk(ng)=.true.
1101 lreadfrc(ng)=.false.
1102 lreadqck(ng)=.false.
1103 END DO
1104#endif
1105!
1106! Report data penalty function.
1107!
1108 DO ng=1,ngrids
1109 IF (master) THEN
1110 DO i=0,nobsvar(ng)
1111 IF (i.eq.0) THEN
1112 string='Total'
1113 ELSE
1114 string=obsname(i)
1115 END IF
1116 IF (fourdvar(ng)%NLPenalty(i).ne.0.0_r8) THEN
1117 WRITE (stdout,30) outer, inner, 'NLM', &
1118 & fourdvar(ng)%NLPenalty(i), &
1119 & trim(string)
1120 END IF
1121 END DO
1122 END IF
1123!
1124! Write out final data penalty function to NetCDF file.
1125!
1126 sourcefile=myfile
1127 SELECT CASE (dav(ng)%IOtype)
1128 CASE (io_nf90)
1129 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1130 & 'NL_fDataPenalty', &
1131 & fourdvar(ng)%NLPenalty(0:), &
1132 & (/1,outer/), &
1133 & (/nobsvar(ng)+1,1/), &
1134 & ncid = dav(ng)%ncid)
1135
1136#if defined PIO_LIB && defined DISTRIBUTE
1137 CASE (io_pio)
1138 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
1139 & 'NL_fDataPenalty', &
1140 & fourdvar(ng)%NLPenalty(0:), &
1141 & (/1,outer/), &
1142 & (/nobsvar(ng)+1,1/), &
1143 & piofile = dav(ng)%pioFile)
1144#endif
1145 END SELECT
1146 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1147!
1148! Clean penalty array before next run of NL model.
1149!
1150 fourdvar(ng)%NLPenalty=0.0_r8
1151 END DO
1152!
1153! Close current forward NetCDF file.
1154!
1155 DO ng=1,ngrids
1156 CALL close_file (ng, inlm, fwd(ng))
1157 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1158!
1159 IF (his(ng)%IOtype.eq.io_nf90) THEN
1160 his(ng)%ncid=-1
1161#if defined PIO_LIB && defined DISTRIBUTE
1162 ELSE IF (his(ng)%IOtype.eq.io_pio) THEN
1163 his(ng)%pioFile%fh=-1
1164#endif
1165 END IF
1166 END DO
1167
1168 END DO outer_loop1
1169!!
1170!! Compute and report model-observation comparison statistics.
1171!!
1172!! DO ng=1,Ngrids
1173!! CALL stats_modobs (ng)
1174!! IF (FoundError(exit_flag, NoError, __LINE__, MyFile)) RETURN
1175!! END DO
1176!
1177!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1178!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1179! Tangent linear of RBL4D-Var
1180!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1181!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1182!
1183! Read observations.
1184!
1185 DO ng=1,ngrids
1186 mstr=nstrobs(ng)
1187
1188 SELECT CASE (obs(ng)%IOtype)
1189 CASE (io_nf90)
1190 CALL netcdf_get_fvar (ng, itlm, obs(ng)%name, &
1191 & 'obs_value', &
1192 & tl_obsval(mstr:), &
1193 & start = (/nstrobs(ng)/), &
1194 & total = (/nobs(ng)/))
1195
1196#if defined PIO_LIB && defined DISTRIBUTE
1197 CASE (io_pio)
1198 CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name, &
1199 & 'obs_value', &
1200 & tl_obsval(mstr:), &
1201 & start = (/nstrobs(ng)/), &
1202 & total = (/nobs(ng)/))
1203#endif
1204 END SELECT
1205 END DO
1206!
1207! Reset value of "Nrun" and clear "TLmodVal".
1208!
1209 nrun=1
1210 tlmodval=0.0_r8
1211!
1212! Start outer loop.
1213!
1214 outer_loop2 : DO my_outer=1,nouter
1215 outer=my_outer
1216 inner=0
1217!
1218! Set basic state trajectory (X_n-1) file to previous outer loop file
1219! (outer-1).
1220!
1221 DO ng=1,ngrids
1222 WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer-1
1223 lstr=len_trim(fwd(ng)%name)
1224 fwd(ng)%base=fwd(ng)%name(1:lstr-3)
1225 END DO
1226!
1227! Clear tangent linear forcing arrays before entering inner-loop.
1228! This is very important since these arrays are non-zero and must
1229! be zero when running the tangent linear model.
1230!
1231 DO ng=1,ngrids
1232 DO tile=first_tile(ng),last_tile(ng),+1
1233 CALL initialize_forces (ng, tile, itlm)
1234 END DO
1235 END DO
1236!
1237 inner_loop2 : DO my_inner=0,ninner
1238 inner=my_inner
1239
1240 IF (master) THEN
1241 WRITE (stdout,50) 'Tangent linear of', &
1242 & uppercase('rbl4dvar'), outer, inner
1243 END IF
1244!
1245! Initialize tangent linear conjugate gradient algorithm depending on
1246! hot start or outer loop index.
1247!
1248 IF (inner.eq.0) THEN
1249 lcgini=.true.
1250 DO ng=1,ngrids
1251 CALL tl_congrad (ng, irpm, outer, inner, ninner, lcgini)
1252 END DO
1253 END IF
1254!
1255! If initialization step, skip the inner-loop computations.
1256!
1257 linner=.false.
1258 IF ((inner.ne.0).or.(nrun.ne.1)) THEN
1259 IF (((inner.eq.0).and.lhotstart).or.(inner.ne.0)) THEN
1260 linner=.true.
1261 END IF
1262 END IF
1263!
1264! Start inner loop computations.
1265!
1266 inner_compute2 : IF (linner) THEN
1267!
1268!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1269! Integrate adjoint model forced with any vector PSI at the observation
1270! locations and generate adjoint trajectory, Lambda_n(t).
1271!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1272!
1273! Initialize the adjoint model from rest.
1274!
1275 DO ng=1,ngrids
1276 CALL ad_initial (ng)
1278 & __line__, myfile)) RETURN
1279 wrtmisfit(ng)=.false.
1280 END DO
1281!
1282! Set adjoint history NetCDF parameters. Define adjoint history
1283! file only once to avoid opening too many files.
1284!
1285 DO ng=1,ngrids
1286 wrtforce(ng)=.true.
1287 IF (nrun.gt.1) ldefadj(ng)=.false.
1288 fcount=adm(ng)%load
1289 adm(ng)%Nrec(fcount)=0
1290 adm(ng)%Rindex=0
1291 END DO
1292!
1293! Time-step adjoint model backwards forced with current PSI vector.
1294!
1295 DO ng=1,ngrids
1296 IF (master) THEN
1297 WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
1298 END IF
1299 END DO
1300!
1301#ifdef SOLVE3D
1302 CALL ad_main3d (runinterval)
1303#else
1304 CALL ad_main2d (runinterval)
1305#endif
1306 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1307!
1308! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
1309! record into the adjoint history file. Note that the weak-constraint
1310! forcing is delayed by nADJ time-steps.
1311!
1312 DO ng=1,ngrids
1313#ifdef DISTRIBUTE
1314 CALL ad_wrt_his (ng, myrank)
1315#else
1316 CALL ad_wrt_his (ng, -1)
1317#endif
1319 & __line__, myfile)) RETURN
1320 END DO
1321!
1322! Write out adjoint initial condition record into the adjoint
1323! history file.
1324!
1325 DO ng=1,ngrids
1326 wrtforce(ng)=.false.
1327#ifdef DISTRIBUTE
1328 CALL ad_wrt_his (ng, myrank)
1329#else
1330 CALL ad_wrt_his (ng, -1)
1331#endif
1333 & __line__, myfile)) RETURN
1334 END DO
1335!
1336! Convolve adjoint trajectory with error covariances.
1337!
1338 lposterior=.false.
1339 CALL error_covariance (itlm, driver, outer, inner, &
1340 & lbck, lini, lold, lnew, &
1341 & rec1, rec2, lposterior)
1342 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1343!
1344! Convert the current adjoint solution in ADM(ng)%name to impulse
1345! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
1346! To facilitate the forcing to the TLM and RPM, the forcing is
1347! processed and written in increasing time coordinates (recall that
1348! the adjoint solution in ADM(ng)%name is backwards in time).
1349!
1350 IF (master) THEN
1351 WRITE (stdout,40) outer, inner
1352 END IF
1353 DO ng=1,ngrids
1354 tlf(ng)%Rindex=0
1355#ifdef DISTRIBUTE
1356 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1357#else
1358 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1359#endif
1361 & __line__, myfile)) RETURN
1362 END DO
1363!
1364!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1365! Integrate tangent linear model forced by the convolved adjoint
1366! trajectory (impulse forcing) to compute R_n * PSI at observation
1367! points.
1368!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1369!
1370! Initialize tangent linear model from initial impulse which is now
1371! stored in file ITL(ng)%name.
1372!
1373 DO ng=1,ngrids
1374 wrtnlmod(ng)=.false.
1375 wrttlmod(ng)=.true.
1376 END DO
1377!
1378! If weak constraint, the impulses are time-interpolated at each
1379! time-steps.
1380!
1381 DO ng=1,ngrids
1382 IF (frcrec(ng).gt.3) THEN
1383 frequentimpulse(ng)=.true.
1384 END IF
1385 END DO
1386!
1387! Initialize tangent linear model from ITL(ng)%name, record Rec1.
1388!
1389 DO ng=1,ngrids
1390 itl(ng)%Rindex=rec1
1391 CALL tl_initial (ng)
1393 & __line__, myfile)) RETURN
1394 END DO
1395!
1396! Activate switch to write out initial misfit between model and
1397! observations.
1398!
1399 IF ((outer.eq.1).and.(inner.eq.1)) THEN
1400 DO ng=1,ngrids
1401 wrtmisfit(ng)=.true.
1402 END DO
1403 END IF
1404!
1405! Set tangent linear history NetCDF parameters. Define tangent linear
1406! history file at the beggining of each inner loop to avoid opening
1407! too many NetCDF files.
1408!
1409 DO ng=1,ngrids
1410 IF (inner.gt.1) ldeftlm(ng)=.false.
1411 fcount=tlm(ng)%load
1412 tlm(ng)%Nrec(fcount)=0
1413 tlm(ng)%Rindex=0
1414 END DO
1415!
1416! Run tangent linear model forward and force with convolved adjoint
1417! trajectory impulses. Compute (H M B M' H')_n * PSI at observation
1418! points which are used in the conjugate gradient algorithm.
1419!
1420 DO ng=1,ngrids
1421 IF (master) THEN
1422 WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
1423 END IF
1424 END DO
1425!
1426#ifdef SOLVE3D
1427 CALL tl_main3d (runinterval)
1428#else
1429 CALL tl_main2d (runinterval)
1430#endif
1431 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1432
1433 DO ng=1,ngrids
1434 wrtnlmod(ng)=.false.
1435 wrttlmod(ng)=.false.
1436 END DO
1437!
1438!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1439! Use conjugate gradient algorithm to find a better approximation
1440! PSI to coefficients Beta_n.
1441!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1442!
1443 nrun=nrun+1
1444 lcgini=.false.
1445 DO ng=1,ngrids
1446 CALL tl_congrad (ng, itlm, outer, inner, ninner, lcgini)
1448 & __line__, myfile)) RETURN
1449 END DO
1450
1451 END IF inner_compute2
1452
1453 END DO inner_loop2
1454!
1455!-----------------------------------------------------------------------
1456! Once the w_n, have been approximated with sufficient accuracy,
1457! compute estimates of Lambda_n and Xhat_n by carrying out one
1458! backward intergration of the adjoint model and one forward
1459! itegration of the nonlinear model.
1460!-----------------------------------------------------------------------
1461!
1462! Initialize the adjoint model always from rest.
1463!
1464 DO ng=1,ngrids
1465 CALL ad_initial (ng)
1466 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1467 END DO
1468!
1469! Set adjoint history NetCDF parameters. Define adjoint history
1470! file one to avoid opening to many files.
1471!
1472 DO ng=1,ngrids
1473 wrtforce(ng)=.true.
1474 IF (nrun.gt.1) ldefadj(ng)=.false.
1475 fcount=adm(ng)%load
1476 adm(ng)%Nrec(fcount)=0
1477 adm(ng)%Rindex=0
1478 END DO
1479!
1480! Time-step adjoint model backwards forced with estimated representer
1481! coefficients, Beta_n.
1482!
1483 DO ng=1,ngrids
1484 IF (master) THEN
1485 WRITE (stdout,20) 'AD', ng, ntstart(ng), ntend(ng)
1486 END IF
1487 END DO
1488!
1489#ifdef SOLVE3D
1490 CALL ad_main3d (runinterval)
1491#else
1492 CALL ad_main2d (runinterval)
1493#endif
1494 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1495!
1496! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
1497! record into the adjoint history file. Note that the weak-constraint
1498! forcing is delayed by nADJ time-steps.
1499!
1500 DO ng=1,ngrids
1501#ifdef DISTRIBUTE
1502 CALL ad_wrt_his (ng, myrank)
1503#else
1504 CALL ad_wrt_his (ng, -1)
1505#endif
1506 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1507 END DO
1508!
1509! Write out adjoint initial condition record into the adjoint
1510! history file.
1511!
1512 DO ng=1,ngrids
1513 wrtforce(ng)=.false.
1514#ifdef DISTRIBUTE
1515 CALL ad_wrt_his (ng, myrank)
1516#else
1517 CALL ad_wrt_his (ng, -1)
1518#endif
1519 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1520 END DO
1521!
1522! Convolve adjoint trajectory with error covariances.
1523!
1524 lposterior=.false.
1525 CALL error_covariance (itlm, driver, outer, inner, &
1526 & lbck, lini, lold, lnew, &
1527 & rec1, rec2, lposterior)
1528 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1529!
1530! Convert the current adjoint solution in ADM(ng)%name to impulse
1531! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
1532! To facilitate the forcing to the TLM and RPM, the forcing is
1533! processed and written in increasing time coordinates (recall that
1534! the adjoint solution in ADM(ng)%name is backwards in time).
1535!
1536 DO ng=1,ngrids
1537 tlf(ng)%Rindex=0
1538#ifdef DISTRIBUTE
1539 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1540#else
1541 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1542#endif
1543 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1544 END DO
1545!
1546!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1547! Integrate tangent linear model forced by the convolved adjoint
1548! trajectory (impulse forcing) to compute R_n * PSI at observation
1549! points.
1550!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1551!
1552 DO ng=1,ngrids
1553 WRITE (fwd(ng)%name,10) trim(fwd(ng)%head), outer
1554 lstr=len_trim(fwd(ng)%name)
1555 fwd(ng)%base=fwd(ng)%name(1:lstr-3)
1556 END DO
1557!
1558! Initialize tangent linear model from initial impulse which is now
1559! stored in file ITL(ng)%name.
1560!
1561 DO ng=1,ngrids
1562 wrtnlmod(ng)=.false.
1563 wrttlmod(ng)=.true.
1564 ldeftlm(ng)=.true.
1565 lwrttlm(ng)=.true.
1566 END DO
1567!
1568! If weak constraint, the impulses are time-interpolated at each
1569! time-steps.
1570!
1571 DO ng=1,ngrids
1572 IF (frcrec(ng).gt.3) THEN
1573 frequentimpulse(ng)=.true.
1574 END IF
1575 END DO
1576!
1577! Initialize tangent linear model from ITL(ng)%name, record Rec1.
1578!
1579 DO ng=1,ngrids
1580 itl(ng)%Rindex=rec1
1581 CALL tl_initial (ng)
1582 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1583 END DO
1584!
1585! Set tangent linear history NetCDF parameters. Define tangent linear
1586! history file at the beggining of each inner loop to avoid opening
1587! too many NetCDF files.
1588!
1589 DO ng=1,ngrids
1590!! AMM IF (inner.gt.1) LdefTLM(ng)=.FALSE.
1591 fcount=tlm(ng)%load
1592 tlm(ng)%Nrec(fcount)=0
1593 tlm(ng)%Rindex=0
1594 END DO
1595!
1596! Run tangent linear model forward and force with convolved adjoint
1597! trajectory impulses. Compute (HMBM'H')_n * PSI at observation points
1598! which are used in the conjugate gradient algorithm.
1599!
1600 DO ng=1,ngrids
1601 IF (master) THEN
1602 WRITE (stdout,20) 'TL', ng, ntstart(ng), ntend(ng)
1603 END IF
1604 END DO
1605!
1606#ifdef SOLVE3D
1607 CALL tl_main3d (runinterval)
1608#else
1609 CALL tl_main2d (runinterval)
1610#endif
1611 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1612
1613 DO ng=1,ngrids
1614 wrtnlmod(ng)=.false.
1615 wrttlmod(ng)=.false.
1616 END DO
1617!
1618! Close tangent linear NetCDF file.
1619!
1620 sourcefile=myfile
1621 DO ng=1,ngrids
1622 CALL close_file (ng, itlm, tlm(ng))
1623 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1624 END DO
1625!
1626! Close current forward NetCDF file.
1627!
1628 DO ng=1,ngrids
1629 CALL close_file (ng, inlm, fwd(ng))
1630 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1631!
1632 IF (his(ng)%IOtype.eq.io_nf90) THEN
1633 his(ng)%ncid=-1
1634#if defined PIO_LIB && defined DISTRIBUTE
1635 ELSE IF (his(ng)%IOtype.eq.io_pio) THEN
1636 his(ng)%pioFile%fh=-1
1637#endif
1638 END IF
1639 END DO
1640
1641 END DO outer_loop2
1642!
1643 10 FORMAT (a,'_',i3.3,'.nc')
1644 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
1645 & ' (Grid: ',i2.2,' TimeSteps: ',i0,' - ',i0,')',/)
1646 30 FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ', &
1647 & 1p,e17.10,0p,t68,a)
1648 40 FORMAT (/,' Converting Convolved Adjoint Trajectory to', &
1649 & ' Impulses: Outer = ',i0,' Inner = ',i0,/)
1650 50 FORMAT (/,'ROMS: ',a,1x,a,', Outer = ',i0,' Inner = ',i0,/)
1651
1652 RETURN
1653 END SUBROUTINE roms_run
1654!
1655 SUBROUTINE roms_finalize
1656!
1657!=======================================================================
1658! !
1659! This routine terminates ROMS nonlinear, tangent linear, and !
1660! adjoint models execution. !
1661! !
1662!=======================================================================
1663!
1664! Local variable declarations.
1665!
1666 integer :: Fcount, ng, thread
1667!
1668 character (len=*), parameter :: MyFile = &
1669 & __FILE__//", ROMS_finalize"
1670!
1671!-----------------------------------------------------------------------
1672! If blowing-up, save latest model state into RESTART NetCDF file.
1673!-----------------------------------------------------------------------
1674!
1675! If cycling restart records, write solution into record 3.
1676!
1677 IF (exit_flag.eq.1) THEN
1678 DO ng=1,ngrids
1679 IF (lwrtrst(ng)) THEN
1680 IF (master) WRITE (stdout,10)
1681 10 FORMAT (/,' Blowing-up: Saving latest model state into ', &
1682 & ' RESTART file',/)
1683 fcount=rst(ng)%load
1684 IF (lcyclerst(ng).and.(rst(ng)%Nrec(fcount).ge.2)) THEN
1685 rst(ng)%Rindex=2
1686 lcyclerst(ng)=.false.
1687 END IF
1690#ifdef DISTRIBUTE
1691 CALL wrt_rst (ng, myrank)
1692#else
1693 CALL wrt_rst (ng, -1)
1694#endif
1695 END IF
1696 END DO
1697 END IF
1698!
1699!-----------------------------------------------------------------------
1700! Stop model and time profiling clocks, report memory requirements, and
1701! close output NetCDF files.
1702!-----------------------------------------------------------------------
1703!
1704! Stop time clocks.
1705!
1706 IF (master) THEN
1707 WRITE (stdout,20)
1708 20 FORMAT (/,'Elapsed wall CPU time for each process (seconds):',/)
1709 END IF
1710!
1711 DO ng=1,ngrids
1712 DO thread=thread_range
1713 CALL wclock_off (ng, inlm, 0, __line__, myfile)
1714 END DO
1715 END DO
1716!
1717! Report dynamic memory and automatic memory requirements.
1718!
1719 CALL memory
1720!
1721! Close IO files.
1722!
1723 DO ng=1,ngrids
1724 CALL close_inp (ng, inlm)
1725 END DO
1726 CALL close_out
1727!
1728 RETURN
1729 END SUBROUTINE roms_finalize
1730
1731 END MODULE roms_kernel_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 memory
Definition memory.F:3
subroutine, public ad_wrt_his(ng, tile)
Definition ad_wrt_his.F:66
subroutine, public close_out
Definition close_io.F:175
subroutine, public close_file(ng, model, s, ncname, lupdate)
Definition close_io.F:43
subroutine, public close_inp(ng, model)
Definition close_io.F:92
subroutine, public congrad(ng, model, outloop, innloop, ninnloop, lcgini)
Definition congrad.F:163
subroutine, public error_covariance(model, driver, outloop, innloop, rbck, rini, rold, rnew, rec1, rec2, lposterior)
Definition convolve.F:215
subroutine, public def_impulse(ng)
Definition def_impulse.F:49
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
subroutine, public inp_par(model)
Definition inp_par.F:56
subroutine, public roms_initialize_arrays
Definition mod_arrays.F:351
subroutine, public roms_allocate_arrays(allocate_vars)
Definition mod_arrays.F:114
subroutine, public initialize_boundary(ng, tile, model)
subroutine, public initialize_forces(ng, tile, model)
type(t_fourdvar), dimension(:), allocatable fourdvar
real(r8), dimension(:), allocatable tl_obsval
logical, dimension(:), allocatable wrtmisfit
logical, dimension(:), allocatable wrtrpmod
integer, dimension(:), allocatable nobs
integer, dimension(:), allocatable nobsvar
logical, dimension(:), allocatable wrttlmod
logical, dimension(:), allocatable wrtnlmod
logical, dimension(:), allocatable wrtforce
logical lhotstart
logical, dimension(:), allocatable wrtobsscale
character(len=40), dimension(:), allocatable obsname
logical, dimension(:), allocatable wrtzetaref
integer, dimension(:), allocatable nstrobs
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 obs
type(t_io), dimension(:), allocatable tlf
type(t_io), dimension(:), allocatable tlm
type(t_io), dimension(:), allocatable itl
type(t_io), dimension(:), allocatable dav
type(t_io), dimension(:), allocatable fwd
type(t_io), dimension(:), allocatable rst
type(t_io), dimension(:), allocatable ini
integer stdout
character(len=256) sourcefile
integer, parameter io_nf90
Definition mod_ncparam.F:95
integer, parameter io_pio
Definition mod_ncparam.F:96
integer isfsur
integer, dimension(:), allocatable idefhis
subroutine, public initialize_ocean(ng, tile, model)
Definition mod_ocean.F:1526
subroutine, public initialize_parallel
integer numthreads
integer, dimension(:), allocatable first_tile
integer mythread
logical master
integer, dimension(:), allocatable last_tile
integer ocn_comm_world
integer, parameter inlm
Definition mod_param.F:662
integer, parameter irpm
Definition mod_param.F:664
integer, dimension(:), allocatable ntilex
Definition mod_param.F:685
integer, parameter iadm
Definition mod_param.F:665
integer ngrids
Definition mod_param.F:113
integer, dimension(:), allocatable ntilee
Definition mod_param.F:686
integer, parameter itlm
Definition mod_param.F:663
integer nsa
Definition mod_param.F:612
integer ninner
logical, dimension(:,:), allocatable lwrtnrm
logical, dimension(:), allocatable lreadqck
integer nouter
logical, dimension(:), allocatable ldefitl
integer blowup
logical, dimension(:), allocatable balance
logical, dimension(:), allocatable ldeftlf
integer erend
integer, dimension(:), allocatable frcrec
logical, dimension(:), allocatable lreadfrc
logical, dimension(:), allocatable ldefadj
logical, dimension(:), allocatable frequentimpulse
logical, dimension(:), allocatable ldefhis
integer, dimension(:), allocatable ntend
logical, dimension(:), allocatable ldefmod
logical, dimension(:,:), allocatable ldefnrm
integer exit_flag
logical, dimension(:), allocatable sporadicimpulse
integer erstr
logical, dimension(:), allocatable lwrthis
logical, dimension(:), allocatable lwrttlm
logical, dimension(:), allocatable lwrtrst
logical, dimension(:), allocatable ldeftlm
integer, dimension(:), allocatable ntstart
integer nrun
logical, dimension(:), allocatable lreadfwd
integer inner
integer noerror
logical, dimension(:), allocatable lcyclerst
integer outer
logical, dimension(:), allocatable lreadblk
integer, dimension(:), allocatable lold
integer, dimension(:), allocatable lbout
integer, dimension(:), allocatable lfinp
integer, dimension(:), allocatable lbinp
integer, dimension(:), allocatable lnew
integer, dimension(:), allocatable lfout
subroutine, public normalization(ng, tile, ifac)
subroutine, public roms_finalize
Definition ad_roms.h:283
subroutine, public roms_run(runinterval)
Definition ad_roms.h:239
subroutine, public roms_initialize(first, mpicomm)
Definition ad_roms.h:52
integer function, public stdout_unit(mymaster)
Definition stdout_mod.F:48
logical, save set_stdoutunit
Definition stdout_mod.F:41
character(len(sinp)) function, public uppercase(sinp)
Definition strings.F:582
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
subroutine, public tl_def_ini(ng)
Definition tl_def_ini.F:43
subroutine, public wrt_impulse(ng, tile, model, inpncname)
Definition wrt_impulse.F:59
subroutine, public wrt_ini(ng, tile, tindex, outrec)
Definition wrt_ini.F:89
subroutine, public wrt_rst(ng, tile)
Definition wrt_rst.F:63
subroutine, public biconj(ng, tile, model, lbck)
subroutine, public balance_ref(ng, tile, lbck)
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_congrad
Definition tl_congrad.F:943
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