ROMS
Loading...
Searching...
No Matches
r4dvar_mod Module Reference

Functions/Subroutines

subroutine, public background (my_outer, runinterval)
 
subroutine, public increment (my_outer, runinterval)
 
subroutine, public analysis (my_outer, runinterval)
 
subroutine, public prior_error (ng)
 
subroutine, public posterior_error (runinterval)
 

Variables

logical ldone
 
integer lini = 1
 
integer lbck = 2
 
integer rec1 = 1
 
integer rec2 = 2
 

Function/Subroutine Documentation

◆ analysis()

subroutine, public r4dvar_mod::analysis ( integer, intent(in) my_outer,
real(dp), intent(in) runinterval )

Definition at line 1383 of file r4dvar.F.

1384!
1385!=======================================================================
1386! !
1387! This routine computes 4D-Var data assimilation analysis, Xa. The !
1388! nonlinear model initial conditions are computed by adding the !
1389! 4D-Var increments to the current background: Xa = Xb + dXa. !
1390! !
1391! On Input: !
1392! !
1393! my_outer Outer-loop counter (integer) !
1394! RunInterval NLM kernel time stepping window (seconds) !
1395! !
1396!=======================================================================
1397!
1398! Imported variable declarations
1399!
1400 integer, intent(in) :: my_outer
1401!
1402 real(dp), intent(in) :: RunInterval
1403!
1404! Local variable declarations.
1405!
1406 integer :: i, ifile, lstr, ng
1407 integer :: InpRec
1408# ifdef PROFILE
1409 integer :: thread
1410# endif
1411# ifdef SPLIT_4DVAR
1412!
1413 real(dp) :: stime
1414# endif
1415!
1416 character (len=10) :: suffix
1417 character (len=20) :: string
1418
1419 character (len=*), parameter :: MyFile = &
1420 & __FILE__//", analysis"
1421!
1422 sourcefile=myfile
1423!
1424!-----------------------------------------------------------------------
1425! Run representer model and compute a "new estimate" of the state
1426! trajectory, X_n(t).
1427!-----------------------------------------------------------------------
1428
1429# ifdef PROFILE
1430!
1431! Start profile clock.
1432!
1433 DO ng=1,ngrids
1434 DO thread=thread_range
1435 CALL wclock_on (ng, inlm, 88, __line__, myfile)
1436 END DO
1437 END DO
1438# endif
1439
1440# ifdef SPLIT_4DVAR
1441!
1442!-----------------------------------------------------------------------
1443! If split 4D-Var algorithm, set several variables that computed or
1444! assigned in other 4D-Var phase executable.
1445!-----------------------------------------------------------------------
1446!
1447! Set Nrun>1, to read in surface forcing and open boundary conditions
1448! increments in "initial", if appropriate.
1449!
1450 nrun=ninner+1
1451!
1452! Set ERstr=Nrun, to set the open boundary condition (OBC_time) and
1453! surface forcing (SF_time) adjustment times, if needed.
1454!
1455 erstr=nrun
1456!
1457! Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
1458!
1459 DO ng=1,ngrids
1460 ldefmod(ng)=.false.
1461 CALL def_mod (ng)
1462 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1463 END DO
1464!
1465! In the split 4D-Var algorithm, we need to read the values of ObsScale
1466! computed in the background phase from the DAV NetCDF file. Such data
1467! is in memory in the unsplit algorithm.
1468!
1469! HGA: What to do in 4D-Var with nested grids?
1470!
1471 ng=1
1472 SELECT CASE (dav(ng)%IOtype)
1473 CASE (io_nf90)
1474 CALL netcdf_get_fvar (ng, irpm, dav(ng)%name, &
1475 & vname(1,idobss), obsscale, &
1476 & ncid = dav(ng)%ncid)
1477
1478# if defined PIO_LIB && defined DISTRIBUTE
1479 CASE (io_pio)
1480 CALL pio_netcdf_get_fvar (ng, irpm, dav(ng)%name, &
1481 & vname(1,idobss), obsscale, &
1482 & piofile = dav(ng)%pioFile)
1483# endif
1484 END SELECT
1485 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1486!
1487! Set nonlinear output history file to be used as the basic state
1488! trajectory. The nonlinear model is run outside of the outer loops,
1489! so we need outer=0 when updating HIS structure.
1490!
1491 DO ng=1,ngrids
1492 WRITE (his(ng)%name,10) trim(fwd(ng)%head), 0
1493 lstr=len_trim(his(ng)%name)
1494 his(ng)%base=his(ng)%name(1:lstr-3)
1495 IF (his(ng)%Nfiles.gt.1) THEN
1496 DO ifile=1,his(ng)%Nfiles
1497 WRITE (suffix,"('_',i4.4,'.nc')") ifile
1498 his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
1499 END DO
1500 his(ng)%name=trim(his(ng)%files(1))
1501 ELSE
1502 his(ng)%files(1)=trim(his(ng)%name)
1503 END IF
1504 END DO
1505!
1506! Read in 4D-Var starting time (sec) from forward trajectory NetCDF
1507! file and set INItime. The starting time is need for boundary and
1508! surface forcing adjustments, if any.
1509!
1510 inprec=1
1511 DO ng=1,ngrids
1512 SELECT CASE (his(ng)%IOtype)
1513 CASE (io_nf90)
1514 CALL netcdf_get_fvar (ng, irpm, his(ng)%name, &
1515 & vname(1,idtime), stime, &
1516 & start = (/inprec/), &
1517 & total = (/1/))
1518
1519# if defined PIO_LIB && defined DISTRIBUTE
1520 CASE (io_pio)
1521 CALL pio_netcdf_get_fvar (ng, irpm, his(ng)%name, &
1522 & vname(1,idtime), stime, &
1523 & start = (/inprec/), &
1524 & total = (/1/))
1525# endif
1526 END SELECT
1527 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1528 initime(ng)=stime
1529
1530# ifdef ADJUST_BOUNDARY
1531!
1532! Set time (sec) for the open boundary adjustment.
1533!
1534 obc_time(1,ng)=stime
1535 DO i=2,nbrec(ng)
1536 obc_time(i,ng)=obc_time(i-1,ng)+nobc(ng)*dt(ng)
1537 END DO
1538# endif
1539
1540# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
1541!
1542! Set time (sec) for the surface forcing adjustment.
1543!
1544 sf_time(1,ng)=stime
1545 DO i=2,nfrec(ng)
1546 sf_time(i,ng)=sf_time(i-1,ng)+nsff(ng)*dt(ng)
1547 END DO
1548# endif
1549 END DO
1550!
1551! If outer>1, set previous outer loop representer model file.
1552!
1553 IF (my_outer.gt.1) THEN
1554 DO ng=1,ngrids
1555 WRITE (tlm(ng)%name,10) trim(fwd(ng)%head), my_outer-1
1556 lstr=len_trim(tlm(ng)%name)
1557 tlm(ng)%base=tlm(ng)%name(1:lstr-3)
1558 END DO
1559 END IF
1560!
1561! Set basic state trajectory to either Nonlinear model (outer=1) or
1562! previous outer loop representer model (outer>1). Also, set switches
1563! to process the FWD structure in routine "check_multifile". Notice
1564! that it is possible to split solution into multiple NetCDF files to
1565! reduce their size.
1566!
1567 IF (my_outer.eq.1) THEN
1568 CALL edit_multifile ('HIS2FWD')
1569 ELSE
1570 CALL edit_multifile ('TLM2FWD')
1571 END IF
1572 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1573 DO ng=1,ngrids
1574 lreadfwd(ng)=.true.
1575 END DO
1576
1577# ifdef FORWARD_FLUXES
1578!
1579! Set the BLK structure to contain the nonlinear model surface fluxes
1580! needed by the tangent linear and adjoint models. Also, set switches
1581! to process that structure in routine "check_multifile". Notice that
1582! it is possible to split the solution into multiple NetCDF files to
1583! reduce their size.
1584!
1585! The switch LreadFRC is deactivated because all the atmospheric
1586! forcing, including shortwave radiation, is read from the NLM
1587! surface fluxes or is assigned during ESM coupling. Such fluxes
1588! are available from the QCK structure. There is no need for reading
1589! and processing from the FRC structure input forcing-files.
1590!
1591 CALL edit_multifile ('QCK2BLK')
1592 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1593 DO ng=1,ngrids
1594 lreadblk(ng)=.true.
1595 lreadfrc(ng)=.false.
1596 END DO
1597# endif
1598# endif
1599!
1600!-----------------------------------------------------------------------
1601! Run representer model with the new analysis initial conditions
1602! cpmputed in the "increment" phase.
1603!-----------------------------------------------------------------------
1604!
1605! Set new basic state trajectory for next outer loop.
1606!
1607 DO ng=1,ngrids
1608 ideftlm(ng)=-1
1609 ldeftlm(ng)=.true.
1610 lwrttlm(ng)=.true.
1611 wrtnlmod(ng)=.false.
1612 wrttlmod(ng)=.true.
1613 wrtrpmod(ng)=.true.
1614 WRITE (tlm(ng)%name,10) trim(fwd(ng)%head), my_outer
1615 lstr=len_trim(tlm(ng)%name)
1616 tlm(ng)%base=tlm(ng)%name(1:lstr-3)
1617 END DO
1618!
1619! If weak constraint, the impulses are time-interpolated at each
1620! time-steps.
1621!
1622 DO ng=1,ngrids
1623 IF (frcrec(ng).gt.3) THEN
1624 frequentimpulse(ng)=.true.
1625 END IF
1626 END DO
1627!
1628! Initialize representer model IRP(ng)%name file, record Rec2.
1629!
1630 DO ng=1,ngrids
1631# ifdef DATALESS_LOOPS
1632 irp(ng)%Rindex=rec1
1633# else
1634 irp(ng)%Rindex=rec2
1635# endif
1636 CALL rp_initial (ng)
1637 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1638 END DO
1639!
1640! Activate switch to write out final misfit between model and
1641! observations.
1642!
1643 IF (my_outer.eq.nouter) THEN
1644 DO ng=1,ngrids
1645 wrtmisfit(ng)=.true.
1646 END DO
1647 END IF
1648!
1649! Run representer model using previous linearized trajectory, X_n-1, as
1650! basic state and forced with convolved adjoint trajectory impulses.
1651!
1652 DO ng=1,ngrids
1653# ifdef RP_AVERAGES
1654 idefavg(ng)=-1
1655 ldefavg(ng)=.true.
1656 lwrtavg(ng)=.true.
1657 WRITE (avg(ng)%name,10) trim(avg(ng)%head), my_outer
1658 lstr=len_trim(avg(ng)%name)
1659 avg(ng)%base=avg(ng)%name(1:lstr-3)
1660# endif
1661 IF (master) THEN
1662 WRITE (stdout,20) 'RP', ng, my_outer, ninner, &
1663 & ntstart(ng), ntend(ng)
1664 END IF
1665 END DO
1666!
1667# ifdef SOLVE3D
1668 CALL rp_main3d (runinterval)
1669# else
1670 CALL rp_main2d (runinterval)
1671# endif
1672 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1673!
1674# ifdef RP_AVERAGES
1675 DO ng=1,ngrids
1676 ldefavg(ng)=.false.
1677 lwrtavg(ng)=.false.
1678 END DO
1679# endif
1680!
1681! Report data penalty function.
1682!
1683 DO ng=1,ngrids
1684 IF (master) THEN
1685 DO i=0,nstatevar(ng)
1686 IF (i.eq.0) THEN
1687 string='Total'
1688 ELSE
1689 string=vname(1,idsvar(i))
1690 END IF
1691 IF (fourdvar(ng)%DataPenalty(i).ne.0.0_r8) THEN
1692 WRITE (stdout,30) my_outer, ninner, 'RPM', &
1693 & fourdvar(ng)%DataPenalty(i), &
1694 & trim(string)
1695# ifdef DATALESS_LOOPS
1696 WRITE (stdout,30) my_outer, ninner, 'NLM', &
1697 & fourdvar(ng)%NLPenalty(i), &
1698 & trim(string)
1699# endif
1700 END IF
1701 END DO
1702 END IF
1703 END DO
1704!
1705! Write out final data penalty function to NetCDF file.
1706!
1707 sourcefile=myfile
1708 DO ng=1,ngrids
1709 SELECT CASE (dav(ng)%IOtype)
1710 CASE (io_nf90)
1711 CALL netcdf_put_fvar (ng, irpm, dav(ng)%name, &
1712 & 'RP_fDataPenalty', &
1713 & fourdvar(ng)%DataPenalty(0:), &
1714 & (/1,my_outer/), &
1715 & (/nstatevar(ng)+1,1/), &
1716 & ncid = dav(ng)%ncid)
1717
1718# if defined PIO_LIB && defined DISTRIBUTE
1719 CASE (io_pio)
1720 CALL pio_netcdf_put_fvar (ng, irpm, dav(ng)%name, &
1721 & 'RP_fDataPenalty', &
1722 & fourdvar(ng)%DataPenalty(0:), &
1723 & (/1,my_outer/), &
1724 & (/nstatevar(ng)+1,1/), &
1725 & piofile = dav(ng)%pioFile)
1726# endif
1727 END SELECT
1728 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1729 END DO
1730!
1731! Clean array before next run of RP model.
1732!
1733 DO ng=1,ngrids
1734 fourdvar(ng)%DataPenalty=0.0_r8
1735# ifdef DATALESS_LOOPS
1736 fourdvar(ng)%NLPenalty=0.0_r8
1737# endif
1738 wrtnlmod(ng)=.false.
1739 wrttlmod(ng)=.false.
1740 END DO
1741!
1742! Close current forward NetCDF file. If last outer loop, set history
1743! file ID to closed state since we manipulated its indices with the
1744! forward file ID, which it was closed.
1745!
1746 DO ng=1,ngrids
1747 CALL close_file (ng, irpm, fwd(ng))
1748 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1749 IF (my_outer.eq.nouter) THEN
1750 IF (his(ng)%IOtype.eq.io_nf90) THEN
1751 his(ng)%ncid=-1
1752# if defined PIO_LIB && defined DISTRIBUTE
1753 ELSE IF (his(ng)%IOtype.eq.io_pio) THEN
1754 his(ng)%pioFile%fh=-1
1755# endif
1756 END IF
1757 END IF
1758 END DO
1759
1760# ifdef PROFILE
1761!
1762! Stop profile clock
1763!
1764 DO ng=1,ngrids
1765 DO thread=thread_range
1766 CALL wclock_off (ng, inlm, 88, __line__, myfile)
1767 END DO
1768 END DO
1769# endif
1770!
1771 10 FORMAT (a,'_outer',i0,'.nc')
1772 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
1773 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
1774 ', TimeSteps: ',i0,' - ',i0,')',/)
1775 30 FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ', &
1776 & 1p,e17.10,0p,t68,a)
1777!
1778 RETURN
subroutine edit_multifile(task)
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

