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