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