References mod_iounits::avg, close_io_mod::close_file(), mod_iounits::dav, def_mod_mod::def_mod(), mod_scalars::dt, edit_multifile(), mod_scalars::erstr, mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_scalars::frcrec, mod_scalars::frequentimpulse, mod_iounits::fwd, mod_iounits::his, mod_ncparam::idefavg, mod_ncparam::ideftlm, mod_ncparam::idobss, mod_ncparam::idsvar, mod_ncparam::idtime, mod_scalars::initime, mod_param::inlm, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_iounits::irp, mod_param::irpm, mod_scalars::ldefavg, mod_scalars::ldefmod, mod_scalars::ldeftlm, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadfwd, mod_scalars::lwrtavg, mod_scalars::lwrttlm, mod_parallel::master, mod_scalars::nbrec, mod_scalars::nfrec, mod_param::ngrids, mod_scalars::ninner, mod_scalars::nobc, mod_scalars::noerror, mod_scalars::nouter, mod_scalars::nrun, mod_scalars::nsff, mod_fourdvar::nstatevar, mod_scalars::ntend, mod_scalars::ntstart, mod_scalars::obc_time, mod_fourdvar::obsscale, rec1, rec2, rp_initial(), rp_main2d(), rp_main3d(), mod_scalars::sf_time, mod_iounits::sourcefile, mod_iounits::stdout, mod_iounits::tlm, mod_ncparam::vname, wclock_off(), wclock_on(), mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtrpmod, and mod_fourdvar::wrttlmod.

Here is the call graph for this function:

◆ background()

subroutine, public r4dvar_mod::background ( integer, intent(in) my_outer,
real(dp), intent(in) runinterval )

Definition at line 157 of file r4dvar.F.

