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