158!
159!=======================================================================
160! !
161! This routine computes the backgound state trajectory, Xb_n-1(t), !
162! used to linearize the tangent linear and adjoint models in the !
163! inner loops. It interpolates the background at the observations !
164! locations, and computes the accept/reject quality control flag, !
165! ObsScale. !
166! !
167! On Input: !
168! !
169! my_outer Outer-loop counter (integer) !
170! RunInterval NLM kernel time stepping window (seconds) !
171! !
172!=======================================================================
173!
174! Imported variable declarations
175!
176 integer, intent(in) :: my_outer
177!
178 real(dp), intent(in) :: RunInterval
179!
180! Local variable declarations.
181!
182 integer :: lstr, ng
183 integer :: Fcount, Tindex
184# ifdef PROFILE
185 integer :: thread
186# endif
187# if defined MODEL_COUPLING && !defined MCT_LIB
188 integer :: NstrStep, NendStep, extra
189!
190 real(dp) :: ENDtime, NEXTtime
191# endif
192!
193 character (len=*), parameter :: MyFile = &
194 & __FILE__//", background"
195!
196 sourcefile=myfile
197!
198!-----------------------------------------------------------------------
199! Initialize and set nonlinear model initial conditions.
200!-----------------------------------------------------------------------
201
202# ifdef PROFILE
203!
204! Start profile clock.
205!
206 DO ng=1,ngrids
207 DO thread=thread_range
208 CALL wclock_on (ng, inlm, 86, __line__, myfile)
209 END DO
210 END DO
211# endif
212!
213! Initialize the switch to gather weak constraint forcing.
214!
215 DO ng=1,ngrids
216 wrtforce(ng)=.false.
217 END DO
218!
219! Initialize and set nonlinear model initial conditions.
220!
221 DO ng=1,ngrids
222 wrtnlmod(ng)=.true.
223 wrtrpmod(ng)=.false.
224 wrttlmod(ng)=.false.
225 END DO
226!
227 CALL initial
228 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
229!
230! Save nonlinear initial conditions (currently in time index 1,
231! background) into record "Lini" of INI(ng)%name NetCDF file. The
232! record "Lbck" becomes the background state record and the record
233! "Lini" becomes current nonlinear initial conditions.
234!
235 tindex=1
236 DO ng=1,ngrids
237 ini(ng)%Rindex=1
238 fcount=ini(ng)%load
239 ini(ng)%Nrec(fcount)=1
240# ifdef DISTRIBUTE
241 CALL wrt_ini (ng, myrank, tindex)
242# else
243 CALL wrt_ini (ng, -1, tindex)
244# endif
245 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
246 END DO
247
248# ifdef STD_MODEL
249!
250! Estimate initial conditions standard deviation fields for the
251! background error covariance using the method of Mogensen et al.
252! (2012), which assumes that the background errors are proportional
253! to the vertical derivatives of the state. The error have simmilar
254! shape as the prior profile but iwth a displacement over the water
255! column.
256!
257 IF (my_outer.eq.0) THEN
258 lstd=1 ! time index in error arrays
259 DO ng=1,ngrids
260 DO tile=first_tile(ng),last_tile(ng),+1
261 CALL background_std (ng, tile, tindex, lstd)
262 END DO
263 END DO
264!
265! Write out estimated background standard deviation fields.
266!
267 DO ng=1,ngrids
268 IF (lwrtstd(ng)) THEN
269 std(5,ng)%Rindex=0
270 fcount=std(5,ng)%load
271 std(5,ng)%Nrec(fcount)=1
272# ifdef DISTRIBUTE
273 CALL wrt_std (ng, myrank, lstd)
274# else
275 CALL wrt_std (ng, -1, lstd)
276# endif
277 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
278 END IF
279 END DO
280 END IF
281# endif
282!
283! Set nonlinear output history file as the initial basic state
284! trajectory.
285!
286 DO ng=1,ngrids
287 ldefhis(ng)=.true.
288 lwrthis(ng)=.true.
289 WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer
290 lstr=len_trim(his(ng)%name)
291 his(ng)%base=his(ng)%name(1:lstr-3)
292 END DO
293!
294! Define output 4DVAR NetCDF file containing all processed data
295! at observation locations.
296!
297 DO ng=1,ngrids
298 ldefmod(ng)=.true.
299 CALL def_mod (ng)
300 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
301 END DO
302!
303!-----------------------------------------------------------------------
304! Run nonlinear model and compute basic state trajectory. It processes
305! and writes the observations accept/reject flag (ObsScale) once to
306! allow background quality control, if any.
307# if defined MODEL_COUPLING && !defined MCT_LIB
308! Since the ROMS kernel has a delayed output and line diagnostics by
309! one timestep, subtact an extra value to the report of starting and
310! ending timestep for clarity. Usually, the model coupling interval
311! is of the same size as ROMS timestep.
312# endif
313!-----------------------------------------------------------------------
314!
315 myruninterval=runinterval
316 IF (master) WRITE (stdout,'(1x)')
317!
318 DO ng=1,ngrids
319# ifdef RP_AVERAGES
320 idefavg(ng)=-1
321 ldefavg(ng)=.false.
322 lwrtavg(ng)=.false.
323 WRITE (avg(ng)%name,10) trim(avg(ng)%head), my_outer
324 lstr=len_trim(avg(ng)%name)
325 avg(ng)%base=avg(ng)%name(1:lstr-3)
326# endif
327 wrtobsscale(ng)=.true.
328# if defined MODEL_COUPLING && !defined MCT_LIB
329!
330 nexttime=time(ng)+runinterval
331 endtime=initime(ng)+(ntimes(ng)-1)*dt(ng)
332 IF ((nexttime.eq.endtime).and.(ng.eq.1)) THEN
333 extra=0 ! last time interval
334 ELSE
335 extra=1
336 END IF
337 step_counter(ng)=0
338 nstrstep=iic(ng)
339 nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
340 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0, &
341 & nstrstep, nendstep
342# else
343 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0, &
344 & ntstart(ng), ntend(ng)
345# endif
346 END DO
347!
348# ifdef SOLVE3D
349 CALL main3d (runinterval)
350# else
351 CALL main2d (runinterval)
352# endif
353 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
354!
355 DO ng=1,ngrids
356 wrtnlmod(ng)=.false.
357 wrtobsscale(ng)=.false.
358 END DO
359
360# ifdef PROFILE
361!
362! Stop profile clock
363!
364 DO ng=1,ngrids
365 DO thread=thread_range
366 CALL wclock_off (ng, inlm, 86, __line__, myfile)
367 END DO
368 END DO
369# endif
370!
371 10 FORMAT (a,'_outer',i0,'.nc')
372 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
373 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
374 ', TimeSteps: ',i0,' - ',i0,')',/)
375!
376 RETURN
subroutine initial
Definition initial.F:3
subroutine main2d
Definition main2d.F:746
subroutine main3d(runinterval)
Definition main3d.F:4

References mod_iounits::avg, def_mod_mod::def_mod(), mod_scalars::dt, mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_iounits::fwd, mod_iounits::his, mod_ncparam::idefavg, mod_scalars::iic, mod_iounits::ini, initial(), mod_scalars::initime, mod_param::inlm, mod_parallel::last_tile, mod_scalars::ldefavg, mod_scalars::ldefhis, mod_scalars::ldefmod, mod_scalars::lwrtavg, mod_scalars::lwrthis, main2d(), main3d(), mod_parallel::master, mod_parallel::myrank, mod_scalars::myruninterval, mod_param::ngrids, mod_scalars::noerror, mod_scalars::ntend, mod_scalars::ntimes, mod_scalars::ntstart, mod_iounits::sourcefile, mod_iounits::std, mod_iounits::stdout, mod_scalars::step_counter, mod_scalars::time, wclock_off(), wclock_on(), wrt_ini_mod::wrt_ini(), mod_fourdvar::wrtforce, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtobsscale, mod_fourdvar::wrtrpmod, and mod_fourdvar::wrttlmod.

Here is the call graph for this function:

◆ increment()

subroutine, public r4dvar_mod::increment ( integer, intent(in) my_outer,
real(dp), intent(in) runinterval )

Definition at line 379 of file r4dvar.F.

380!
381!=======================================================================
382! !
383! This routine computes the 4D-Var data assimilation state increment, !
384! dXa, by iterating the inner loops and minimizing the cost function. !
385! !
386! On Input: !
387! !
388! my_outer Outer-loop counter (integer) !
389! RunInterval TLM/ADM kernels time stepping window (seconds) !
390! !
391! !
392! It solves the system: !
393! !
394! (R_n + Cobs) * Beta_n = h_n !
395! !
396! h_n = Xo - H * X_n !
397! !
398! where R_n is the representer matrix, Cobs is the observation-error !
399! covariance, Beta_n are the representer coefficients, h_n is the !
400! misfit between observations (Xo) and model (H * X_n), and H is !
401! the linearized observation operator. Here, _n denotes a sequence !
402! of estimates. !
403! !
404! The system does not need to be solved explicitly by inverting the !
405! symmetric stabilized representer matrix, P_n: !
406! !
407! P_n = R_n + Cobs !
408! !
409! but by computing the action of P_n on any vector PSI, such that !
410! !
411! P_n * PSI = R_n * PSI + Cobs * PSI !
412! !
413! The representer matrix is not explicitly computed but evaluated by !
414! one integration backward of the adjoint model and one integration !
415! forward of the tangent linear model for any forcing vector PSI. !
416! !
417! A preconditioned conjugate gradient algorithm is used to compute !
418! an approximation PSI for Beta_n. !
419! !
420!=======================================================================
421!
422! Imported variable declarations
423!
424 integer, intent(in) :: my_outer
425!
426 real(dp), intent(in) :: RunInterval
427!
428! Local variable declarations.
429!
430 logical :: Lcgini, Linner, Lposterior, add
431!
432 integer :: i, ifile, lstr, my_inner, ng, tile
433 integer :: Fcount, InpRec, Tindex
434# ifdef PROFILE
435 integer :: thread
436# endif
437# ifdef SPLIT_4DVAR
438!
439 real(dp) :: stime
440# endif
441!
442 character (len=6 ) :: driver = 'r4dvar'
443 character (len=10 ) :: suffix
444 character (len=20 ) :: string
445# ifdef SPLIT_4DVAR
446 character (len=256) :: ncname
447# endif
448 character (len=*), parameter :: MyFile = &
449 & __FILE__//", increment"
450!
451 sourcefile=myfile
452!
453!-----------------------------------------------------------------------
454! Run representer model and compute a "prior estimate" state
455! trajectory, X_n(t). Use linearized state trajectory (X_n-1) as
456! basic state.
457!-----------------------------------------------------------------------
458
459# ifdef PROFILE
460!
461! Start profile clock.
462!
463 DO ng=1,ngrids
464 DO thread=thread_range
465 CALL wclock_on (ng, itlm, 87, __line__, myfile)
466 END DO
467 END DO
468# endif
469
470# ifdef SPLIT_4DVAR
471!
472!-----------------------------------------------------------------------
473! If split 4D-Var algorithm, set several variables that are computed
474! or assigned in other 4D-Var phase executable.
475!-----------------------------------------------------------------------
476!
477! Reset Nrun counter to a value greater than one.
478!
479 IF (my_outer.gt.1) THEN
480 nrun=1+(my_outer-1)*ninner
481 END IF
482!
483! Open Nonlinear model initial conditions NetCDF file (INI struc) and
484! inquire about its variables IDs.
485!
486 DO ng=1,ngrids
487 ldefini(ng)=.false.
488 CALL def_ini (ng)
489 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
490 END DO
491!
492! If outer>1, open tangent linear initial conditions NetCDF file
493! (ITL struc) and inquire about its variables IDs.
494!
495 IF (my_outer.gt.1) THEN
496 DO ng=1,ngrids
497 ldefitl(ng)=.false.
498 CALL tl_def_ini (ng)
499 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
500 END DO
501 END IF
502!
503! If outer>1, open representer model initial conditions NetCDF file
504! (IRP struc) and inquire about its variables IDs.
505!
506 IF (my_outer.gt.1) THEN
507 DO ng=1,ngrids
508 ldefirp(ng)=.false.
509 CALL rp_def_ini (ng)
510 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
511 END DO
512 END IF
513!
514! If outer>1, open impulse forcing NetCDF file (TLF struc) and inquire
515! about its variables IDs.
516!
517 IF (my_outer.gt.1) THEN
518 DO ng=1,ngrids
519 ldeftlf(ng)=.false.
520 CALL def_impulse (ng)
521 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
522 END DO
523 END IF
524!
525! If outer>1, open adjoint history NetCDF file (ADM struc) and inquire
526! about its variables IDs and set "FrcRec".
527!
528 IF (my_outer.gt.1) THEN
529 DO ng=1,ngrids
530 ldefadj(ng)=.false.
531 CALL ad_def_his (ng, ldefadj(ng))
532 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
533 fcount=adm(ng)%Fcount
534 frcrec(ng)=adm(ng)%Nrec(fcount)
535 END DO
536 END IF
537!
538! Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
539! Deactivate "haveNLmod" since its values are read below. Otherwise,
540! the TLM will read zero values when calling "obs_read" for outer>1.
541! The DAV file does not have values for NLmodel_value(:,outer) yet.
542!
543 DO ng=1,ngrids
544 ldefmod(ng)=.false.
545 CALL def_mod (ng)
546 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
547 havenlmod(ng)=.false. ! because is activated in "def_mod"
548 END DO
549!
550! If split 4D-Var and outer>1, read several variables from 4D-VAR
551! NetCDF (DAV struc) file needed in the conjugate gradient algorithm.
552! In the unsplit case, such values are available in memory.
553!
554 IF (my_outer.gt.1) THEN
555 DO ng=1,ngrids
556 CALL cg_read_congrad (ng, itlm, my_outer)
557 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
558 END DO
559 END IF
560!
561! In the split 4D-Var algorithm, we need to read in initial NLmodVal
562! ("NLmodel_initial") and ObsScale ("obs_scale"), computed in the
563! background phase, from the DAV NetCDF file. Such data is in memory
564! in the unsplit algorithm.
565!
566! HGA: What to do in 4D-Var with nested grids?
567!
568 IF (my_outer.eq.1) THEN
569 ng=1 ! initial values from "background"
570 SELECT CASE (dav(ng)%IOtype)
571 CASE (io_nf90)
572 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
573 & vname(1,idnlmi), nlmodval, &
574 & ncid = dav(ng)%ncid)
575
576# if defined PIO_LIB && defined DISTRIBUTE
577 CASE (io_pio)
578 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
579 & vname(1,idnlmi), nlmodval, &
580 & piofile = dav(ng)%pioFile)
581# endif
582 END SELECT
583 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
584 END IF
585!
586 ng=1
587 SELECT CASE (dav(ng)%IOtype)
588 CASE (io_nf90)
589 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
590 & vname(1,idobss), obsscale, &
591 & ncid = dav(ng)%ncid)
592
593# if defined PIO_LIB && defined DISTRIBUTE
594 CASE (io_pio)
595 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
596 & vname(1,idobss), obsscale, &
597 & piofile = dav(ng)%pioFile)
598# endif
599 END SELECT
600 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
601!
602! In the split 4D-Var alorithm, we also need to read in ObsVal
603! ("obs_value") and ObsErr ("obs_error") from the observation file.
604! They are used in the initialization of the conjugate gradient
605! algorithm.
606!
607 ng=1
608 SELECT CASE (obs(ng)%IOtype)
609 CASE (io_nf90)
610 CALL netcdf_get_fvar (ng, itlm, obs(ng)%name, &
611 & vname(1,idoval), obsval)
612 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
613!
614 CALL netcdf_get_fvar (ng, itlm, obs(ng)%name, &
615 & vname(1,idoerr), obserr)
616 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
617
618# if defined PIO_LIB && defined DISTRIBUTE
619 CASE (io_pio)
620 CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name, &
621 & vname(1,idoval), obsval)
622 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
623!
624 CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name, &
625 & vname(1,idoerr), obserr)
626 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
627# endif
628 END SELECT
629!
630! Set nonlinear output history file to be used as the basic state
631! trajectory. The nonlinear model is run outside of the outer loops,
632! so we need outer=0 when updating HIS structure. Notice that this
633! structure is needed to set the BLK structure below.
634!
635 DO ng=1,ngrids
636 WRITE (his(ng)%name,10) trim(fwd(ng)%head), 0
637 lstr=len_trim(his(ng)%name)
638 his(ng)%base=his(ng)%name(1:lstr-3)
639 IF (his(ng)%Nfiles.gt.1) THEN
640 DO ifile=1,his(ng)%Nfiles
641 WRITE (suffix,"('_',i4.4,'.nc')") ifile
642 his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
643 END DO
644 his(ng)%name=trim(his(ng)%files(1))
645 ELSE
646 his(ng)%files(1)=trim(his(ng)%name)
647 END IF
648 END DO
649!
650! Read in 4D-Var starting time (sec) from forward trajectory NetCDF
651! file and set INItime. The starting time is need for boundary and
652! surface forcing adjustments, if any.
653!
654 inprec=1
655 DO ng=1,ngrids
656 SELECT CASE (his(ng)%IOtype)
657 CASE (io_nf90)
658 CALL netcdf_get_fvar (ng, itlm, his(ng)%name, &
659 & vname(1,idtime), stime, &
660 & start = (/inprec/), &
661 & total = (/1/))
662
663# if defined PIO_LIB && defined DISTRIBUTE
664 CASE (io_pio)
665 CALL pio_netcdf_get_fvar (ng, itlm, his(ng)%name, &
666 & vname(1,idtime), stime, &
667 & start = (/inprec/), &
668 & total = (/1/))
669# endif
670 END SELECT
671 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
672 initime(ng)=stime
673
674# ifdef ADJUST_BOUNDARY
675!
676! Set time (sec) for the open boundary adjustment.
677!
678 obc_time(1,ng)=stime
679 DO i=2,nbrec(ng)
680 obc_time(i,ng)=obc_time(i-1,ng)+nobc(ng)*dt(ng)
681 END DO
682# endif
683
684# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
685!
686! Set time (sec) for the surface forcing adjustment.
687!
688 sf_time(1,ng)=stime
689 DO i=2,nfrec(ng)
690 sf_time(i,ng)=sf_time(i-1,ng)+nsff(ng)*dt(ng)
691 END DO
692# endif
693 END DO
694
695# if defined ADJUST_BOUNDARY || \
696 defined adjust_stflux || defined adjust_wstress
697!
698! If outer>1, read in adjust forcing arrays from IRP file (Rec2),
699! which are computed in the previous outer loop when calling the
700! "error_covarinace" after the inner loops. Such values are in
701! memory in the unsplit algorithm. Use "iTLM" descriptor so the
702! switch get_adjust=.TRUE. Notice that all the TL arrays are
703! cleared in the next call to "rp_initial" except the "tl_*_obc"
704! arrays used when adjusting OBCs.
705!
706 IF (my_outer.gt.1) THEN
707 DO ng=1,ngrids
708 inprec=rec2
709 tindex=1
710 CALL get_state (ng, itlm, 1, irp(ng), inprec, tindex)
711 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
712 END DO
713 END IF
714# endif
715!
716! If outer>1, set previous outer loop representer model file.
717!
718 IF (my_outer.gt.1) THEN
719 DO ng=1,ngrids
720 WRITE (tlm(ng)%name,10) trim(fwd(ng)%head), my_outer-1
721 lstr=len_trim(tlm(ng)%name)
722 tlm(ng)%base=tlm(ng)%name(1:lstr-3)
723 END DO
724 END IF
725# endif
726!
727!-----------------------------------------------------------------------
728! On first pass (outer=1), create NetCDF files and initialize all the
729! records in the ITL file to zero.
730!-----------------------------------------------------------------------
731!
732 check_outer1 : IF (my_outer.eq.1) THEN
733!
734! If outer=1, define tangent linear initial conditions file.
735!
736 DO ng=1,ngrids
737 ldefitl(ng)=.true.
738 CALL tl_def_ini (ng)
739 ldefitl(ng)=.false.
740 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
741 END DO
742!
743! If outer=1, define TLM/RPM impulse forcing NetCDF file.
744!
745 DO ng=1,ngrids
746 ldeftlf(ng)=.true.
747 CALL def_impulse (ng)
748 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
749 END DO
750
751# if defined POSTERIOR_EOFS || defined POSTERIOR_ERROR_I || \
752 defined posterior_error_f
753!
754! Define output Hessian NetCDF file that will eventually contain
755! the intermediate posterior analysis error covariance matrix
756! fields or its EOFs.
757!
758 DO ng=1,ngrids
759 ldefhss(ng)=.true.
760 CALL def_hessian (ng)
761 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
762 END DO
763# endif
764
765# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
766!
767! Define output initial or final full posterior error covariance
768! (diagonal) matrix NetCDF.
769!
770 DO ng=1,ngrids
771 ldeferr(ng)=.true.
772 CALL def_error (ng)
773 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
774 END DO
775# endif
776 END IF check_outer1
777!
778!-----------------------------------------------------------------------
779! Run finite amplitude tangent linear (representers) model.
780!-----------------------------------------------------------------------
781!
782! Clear the nonlinear state arrays to make sure that the RHS terms
783! rzeta, rubar, rvbar, ru, and rv are zero. Otherwise, those arrays
784! will have the last computed values when running the nonlinear model
785! if not processing the forward trajectory RHS terms. This needs to
786! be done to get identical solutions with the split schemes.
787!
788 DO ng=1,ngrids
789 DO tile=first_tile(ng),last_tile(ng),+1
790 CALL initialize_ocean (ng, tile, inlm)
791 END DO
792 END DO
793!
794! Set structure for the nonlinear forward trajectory to be processed
795! by the tangent linear and adjoint models. Also, set switches to
796! process the FWD structure in routine "check_multifile". Notice that
797! it is possible to split solution into multiple NetCDF files to reduce
798! their size.
799!
800 IF (my_outer.eq.1) THEN
801 CALL edit_multifile ('HIS2FWD')
802 ELSE
803 CALL edit_multifile ('TLM2FWD')
804 END IF
805 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
806 DO ng=1,ngrids
807 lreadfwd(ng)=.true.
808 END DO
809
810# ifdef FORWARD_FLUXES
811!
812! Set the BLK structure to contain the nonlinear model surface fluxes
813! needed by the tangent linear and adjoint models. Also, set switches
814! to process that structure in routine "check_multifile". Notice that
815! it is possible to split the solution into multiple NetCDF files to
816! reduce their size.
817!
818! The switch LreadFRC is deactivated because all the atmospheric
819! forcing, including shortwave radiation, is read from the NLM
820! surface fluxes or is assigned during ESM coupling. Such fluxes
821! are available from the QCK structure. There is no need for reading
822! and processing from the FRC structure input forcing-files.
823!
824 CALL edit_multifile ('QCK2BLK')
825 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
826 DO ng=1,ngrids
827 lreadblk(ng)=.true.
828 lreadfrc(ng)=.false.
829 END DO
830# endif
831!
832! Set representer model output file name. The strategy is to write
833! the representer solution at the beginning of each outer loop.
834!
835 DO ng=1,ngrids
836 ideftlm(ng)=-1
837 ldeftlm(ng)=.true.
838 lwrttlm(ng)=.true.
839 WRITE (tlm(ng)%name,10) trim(fwd(ng)%head), my_outer
840 lstr=len_trim(tlm(ng)%name)
841 tlm(ng)%base=tlm(ng)%name(1:lstr-3)
842 END DO
843!
844! Activate switch to write the representer model at observation points.
845! Turn off writing into history file and turn off impulse forcing.
846!
847 DO ng=1,ngrids
848 wrtrpmod(ng)=.true.
849 sporadicimpulse(ng)=.false.
850 frequentimpulse(ng)=.false.
851 END DO
852
853# ifndef DATALESS_LOOPS
854!
855! As in the nonlinear model, initialize always the representer model
856! here with the background or reference state (IRP(ng)%name, record
857! Rec1).
858!
859 DO ng=1,ngrids
860 irp(ng)%Rindex=rec1
861 CALL rp_initial (ng)
862 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
863 END DO
864!
865! Run representer model using the nonlinear trajectory as a basic
866! state. Compute model solution at observation points, H * X_n.
867!
868 DO ng=1,ngrids
869# ifdef RP_AVERAGES
870 ldefavg(ng)=.false.
871 lwrtavg(ng)=.false.
872# endif
873 IF (master) THEN
874 WRITE (stdout,20) 'RP', ng, my_outer, 0, &
875 & ntstart(ng), ntend(ng)
876 END IF
877 END DO
878!
879# ifdef SOLVE3D
880 CALL rp_main3d (runinterval)
881# else
882 CALL rp_main2d (runinterval)
883# endif
884 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
885!
886! Report data penalty function.
887!
888 DO ng=1,ngrids
889 IF (master) THEN
890 DO i=0,nstatevar(ng)
891 IF (i.eq.0) THEN
892 string='Total'
893 ELSE
894 string=vname(1,idsvar(i))
895 END IF
896 IF (fourdvar(ng)%DataPenalty(i).ne.0.0_r8) THEN
897 WRITE (stdout,30) my_outer, inner, 'RPM', &
898 & fourdvar(ng)%DataPenalty(i), &
899 & trim(string)
900 END IF
901 END DO
902 END IF
903!
904! Write out initial data penalty function to NetCDF file.
905!
906 sourcefile=myfile
907 SELECT CASE (dav(ng)%IOtype)
908 CASE (io_nf90)
909 CALL netcdf_put_fvar (ng, irpm, dav(ng)%name, &
910 & 'RP_iDataPenalty', &
911 & fourdvar(ng)%DataPenalty(0:), &
912 & (/1,my_outer/), &
913 & (/nstatevar(ng)+1,1/), &
914 & ncid = dav(ng)%ncid)
915
916# if defined PIO_LIB && defined DISTRIBUTE
917 CASE (io_pio)
918 CALL pio_netcdf_put_fvar (ng, irpm, dav(ng)%name, &
919 & 'RP_iDataPenalty', &
920 & fourdvar(ng)%DataPenalty(0:), &
921 & (/1,my_outer/), &
922 & (/nstatevar(ng)+1,1/), &
923 & piofile = dav(ng)%pioFile)
924# endif
925 END SELECT
926 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
927!
928! Clean penalty array before next run of RP model.
929!
930 fourdvar(ng)%DataPenalty=0.0_r8
931 END DO
932!
933! Turn off IO switches.
934!
935 DO ng=1,ngrids
936 ldeftlm(ng)=.false.
937 lwrttlm(ng)=.false.
938 wrtrpmod(ng)=.false.
939 END DO
940!
941! Clear tangent linear forcing arrays before entering inner-loop.
942! This is very important since these arrays are non-zero after
943! running the representer model and must be zero when running the
944! tangent linear model.
945!
946 DO ng=1,ngrids
947 DO tile=first_tile(ng),last_tile(ng),+1
948 CALL initialize_forces (ng, tile, itlm)
949# ifdef ADJUST_BOUNDARY
950 CALL initialize_boundary (ng, tile, itlm)
951# endif
952 END DO
953 END DO
954
955# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
956!
957! Compute the reference zeta and biconjugate gradient arrays
958! required for the balance of free surface.
959!
960 IF (balance(isfsur)) THEN
961 DO ng=1,ngrids
962 CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
963 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
964!
965 DO tile=first_tile(ng),last_tile(ng),+1
966 CALL balance_ref (ng, tile, lini)
967 CALL biconj (ng, tile, inlm, lini)
968 END DO
969 wrtzetaref(ng)=.true.
970 END DO
971 END IF
972# endif
973!
974!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
975! 4D-Var inner loops.
976!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
977!
978 inner_loop : DO my_inner=0,ninner
979 inner=my_inner
980!
981! Initialize conjugate gradient algorithm depending on hot start or
982! outer loop index.
983!
984 IF (my_inner.eq.0) THEN
985 lcgini=.true.
986 DO ng=1,ngrids
987 CALL congrad (ng, irpm, my_outer, my_inner, ninner, lcgini)
988 END DO
989 END IF
990!
991! If initialization step, skip the inner-loop computations.
992!
993 linner=.false.
994 IF ((my_inner.ne.0).or.(nrun.ne.1)) THEN
995 IF (((my_inner.eq.0).and.lhotstart).or.(my_inner.ne.0)) THEN
996 linner=.true.
997 END IF
998 END IF
999!
1000! Start inner loop computations.
1001!
1002 inner_compute : IF (linner) THEN
1003!
1004!-----------------------------------------------------------------------
1005! Integrate adjoint model forced with any vector PSI at the observation
1006! locations and generate adjoint trajectory, Lambda_n(t).
1007!-----------------------------------------------------------------------
1008!
1009! Initialize the adjoint model from rest.
1010!
1011 DO ng=1,ngrids
1012 CALL ad_initial (ng)
1013 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1014 wrtmisfit(ng)=.false.
1015 END DO
1016
1017# ifdef RPM_RELAXATION
1018!
1019! Adjoint of representer relaxation is not applied during the
1020! inner-loops.
1021!
1022 DO ng=1,ngrids
1023 lweakrelax(ng)=.false.
1024 END DO
1025# endif
1026!
1027! Set adjoint history NetCDF parameters. Define adjoint history
1028! file only once to avoid opening too many files.
1029!
1030 DO ng=1,ngrids
1031 wrtforce(ng)=.true.
1032 IF (nrun.gt.1) ldefadj(ng)=.false.
1033 fcount=adm(ng)%load
1034 adm(ng)%Nrec(fcount)=0
1035 adm(ng)%Rindex=0
1036 END DO
1037!
1038! Time-step adjoint model backwards forced with current PSI vector.
1039!
1040 DO ng=1,ngrids
1041 IF (master) THEN
1042 WRITE (stdout,20) 'AD', ng, my_outer, my_inner, &
1043 & ntstart(ng), ntend(ng)
1044 END IF
1045 END DO
1046!
1047# ifdef SOLVE3D
1048 CALL ad_main3d (runinterval)
1049# else
1050 CALL ad_main2d (runinterval)
1051# endif
1052 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1053!
1054! Activate "LwrtState2d" switch to write the correct ad_ubar, ad_vbar,
1055! ad_zeta here instead of ad_ubar_sol, ad_vbar_sol, and ad_zeta_sol in
1056! the calls to "ad_wrt_his" below (HGA).
1057! Here, ad_zeta(:,:,kout)=ad_zeta_sol. However, ad_ubar and ad_vbar
1058! not equal to ad_ubar_sol and ad_vbar_sol, respectively. It is
1059! irrelevant because ubar and vbar are not part of the state in
1060! 3D application. Notice that "LwrtState2d" will be turned off at
1061! the bottom of "error_covariance".
1062!
1063 DO ng=1,ngrids
1064 lwrtstate2d(ng)=.true.
1065 END DO
1066!
1067! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
1068! record into the adjoint history file.
1069!
1070 DO ng=1,ngrids
1071# ifdef DISTRIBUTE
1072 CALL ad_wrt_his (ng, myrank)
1073# else
1074 CALL ad_wrt_his (ng, -1)
1075# endif
1076 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1077 END DO
1078!
1079! Write out adjoint initial condition record into the adjoint
1080! history file.
1081!
1082 DO ng=1,ngrids
1083 wrtforce(ng)=.false.
1084 lwrtstate2d(ng)=.true.
1085# ifdef DISTRIBUTE
1086 CALL ad_wrt_his (ng, myrank)
1087# else
1088 CALL ad_wrt_his (ng, -1)
1089# endif
1090 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1091 END DO
1092!
1093! Convolve adjoint trajectory with error covariances. Notice that we
1094! use iTLM to write the convolved solution into the ITL(ng)%name file
1095! (record Rec1) as the TLM initial conditions for the next inner loop.
1096!
1097# ifdef POSTERIOR_ERROR_I
1098 lposterior=.true.
1099# else
1100 lposterior=.false.
1101# endif
1102 CALL error_covariance (itlm, driver, my_outer, my_inner, &
1103 & lbck, lini, lold, lnew, &
1104 & rec1, rec2, lposterior)
1105 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1106!
1107! Convert the current adjoint solution in ADM(ng)%name to impulse
1108! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
1109! To facilitate the forcing to the TLM and RPM, the forcing is
1110! processed and written in increasing time coordinates (recall that
1111! the adjoint solution in ADM(ng)%name is backwards in time).
1112!
1113 IF (master) THEN
1114 WRITE (stdout,40) my_outer, my_inner
1115 END IF
1116 DO ng=1,ngrids
1117 tlf(ng)%Rindex=0
1118# ifdef DISTRIBUTE
1119 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1120# else
1121 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1122# endif
1123 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1124 END DO
1125!
1126!-----------------------------------------------------------------------
1127! Integrate tangent linear model forced by the convolved adjoint
1128! trajectory (impulse forcing) to compute R_n * PSI at observation
1129! points.
1130!-----------------------------------------------------------------------
1131!
1132 DO ng=1,ngrids
1133 wrtnlmod(ng)=.false.
1134 wrttlmod(ng)=.true.
1135 END DO
1136!
1137! If weak constraint, the impulses are time-interpolated at each
1138! time-steps.
1139!
1140 DO ng=1,ngrids
1141 IF (frcrec(ng).gt.3) THEN
1142 frequentimpulse(ng)=.true.
1143 END IF
1144 END DO
1145!
1146! Initialize tangent linear model from ITL(ng)%name, record Rec1.
1147!
1148 DO ng=1,ngrids
1149 itl(ng)%Rindex=rec1
1150 CALL tl_initial (ng)
1151 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1152 END DO
1153!
1154! Activate switch to write out initial misfit between model and
1155! observations.
1156!
1157 IF ((my_outer.eq.1).and.(my_inner.eq.1)) THEN
1158 DO ng=1,ngrids
1159 wrtmisfit(ng)=.true.
1160 END DO
1161 END IF
1162!
1163! Set tangent linear history NetCDF parameters. Define tangent linear
1164! history file at the beggining of each inner loop to avoid opening
1165! too many NetCDF files.
1166!
1167 DO ng=1,ngrids
1168 IF (my_inner.gt.1) ldeftlm(ng)=.false.
1169 fcount=tlm(ng)%load
1170 tlm(ng)%Nrec(fcount)=0
1171 tlm(ng)%Rindex=0
1172 END DO
1173!
1174! Run tangent linear model forward and force with convolved adjoint
1175! trajectory impulses. Compute R_n * PSI at observation points which
1176! are used in the conjugate gradient algorithm.
1177!
1178 DO ng=1,ngrids
1179 IF (master) THEN
1180 WRITE (stdout,20) 'TL', ng, my_outer, my_inner, &
1181 & ntstart(ng), ntend(ng)
1182 END IF
1183 END DO
1184!
1185# ifdef SOLVE3D
1186 CALL tl_main3d (runinterval)
1187# else
1188 CALL tl_main2d (runinterval)
1189# endif
1190 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1191!
1192 DO ng=1,ngrids
1193 wrtnlmod(ng)=.false.
1194 wrttlmod(ng)=.false.
1195 END DO
1196
1197# ifdef POSTERIOR_ERROR_F
1198!
1199! Copy the final time tl_var(Lold) into ad_var(Lold) so that it can be
1200! written to the Hessian NetCDF file.
1201!
1202 add=.false.
1203 DO ng=1,ngrids
1204 DO tile=first_tile(ng),last_tile(ng),+1
1205 CALL load_tltoad (ng, tile, lold(ng), lold(ng), add)
1206 END DO
1207 END DO
1208!
1209! Write evolved tangent solution into hessian netcdf file for use
1210! later.
1211!
1212 IF (my_inner.ne.0) THEN
1213 DO ng=1,ngrids
1214# ifdef DISTRIBUTE
1215 CALL wrt_hessian (ng, myrank, lold(ng), lold(ng))
1216# else
1217 CALL wrt_hessian (ng, -1, lold(ng), lold(ng))
1218# endif
1219 IF (founderror(exit_flag, noerror, &
1220 & __line__, myfile)) RETURN
1221 END DO
1222 END IF
1223# endif
1224!
1225!-----------------------------------------------------------------------
1226! Use conjugate gradient algorithm to find a better approximation
1227! PSI to representer coefficients Beta_n.
1228!-----------------------------------------------------------------------
1229!
1230 nrun=nrun+1
1231 DO ng=1,ngrids
1232 lcgini=.false.
1233 CALL congrad (ng, irpm, my_outer, my_inner, ninner, lcgini)
1234 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1235 END DO
1236
1237 END IF inner_compute
1238
1239 END DO inner_loop
1240!
1241! Close tangent linear NetCDF file.
1242!
1243 sourcefile=myfile
1244 DO ng=1,ngrids
1245 CALL close_file (ng, itlm, tlm(ng))
1246 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1247 END DO
1248!
1249!-----------------------------------------------------------------------
1250! Once that the representer coefficients, Beta_n, have been
1251! approximated with sufficient accuracy, compute estimates of
1252! Lambda_n and Xhat_n by carrying out one backward intergration
1253! of the adjoint model and one forward itegration of the representer
1254! model.
1255!-----------------------------------------------------------------------
1256!
1257! Initialize the adjoint model always from rest.
1258!
1259 DO ng=1,ngrids
1260 CALL ad_initial (ng)
1261 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1262 END DO
1263
1264# ifdef RPM_RELAXATION
1265!
1266! Adjoint of representer relaxation is applied during the
1267! outer-loops.
1268!
1269 DO ng=1,ngrids
1270 lweakrelax(ng)=.true.
1271 END DO
1272# endif
1273!
1274! Set adjoint history NetCDF parameters. Define adjoint history
1275! file one to avoid opening to many files.
1276!
1277 DO ng=1,ngrids
1278 wrtforce(ng)=.true.
1279 IF (nrun.gt.1) ldefadj(ng)=.false.
1280 fcount=adm(ng)%load
1281 adm(ng)%Nrec(fcount)=0
1282 adm(ng)%Rindex=0
1283 END DO
1284!
1285! Time-step adjoint model backwards forced with estimated representer
1286! coefficients, Beta_n.
1287!
1288 DO ng=1,ngrids
1289 IF (master) THEN
1290 WRITE (stdout,20) 'AD', ng, my_outer, ninner, &
1291 & ntstart(ng), ntend(ng)
1292 END IF
1293 END DO
1294!
1295# ifdef SOLVE3D
1296 CALL ad_main3d (runinterval)
1297# else
1298 CALL ad_main2d (runinterval)
1299# endif
1300
1301 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1302!
1303! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
1304! record into the adjoint history file. Note that the weak-constraint
1305! forcing is delayed by nADJ time-steps.
1306!
1307 DO ng=1,ngrids
1308# ifdef DISTRIBUTE
1309 CALL ad_wrt_his (ng, myrank)
1310# else
1311 CALL ad_wrt_his (ng, -1)
1312# endif
1313 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1314 END DO
1315!
1316! Write out adjoint initial condition record into the adjoint
1317! history file.
1318!
1319 DO ng=1,ngrids
1320 wrtforce(ng)=.false.
1321# ifdef DISTRIBUTE
1322 CALL ad_wrt_his (ng, myrank)
1323# else
1324 CALL ad_wrt_his (ng, -1)
1325# endif
1326 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1327 END DO
1328!
1329! Convolve adjoint trajectory with error covariances. Notice that we
1330! use iRPM to write the convolved solution into the IRP(ng)%name file
1331! (record Rec2) as the RPM initial conditions for the next outer loop.
1332!
1333 lposterior=.false.
1334 CALL error_covariance (irpm, driver, my_outer, inner, &
1335 & lbck, lini, lold, lnew, &
1336 & rec1, rec2, lposterior)
1337 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1338!
1339! Convert the current adjoint solution in ADM(ng)%name to impulse
1340! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
1341! To facilitate the forcing to the TLM and RPM, the forcing is
1342! processed and written in increasing time coordinates (recall that
1343! the adjoint solution in ADM(ng)%name is backwards in time).
1344!
1345 IF (master) THEN
1346 WRITE (stdout,40) my_outer, inner
1347 END IF
1348 DO ng=1,ngrids
1349 tlf(ng)%Rindex=0
1350# ifdef DISTRIBUTE
1351 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1352# else
1353 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1354# endif
1355 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1356 END DO
1357
1358# endif /* !DATALESS_LOOPS */
1359
1360# ifdef PROFILE
1361!
1362! Stop profile clock
1363!
1364 DO ng=1,ngrids
1365 DO thread=thread_range
1366 CALL wclock_off (ng, itlm, 87, __line__, myfile)
1367 END DO
1368 END DO
1369# endif
1370!
1371 10 FORMAT (a,'_outer',i0,'.nc')
1372 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
1373 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
1374 ', TimeSteps: ',i0,' - ',i0,')',/)
1375 30 FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ', &
1376 & 1p,e17.10,0p,t68,a)
1377 40 FORMAT (/,' Converting Convolved Adjoint Trajectory to', &
1378 & ' Impulses: Outer = ',i3.3,' Inner = ',i3.3,/)
1379!
1380 RETURN
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 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

References ad_def_his_mod::ad_def_his(), ad_initial(), ad_main2d(), ad_main3d(), ad_wrt_his_mod::ad_wrt_his(), mod_iounits::adm, mod_scalars::balance, zeta_balance_mod::balance_ref(), zeta_balance_mod::biconj(), congrad_mod::cg_read_congrad(), close_io_mod::close_file(), congrad_mod::congrad(), mod_iounits::dav, def_error_mod::def_error(), def_hessian_mod::def_hessian(), def_impulse_mod::def_impulse(), def_ini_mod::def_ini(), def_mod_mod::def_mod(), mod_scalars::dt, edit_multifile(), convolve_mod::error_covariance(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_scalars::frcrec, mod_scalars::frequentimpulse, mod_iounits::fwd, get_state_mod::get_state(), mod_fourdvar::havenlmod, mod_iounits::his, mod_param::iadm, mod_ncparam::ideftlm, mod_ncparam::idnlmi, mod_ncparam::idobss, mod_ncparam::idoerr, mod_ncparam::idoval, mod_ncparam::idsvar, mod_ncparam::idtime, mod_iounits::ini, mod_boundary::initialize_boundary(), mod_forces::initialize_forces(), mod_ocean::initialize_ocean(), mod_scalars::initime, mod_param::inlm, mod_scalars::inner, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_iounits::irp, mod_param::irpm, mod_ncparam::isfsur, mod_iounits::itl, mod_param::itlm, mod_parallel::last_tile, lbck, mod_scalars::ldefadj, mod_scalars::ldefavg, mod_scalars::ldeferr, mod_scalars::ldefhss, mod_scalars::ldefini, mod_scalars::ldefirp, mod_scalars::ldefitl, mod_scalars::ldefmod, mod_scalars::ldeftlf, mod_scalars::ldeftlm, mod_fourdvar::lhotstart, lini, mod_stepping::lnew, mod_stepping::lold, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadfwd, mod_fourdvar::lweakrelax, mod_scalars::lwrtavg, mod_scalars::lwrtstate2d, mod_scalars::lwrttlm, mod_parallel::master, mod_parallel::myrank, mod_scalars::nbrec, mod_scalars::nfrec, mod_param::ngrids, mod_scalars::ninner, mod_fourdvar::nlmodval, mod_scalars::nobc, mod_scalars::noerror, mod_scalars::nrun, mod_scalars::nsff, mod_fourdvar::nstatevar, mod_scalars::ntend, mod_scalars::ntstart, mod_scalars::obc_time, mod_iounits::obs, mod_fourdvar::obserr, mod_fourdvar::obsscale, mod_fourdvar::obsval, rec1, rec2, rp_def_ini_mod::rp_def_ini(), rp_initial(), rp_main2d(), rp_main3d(), mod_scalars::sf_time, mod_iounits::sourcefile, mod_scalars::sporadicimpulse, mod_iounits::stdout, tl_def_ini_mod::tl_def_ini(), tl_initial(), tl_main2d(), tl_main3d(), mod_iounits::tlf, mod_iounits::tlm, mod_ncparam::vname, wclock_off(), wclock_on(), wrt_hessian_mod::wrt_hessian(), wrt_impulse_mod::wrt_impulse(), mod_fourdvar::wrtforce, mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, mod_fourdvar::wrtrpmod, mod_fourdvar::wrttlmod, and mod_fourdvar::wrtzetaref.

Here is the call graph for this function:

◆ posterior_error()

subroutine, public r4dvar_mod::posterior_error ( real(dp), intent(in) runinterval)

Definition at line 1945 of file r4dvar.F.

1946!
1947!=======================================================================
1948! !
1949! This routine computes posterior analysis error covariance matrix !
1950! using the Lanczos vectors. !
1951! !
1952! On Input: !
1953! !
1954! RunInterval Kernel time stepping window (seconds) !
1955! !
1956! Warning: !
1957! !
1958! Currently, this code only works for a single outer-loop. !
1959! !
1960!=======================================================================
1961!
1962! Imported variable declarations
1963!
1964 real(dp), intent(in) :: RunInterval
1965!
1966! Local variable declarations.
1967!
1968 logical :: Ltrace
1969!
1970 integer :: my_inner, my_outer, ng, tile
1971 integer :: Fcount, Rec
1972!
1973 character (len=6) :: driver = 'w4dvar'
1974
1975 character (len=*), parameter :: MyFile = &
1976 & __FILE__//", posterior_error"
1977!
1978 sourcefile=myfile
1979
1980# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
1981!
1982!-----------------------------------------------------------------------
1983! Compute full (diagonal) posterior analysis error covariance matrix.
1984!-----------------------------------------------------------------------
1985!
1986! Clear tangent and adjoint arrays because they are used as
1987! work arrays below.
1988!
1989 DO ng=1,ngrids
1990 DO tile=first_tile(ng),last_tile(ng),+1
1991 CALL initialize_ocean (ng, tile, iadm)
1992 CALL initialize_ocean (ng, tile, itlm)
1993 CALL initialize_forces (ng, tile, iadm)
1994 CALL initialize_forces (ng, tile, itlm)
1995# ifdef ADJUST_BOUNDARY
1996 CALL initialize_boundary (ng, tile, iadm)
1997 CALL initialize_boundary (ng, tile, itlm)
1998# endif
1999 END DO
2000 END DO
2001!
2002! Compute the diagonal of the posterior/analysis error covariance
2003! matrix. The result is written to record 2 of the ITL netcdf file.
2004!
2005 var_oloop : DO my_outer=1,nouter
2006 outer=my_outer
2007 DO ng=1,ngrids
2008 DO tile=first_tile(ng),last_tile(ng),+1
2009 CALL posterior_var (ng, tile, itlm, my_outer)
2010 END DO
2011 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2012 END DO
2013 END DO var_oloop
2014!
2015! Write out the diagonal of the posterior/analysis covariance matrix
2016! which is in tl_var(Rec1) to 4DVar error NetCDF file.
2017!
2018 DO ng=1,ngrids
2019# ifdef DISTRIBUTE
2020 CALL wrt_error (ng, myrank, rec1, rec1)
2021# else
2022 CALL wrt_error (ng, -1, rec1, rec1)
2023# endif
2024 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2025 END DO
2026!
2027! Clear tangent and adjoint arrays because they are used as
2028! work arrays below.
2029!
2030 DO ng=1,ngrids
2031 DO tile=first_tile(ng),last_tile(ng),+1
2032 CALL initialize_ocean (ng, tile, iadm)
2033 CALL initialize_ocean (ng, tile, itlm)
2034 CALL initialize_forces (ng, tile, iadm)
2035 CALL initialize_forces (ng, tile, itlm)
2036# ifdef ADJUST_BOUNDARY
2037 CALL initialize_boundary (ng, tile, iadm)
2038 CALL initialize_boundary (ng, tile, itlm)
2039# endif
2040 END DO
2041 END DO
2042# endif
2043
2044# ifdef POSTERIOR_EOFS
2045!
2046!-----------------------------------------------------------------------
2047! Compute the posterior analysis error covariance matrix EOFs using a
2048! Lanczos algorithm.
2049!
2050! NOTE: Currently, this code only works for a single outer-loop.
2051!-----------------------------------------------------------------------
2052!
2053 IF (master) WRITE (stdout,10)
2054!
2055! Estimate first the trace of the posterior analysis error
2056! covariance matrix since the evolved and convolved Lanczos
2057! vectors stored in the Hessian NetCDF file will be destroyed
2058! later.
2059!
2060 ltrace=.true.
2061
2062 trace_oloop : DO my_outer=1,nouter
2063 outer=my_outer
2064 inner=0
2065
2066 trace_iloop : DO my_inner=1,nposti
2067 inner=my_inner
2068!
2069! Initialize the tangent linear variables with a random vector
2070! comprised of +1 and -1 elements randomly chosen.
2071!
2072 DO ng=1,ngrids
2073 DO tile=first_tile(ng),last_tile(ng),+1
2074 CALL random_ic (ng, tile, itlm, inner, outer, &
2075 & lold(ng), ltrace)
2076 END DO
2077 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2078 END DO
2079!
2080! Apply horizontal convolution.
2081!
2082 CALL convolve (driver, lini, lold, lnew)
2083 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2084!
2085! Compute Lanczos vector and eigenvectors of the posterior analysis
2086! error covariance matrix.
2087!
2088 DO ng=1,ngrids
2089 DO tile=first_tile(ng),last_tile(ng),+1
2090 CALL posterior (ng, tile, itlm, inner, outer, ltrace)
2091 END DO
2092 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2093 END DO
2094
2095 END DO trace_iloop
2096
2097 END DO trace_oloop
2098!
2099! Estimate posterior analysis error covariance matrix.
2100!
2101 ltrace=.false.
2102
2103 post_oloop : DO my_outer=1,nouter
2104 outer=my_outer
2105 inner=0
2106!
2107! The Lanczos algorithm requires to save all the Lanczos vectors.
2108! They are used to compute the posterior EOFs.
2109!
2110 DO ng=1,ngrids
2111 adm(ng)%Rindex=0
2112 fcount=adm(ng)%load
2113 adm(ng)%Nrec(fcount)=0
2114 END DO
2115
2116 post_iloop : DO my_inner=0,nposti
2117 inner=my_inner
2118!
2119! Read first record of ITL file and apply convolutions.
2120!
2121! NOTE: If inner=0, we would like to use a random starting vector.
2122! For now we can use what ever is in record 1.
2123!
2124 IF (my_inner.ne.0) THEN
2125 DO ng=1,ngrids
2126 rec=1
2127 CALL get_state (ng, itlm, 1, itl(ng), rec, lold(ng))
2128 IF (founderror(exit_flag, noerror, &
2129 & __line__, myfile)) RETURN
2130 END DO
2131 ELSE
2132 DO ng=1,ngrids
2133 DO tile=first_tile(ng),last_tile(ng),+1
2134 CALL random_ic (ng, tile, itlm, inner, outer, &
2135 & lold(ng), ltrace)
2136 END DO
2137 IF (founderror(exit_flag, noerror, &
2138 & __line__, myfile)) RETURN
2139 END DO
2140 END IF
2141!
2142! Apply horizontal convolution.
2143!
2144 CALL convolve (driver, lini, lold, lnew)
2145 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2146!
2147! Compute Lanczos vector and eigenvectors of the posterior analysis
2148! error covariance matrix.
2149!
2150 DO ng=1,ngrids
2151 DO tile=first_tile(ng),last_tile(ng),+1
2152 CALL posterior (ng, tile, itlm, inner, outer, ltrace)
2153 END DO
2154 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2155 END DO
2156!
2157! Write the Lanczos vectors of the posterior error covariance
2158! to the adjoint NetCDF file.
2159!
2160 DO ng=1,ngrids
2161# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
2162 lfout(ng)=lnew(ng)
2163# endif
2164# ifdef ADJUST_BOUNDARY
2165 lbout(ng)=lnew(ng)
2166# endif
2167 kstp(ng)=lnew(ng)
2168# ifdef SOLVE3D
2169 nstp(ng)=lnew(ng)
2170# endif
2171 lwrtstate2d(ng)=.true.
2172# ifdef DISTRIBUTE
2173 CALL ad_wrt_his (ng, myrank)
2174# else
2175 CALL ad_wrt_his (ng, -1)
2176# endif
2177 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2178 lwrtstate2d(ng)=.false.
2179 END DO
2180!
2181! Write out tangent linear model initial conditions and tangent
2182! linear surface forcing adjustments for next inner loop into
2183! ITL(ng)%name (record Rec1).
2184!
2185 DO ng=1,ngrids
2186# ifdef DISTRIBUTE
2187 CALL tl_wrt_ini (ng, myrank, lnew(ng), rec1)
2188# else
2189 CALL tl_wrt_ini (ng, -1, lnew(ng), rec1)
2190# endif
2191 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2192 END DO
2193
2194 END DO post_iloop
2195
2196 END DO post_oloop
2197
2198# endif /* POSTERIOR_EOFS */
2199!
2200! Done. Set history file ID to closed state since we manipulated
2201! its indices with the forward file ID which was closed above.
2202!
2203 DO ng=1,ngrids
2204 his(ng)%ncid=-1
2205 END DO
2206
2207# ifdef POSTERIOR_EOFS
2208!
2209 10 FORMAT (/,' <<<< Posterior Analysis Error Covariance Matrix', &
2210 & ' Estimation >>>>',/)
2211# endif
2212!
2213 RETURN

References ad_wrt_his_mod::ad_wrt_his(), mod_iounits::adm, convolve_mod::convolve(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), get_state_mod::get_state(), mod_iounits::his, mod_param::iadm, mod_boundary::initialize_boundary(), mod_forces::initialize_forces(), mod_ocean::initialize_ocean(), mod_scalars::inner, mod_iounits::itl, mod_param::itlm, mod_stepping::kstp, mod_parallel::last_tile, mod_stepping::lbout, mod_stepping::lfout, lini, mod_stepping::lnew, mod_stepping::lold, mod_scalars::lwrtstate2d, mod_parallel::master, mod_parallel::myrank, mod_param::ngrids, mod_scalars::noerror, mod_scalars::nouter, mod_fourdvar::nposti, mod_stepping::nstp, mod_scalars::outer, posterior_mod::posterior(), posterior_var_mod::posterior_var(), random_ic_mod::random_ic(), rec1, mod_iounits::sourcefile, mod_iounits::stdout, and tl_wrt_ini_mod::tl_wrt_ini().

Here is the call graph for this function:

◆ prior_error()

subroutine, public r4dvar_mod::prior_error ( integer, intent(in) ng)

Definition at line 1781 of file r4dvar.F.

1782!
1783!=======================================================================
1784! !
1785! This routine processes background prior error covariance standard !
1786! deviations and normalization coefficients. !
1787! !
1788! On Input: !
1789! !
1790! ng Nested grid number !
1791! !
1792!=======================================================================
1793!
1794! Imported variable declarations
1795!
1796 integer, intent(in) :: ng
1797!
1798! Local variable declarations.
1799!
1800 integer :: tile
1801 integer :: NRMrec, STDrec, Tindex
1802!
1803 character (len=*), parameter :: MyFile = &
1804 & __FILE__//", prior_error"
1805!
1806 sourcefile=myfile
1807!
1808!-----------------------------------------------------------------------
1809! Set application grid, metrics, and associated variables and
1810! parameters.
1811!-----------------------------------------------------------------------
1812!
1813! The ROMS application grid configuration is done once. It is usually
1814! done in the "initial" kernel routine. However, since we are calling
1815! the "normalization" routine here, we need several grid variables and
1816! parameter. Also, if reading only water points, we need to know the
1817! land/sea mask arrays to unpack.
1818!
1819 IF (setgridconfig(ng)) THEN
1820 CALL set_grid (ng, inlm)
1821 END IF
1822!
1823!-----------------------------------------------------------------------
1824! Read in standard deviation factors for error covariance.
1825!-----------------------------------------------------------------------
1826!
1827! Initial conditions standard deviation. They are loaded in Tindex=1
1828! of the e_var(...,Tindex) state variables.
1829!
1830 IF (lreadstd(ng)) THEN
1831 stdrec=1
1832 tindex=1
1833# ifdef STD_MODEL
1834 CALL get_state (ng, 10, 10, std(5,ng), stdrec, tindex)
1835# else
1836 CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
1837# endif
1838 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1839 END IF
1840!
1841! Model error standard deviation. They are loaded in Tindex=2
1842! of the e_var(...,Tindex) state variables.
1843!
1844 stdrec=1
1845 tindex=2
1846 IF (nsa.eq.2) THEN
1847 CALL get_state (ng, 11, 11, std(2,ng), stdrec, tindex)
1848 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1849 END IF
1850
1851# ifdef ADJUST_BOUNDARY
1852!
1853! Open boundary conditions standard deviation.
1854!
1855 stdrec=1
1856 tindex=1
1857 CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
1858 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1859# endif
1860
1861# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
1862!
1863! Surface forcing standard deviation.
1864!
1865 stdrec=1
1866 tindex=1
1867 CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
1868 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1869# endif
1870
1871# ifdef STD_MODEL
1872!
1873! If computing the standard deviation (STD) directly from the
1874! background (prior) field as an alternative to climatological
1875! read from the input NetCDF file, define the output STD NetCDF
1876! file.
1877!
1878 IF (ldefstd(ng)) THEN
1879 CALL def_std (ng)
1880 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1881 END IF
1882# endif
1883!
1884!-----------------------------------------------------------------------
1885! Error covariance normalization coefficients.
1886!-----------------------------------------------------------------------
1887!
1888! Compute or read in the error covariance normalization factors.
1889! If computing, write out factors to NetCDF. This is an expensive
1890! computation that needs to be computed only once for a particular
1891! application grid and decorrelation scales.
1892!
1893 IF (any(lwrtnrm(:,ng))) THEN
1894 CALL def_norm (ng, inlm, 1)
1895 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1896
1897 IF (nsa.eq.2) THEN
1898 CALL def_norm (ng, inlm, 2)
1899 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1900 END IF
1901
1902# ifdef ADJUST_BOUNDARY
1903 CALL def_norm (ng, inlm, 3)
1904 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1905# endif
1906
1907# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
1908 CALL def_norm (ng, inlm, 4)
1909 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1910# endif
1911!
1912 DO tile=first_tile(ng),last_tile(ng),+1
1913 CALL normalization (ng, tile, 2)
1914 END DO
1915 ldefnrm(1:4,ng)=.false.
1916 lwrtnrm(1:4,ng)=.false.
1917
1918 ELSE
1919 nrmrec=1
1920 CALL get_state (ng, 14, 14, nrm(1,ng), nrmrec, 1)
1921 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1922
1923 IF (nsa.eq.2) THEN
1924 CALL get_state (ng, 15, 15, nrm(2,ng), nrmrec, 2)
1925 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1926 END IF
1927
1928# ifdef ADJUST_BOUNDARY
1929 CALL get_state (ng, 16, 16, nrm(3,ng), nrmrec, 1)
1930 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1931# endif
1932
1933# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
1934 CALL get_state (ng, 17, 17, nrm(4,ng), nrmrec, 1)
1935 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1936# endif
1937 END IF
1938!
1939 RETURN
subroutine set_grid(ng, model)
Definition set_grid.F:3

References def_norm_mod::def_norm(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), get_state_mod::get_state(), mod_param::inlm, mod_parallel::last_tile, mod_scalars::ldefnrm, mod_scalars::lreadstd, mod_scalars::lwrtnrm, mod_scalars::noerror, normalization_mod::normalization(), mod_iounits::nrm, mod_param::nsa, posterior_mod::posterior_eofs(), set_grid(), mod_scalars::setgridconfig, mod_iounits::sourcefile, and mod_iounits::std.

Here is the call graph for this function:

Variable Documentation

◆ lbck

integer r4dvar_mod::lbck = 2

Definition at line 151 of file r4dvar.F.

151 integer :: Lbck = 2 ! background record in INI

Referenced by increment().

◆ ldone

logical r4dvar_mod::ldone

Definition at line 147 of file r4dvar.F.

147 logical :: Ldone

◆ lini

integer r4dvar_mod::lini = 1

Definition at line 150 of file r4dvar.F.

150 integer :: Lini = 1 ! NLM initial conditions record in INI

Referenced by increment(), and posterior_error().

◆ rec1

integer r4dvar_mod::rec1 = 1

Definition at line 152 of file r4dvar.F.

152 integer :: Rec1 = 1

Referenced by analysis(), increment(), and posterior_error().

◆ rec2

integer r4dvar_mod::rec2 = 2

Definition at line 153 of file r4dvar.F.

153 integer :: Rec2 = 2

Referenced by analysis(), and increment().