ROMS
Loading...
Searching...
No Matches
propagator_hso.h
Go to the documentation of this file.
1 MODULE propagator_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! Stochastic Optimals Propagator for white noise forcing. !
11! !
12! Reference: !
13! !
14! Moore, A.M. et al., 2004: A comprehensive ocean prediction and !
15! analysis system based on the tangent linear and adjoint of a !
16! regional ocean model, Ocean Modelling, 7, 227-258. !
17! !
18! Smith, K. D., A. M. Moore, and H. G. Arango, 2015: Estimates of !
19! ocean forecast error covariance derived from Hessian Singular !
20! Vectors, Ocean Modelling, 89, 104-121, !
21! doi:10.1016/j.ocemod.2015.03.003. !
22! !
23!=======================================================================
24!
25 USE mod_kinds
26!
27 implicit none
28!
29 PRIVATE
30 PUBLIC :: propagator_hso
31!
32 CONTAINS
33!
34!***********************************************************************
35 SUBROUTINE propagator_hso (RunInterval, Iter, state, ad_state)
36!***********************************************************************
37!
38 USE mod_param
39 USE mod_parallel
40 USE mod_netcdf
41#ifdef SOLVE3D
42 USE mod_coupling
43#endif
44 USE mod_iounits
45 USE mod_ocean
46 USE mod_scalars
47 USE mod_stepping
48 USE mod_storage
49!
50 USE close_io_mod, ONLY : close_inp
55 USE inner2state_mod, ONLY : ini_c_norm
56#ifdef STOCH_OPT_WHITE
57 USE packing_mod, ONLY : ad_so_pack, ad_unpack, tl_unpack
58#else
59 USE packing_mod, ONLY : ad_so_pack_red, ad_unpack, tl_unpack
60#endif
61#ifdef SOLVE3D
62 USE set_depth_mod, ONLY : set_depth
63#endif
64 USE strings_mod, ONLY : founderror
65!
66! Imported variable declarations.
67!
68 integer :: iter
69!
70 real(dp), intent(in) :: runinterval
71!
72 TYPE (t_gst), intent(in) :: state(ngrids)
73 TYPE (t_gst), intent(inout) :: ad_state(ngrids)
74!
75! Local variable declarations.
76!
77 logical :: soruntl
78#ifdef STOCH_OPT_WHITE
79 logical :: sorunad
80#endif
81!
82 integer :: ng, tile
83 integer :: ktmp, ntmp, lini
84 integer :: kout, nout
85 integer :: fcount, interval, inttrap
86!
87 real(r8) :: statenorm(ngrids)
88 real(r8) :: so_run_time
89!
90 character (len=*), parameter :: myfile = &
91 & __FILE__
92!
93!=======================================================================
94! Forward integration of the tangent linear model.
95!=======================================================================
96!
97 nrun=nrun+1
98 IF (master) THEN
99 DO ng=1,ngrids
100 WRITE (stdout,10) ' PROPAGATOR - Grid: ', ng, &
101 & ', Iteration: ', nrun, &
102 & ', number converged RITZ values: ', &
103 & nconv(ng)
104 END DO
105 END IF
106!
107! Loop over the required numger if trapezoidal intervals in time.
108!
109 interval_loop : DO interval=1,nintervals+1
110
111 soruntl=.true.
112#ifdef STOCH_OPT_WHITE
113 sorunad=.true.
114#endif
115 IF (interval.eq.nintervals+1) THEN
116 soruntl=.false.
117#ifdef STOCH_OPT_WHITE
118 sorunad=.false.
119#endif
120 END IF
121 inttrap=interval
122!
123 IF (master) THEN
124 WRITE (stdout,20) ' Stochastic Optimals Time Interval = ', &
125 & interval
126 END IF
127!
128! Initialize time stepping indices and counters.
129!
130 lini=1
131 DO ng=1,ngrids
132#ifndef STOCH_OPT_WHITE
133 IF (interval.eq.1) THEN
134 soinitial(ng)=.true.
135 ELSE
136 soinitial(ng)=.false.
137 END IF
138#endif
139 IF (interval.eq.1) THEN
140 lwrttlm(ng)=.true.
141 ELSE
142 lwrttlm(ng)=.false.
143 END IF
144 iif(ng)=1
145 indx1(ng)=1
146 kstp(ng)=1
147 krhs(ng)=1
148 knew(ng)=1
149 predictor_2d_step(ng)=.false.
150!
151 iic(ng)=0
152 nstp(ng)=1
153 nrhs(ng)=1
154 nnew(ng)=1
155!
156 synchro_flag(ng)=.true.
157 tdays(ng)=dstart+real(ntimes(ng),r8)*real(interval-1,r8)* &
158 & dt(ng)*sec2day/real(nintervals,r8)
159 time(ng)=tdays(ng)*day2sec
160 ntstart(ng)=int((time(ng)-dstart*day2sec)/dt(ng))+1
161 ntend(ng)=ntimes(ng)
162 ntfirst(ng)=ntstart(ng)
163 so_run_time=dt(ng)*real(ntend(ng)-ntstart(ng)+1,r8)
164 END DO
165!
166! Set switches and counters to manage output adjoint and tangent linear
167! history NetCDF files.
168!
169 DO ng=1,ngrids
170 IF (iter.gt.0) THEN ! Arnoldi iteration loop
171 IF ((iter.eq.1).and.(interval.eq.1)) THEN
172 ldefadj(ng)=.true.
173 ldeftlm(ng)=.true. ! NetCDF files are created
174 ELSE ! on the first pass
175 ldefadj(ng)=.false.
176 ldeftlm(ng)=.false.
177 END IF
178 fcount=adm(ng)%load
179 adm(ng)%Nrec(fcount)=0
180 adm(ng)%Rindex=0
181 fcount=tlm(ng)%load
182 tlm(ng)%Nrec(fcount)=0
183 tlm(ng)%Rindex=0
184 ELSE ! Computing eigenvectors
185 IF (lmultigst.and.(interval.eq.1)) THEN
186 ldeftlm(ng)=.true.
187 ELSE
188 ldeftlm(ng)=.false.
189 END IF
190#ifdef STOCH_OPT_WHITE
191 IF (interval.le.nintervals) THEN
192 fcount=adm(ng)%load
193 adm(ng)%Nrec(fcount)=0
194 adm(ng)%Rindex=0
195 END IF
196#else
197 fcount=adm(ng)%load
198 adm(ng)%Nrec(fcount)=0
199 adm(ng)%Rindex=0
200#endif
201 IF ((lmultigst.or.(abs(iter).eq.1)).and. &
202 & (interval.eq.1)) THEN
203 fcount=tlm(ng)%load
204 tlm(ng)%Nrec(fcount)=0
205 tlm(ng)%Rindex=0
206 END IF
207 END IF
208 END DO
209!
210!-----------------------------------------------------------------------
211! Clear tangent linear state variables. There is not need to clean
212! the basic state arrays since they were zeroth out at initialization
213! and bottom of previous iteration.
214!-----------------------------------------------------------------------
215!
216 DO ng=1,ngrids
217 DO tile=first_tile(ng),last_tile(ng),+1
218 CALL initialize_ocean (ng, tile, itlm)
219 END DO
220 END DO
221
222#ifdef SOLVE3D
223!
224!-----------------------------------------------------------------------
225! Compute basic state initial level thicknesses used for state norm
226! scaling. It uses zero time-averaged free-surface (rest state).
227! Therefore, the norm scaling is time invariant.
228!-----------------------------------------------------------------------
229!
230 DO ng=1,ngrids
231 DO tile=last_tile(ng),first_tile(ng),-1
232 CALL set_depth (ng, tile, itlm)
233 END DO
234 END DO
235#endif
236!
237!-----------------------------------------------------------------------
238! Compute tangent linear initial conditions from state vector.
239!-----------------------------------------------------------------------
240!
241 DO ng=1,ngrids
242 DO tile=first_tile(ng),last_tile(ng),+1
243 CALL tl_inner2state (ng, tile, lini, state(ng)%vector)
244 END DO
245 END DO
246!
247!-----------------------------------------------------------------------
248! Compute initial tangent linear state analysis error norm.
249!-----------------------------------------------------------------------
250!
251 IF (interval.eq.1) THEN
252 DO ng=1,ngrids
253 DO tile=last_tile(ng),first_tile(ng),-1
254 CALL ini_c_norm (ng, tile, kstp(ng), nstp(ng), &
255 & statenorm(ng))
256 END DO
257 IF (master) THEN
258 WRITE (stdout,30) ' PROPAGATOR - Grid: ', ng, &
259 & ', Tangent Initial Norm: ', &
260 & statenorm(ng)
261 END IF
262 END DO
263 END IF
264!
265!-----------------------------------------------------------------------
266! Read in initial forcing, climatology and assimilation data from
267! input NetCDF files. It loads the first relevant data record for
268! the time-interpolation between snapshots.
269!-----------------------------------------------------------------------
270!
271 IF (soruntl) THEN ! do not run TLM on last interval
272 DO ng=1,ngrids
273 CALL close_inp (ng, itlm)
274 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
275
276 CALL tl_get_idata (ng)
277 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
278
279 CALL tl_get_data (ng)
280 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
281 END DO
282!
283!-----------------------------------------------------------------------
284! Time-step the tangent linear model.
285!-----------------------------------------------------------------------
286!
287 DO ng=1,ngrids
288 IF (master) THEN
289 WRITE (stdout,40) 'TL', ng, ntstart(ng), ntend(ng)
290 END IF
291 time(ng)=time(ng)-dt(ng)
292 iic(ng)=ntstart(ng)-1
293 END DO
294
295#ifdef SOLVE3D
296 CALL tl_main3d (so_run_time)
297#else
298 CALL tl_main2d (so_run_time)
299#endif
300 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
301 END IF
302!
303!-----------------------------------------------------------------------
304! Clear nonlinear (basic state) and adjoint state variables.
305!-----------------------------------------------------------------------
306!
307 DO ng=1,ngrids
308 DO tile=first_tile(ng),last_tile(ng),+1
309 CALL initialize_ocean (ng, tile, inlm)
310 CALL initialize_ocean (ng, tile, iadm)
311#ifdef SOLVE3D
312 CALL initialize_coupling (ng, tile, 0)
313#endif
314 END DO
315 END DO
316
317#ifdef SOLVE3D
318!
319!-----------------------------------------------------------------------
320! Compute basic state final level thicknesses used for state norm
321! scaling. It uses zero time-averaged free-surface (rest state).
322! Therefore, the norm scaling is time invariant.
323!-----------------------------------------------------------------------
324!
325 DO ng=1,ngrids
326 DO tile=last_tile(ng),first_tile(ng),-1
327 CALL set_depth (ng, tile, itlm)
328 END DO
329 END DO
330#endif
331!
332!-----------------------------------------------------------------------
333! Compute final tangent linear energy norm.
334!-----------------------------------------------------------------------
335!
336 IF (interval.eq.nintervals+1) THEN
337 DO ng=1,ngrids
338 DO tile=first_tile(ng),last_tile(ng),+1
339 CALL tl_statenorm (ng, tile, kstp(ng), nstp(ng), &
340 & statenorm(ng))
341 END DO
342 IF (master) THEN
343 WRITE (stdout,30) ' PROPAGATOR - Grid: ', ng, &
344 & ', Tangent Final Norm: ', &
345 & statenorm(ng)
346 END IF
347 END DO
348 END IF
349!
350!=======================================================================
351! Backward integration with the adjoint model.
352!=======================================================================
353!
354! Initialize time stepping indices and counters.
355!
356 DO ng=1,ngrids
357#ifndef STOCH_OPT_WHITE
358 lwrtstate2d(ng)=.false.
359 lwrtstate3d(ng)=.false.
360#endif
361 iif(ng)=1
362 indx1(ng)=1
363 ktmp=knew(ng)
364 IF (inttrap.eq.nintervals+1) THEN
365 ktmp=kstp(ng)
366 END IF
367 kstp(ng)=1
368 krhs(ng)=3
369 knew(ng)=2
370 kout=knew(ng)
371#ifdef STOCH_OPT_WHITE
372 IF (inttrap.eq.nintervals+1) THEN
373 kout=kstp(ng)
374 END IF
375#endif
376 predictor_2d_step(ng)=.false.
377!
378 iic(ng)=0
379 ntmp=nstp(ng)
380 nstp(ng)=1
381 nrhs(ng)=1
382 nnew(ng)=2
383!
384 synchro_flag(ng)=.true.
385 tdays(ng)=dstart+dt(ng)*real(ntimes(ng),r8)*sec2day
386 time(ng)=tdays(ng)*day2sec
387 ntstart(ng)=ntimes(ng)+1
388# ifdef STOCH_OPT_WHITE
389 ntend(ng)=1+(interval-1)*ntimes(ng)/nintervals
390# else
391 ntend(ng)=1
392# endif
393 ntfirst(ng)=ntend(ng)
394 END DO
395!
396!-----------------------------------------------------------------------
397! Initialize adjoint model with the final tangent linear solution
398! scaled by the energy norm.
399!-----------------------------------------------------------------------
400!
401 DO ng=1,ngrids
402 DO tile=last_tile(ng),first_tile(ng),-1
403 CALL ad_ini_perturb (ng, tile, &
404 & ktmp, kout, ntmp, nstp(ng))
405 END DO
406 END DO
407!
408!-----------------------------------------------------------------------
409! Read in initial forcing, climatology and assimilation data from
410! input NetCDF files. It loads the first relevant data record for
411! the time-interpolation between snapshots.
412!-----------------------------------------------------------------------
413!
414#ifdef STOCH_OPT_WHITE
415 IF (sorunad) THEN ! do not run ADM on last interval
416#endif
417 DO ng=1,ngrids
418 CALL close_inp (ng, iadm)
419 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
420
421 CALL ad_get_idata (ng)
422 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
423
424 CALL ad_get_data (ng)
425 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
426 END DO
427!
428!-----------------------------------------------------------------------
429! Time-step the adjoint model backwards.
430!-----------------------------------------------------------------------
431!
432 DO ng=1,ngrids
433 IF (master) THEN
434 WRITE (stdout,40) 'AD', ng, ntstart(ng), ntend(ng)
435 END IF
436 time(ng)=time(ng)+dt(ng)
437 iic(ng)=ntstart(ng)+1
438 END DO
439
440#ifdef SOLVE3D
441# ifdef STOCH_OPT_WHITE
442 CALL ad_main3d (so_run_time)
443# else
444 CALL ad_main3d (runinterval)
445# endif
446#else
447# ifdef STOCH_OPT_WHITE
448 CALL ad_main2d (so_run_time)
449# else
450 CALL ad_main2d (runinterval)
451# endif
452#endif
453 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
454#ifdef STOCH_OPT_WHITE
455 END IF
456#endif
457!
458!-----------------------------------------------------------------------
459! Clear nonlinear state (basic state) variables for next iteration
460! and to insure a rest state time averaged free-surface before adjoint
461! state norm scaling.
462!-----------------------------------------------------------------------
463!
464 DO ng=1,ngrids
465 DO tile=first_tile(ng),last_tile(ng),+1
466 CALL initialize_ocean (ng, tile, inlm)
467#ifdef SOLVE3D
468 CALL initialize_coupling (ng, tile, 0)
469#endif
470 END DO
471 END DO
472
473#ifdef SOLVE3D
474!
475!-----------------------------------------------------------------------
476! Compute basic state initial level thicknesses used for state norm
477! scaling. It uses zero free-surface (rest state). Therefore, the
478! norm scaling is time invariant.
479!-----------------------------------------------------------------------
480!
481 DO ng=1,ngrids
482 DO tile=last_tile(ng),first_tile(ng),-1
483 CALL set_depth (ng, tile, iadm)
484 END DO
485 END DO
486#endif
487!
488!-----------------------------------------------------------------------
489! Pack final adjoint solution into adjoint state vector.
490!-----------------------------------------------------------------------
491!
492 DO ng=1,ngrids
493 DO tile=first_tile(ng),last_tile(ng),+1
494# ifdef STOCH_OPT_WHITE
495 CALL ad_so_pack (ng, tile, nstr(ng), nend(ng), inttrap, &
496 & storage(ng)%my_state)
497# else
498 CALL ad_so_pack_red (ng, tile, nstr(ng), nend(ng), &
499 & inttrap, storage(ng)%my_state)
500# endif
501 END DO
502 END DO
503 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
504!
505!-----------------------------------------------------------------------
506! Clear forcing variables for next iteration.
507!-----------------------------------------------------------------------
508!
509 DO ng=1,ngrids
510 DO tile=last_tile(ng),first_tile(ng),-1
511 CALL initialize_forces (ng, tile, itlm)
512 CALL initialize_forces (ng, tile, iadm)
513 END DO
514 END DO
515
516 END DO interval_loop
517!
518!-----------------------------------------------------------------------
519! Finally, unpack "my_state", and compute the adjoint state vector.
520!-----------------------------------------------------------------------
521!
522 DO ng=1,ngrids
523 DO tile=first_tile(ng),last_tile(ng),+1
524 CALL ad_unpack (ng, tile, nstr(ng), nend(ng), &
525 & storage(ng)%my_state)
526 END DO
527 END DO
528 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
529!
530 DO ng=1,ngrids
531 DO tile=last_tile(ng),first_tile(ng),-1
532 CALL ad_inner2state (ng, tile, lini, ad_state(ng)%vector)
533 END DO
534 END DO
535 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
536!
537 10 FORMAT (/,a,i2.2,a,i3.3,a,i3.3/)
538 20 FORMAT (/,a,i2.2)
539 30 FORMAT (/,a,i2.2,a,1p,e15.6,/)
540 40 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
541 & ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')')
542!
543 RETURN
544 END SUBROUTINE propagator_hso
545
546 END MODULE propagator_mod
subroutine ad_get_data(ng)
Definition ad_get_data.F:4
subroutine ad_get_idata(ng)
Definition ad_get_idata.F:4
subroutine ad_main2d
Definition ad_main2d.F:586
subroutine ad_main3d(runinterval)
Definition ad_main3d.F:4
subroutine, public close_inp(ng, model)
Definition close_io.F:92
subroutine, public tl_statenorm(ng, tile, kinp, ninp, statenorm)
Definition dotproduct.F:699
subroutine, public ad_ini_perturb(ng, tile, kinp, kout, ninp, nout)
subroutine, public tl_inner2state(ng, tile, lini, state)
Definition inner2state.F:63
subroutine, public ad_inner2state(ng, tile, lini, ad_state)
subroutine, public ini_c_norm(ng, tile, kinp, ninp, statenorm)
subroutine, public initialize_coupling(ng, tile, model)
subroutine, public initialize_forces(ng, tile, model)
type(t_io), dimension(:), allocatable adm
type(t_io), dimension(:), allocatable tlm
integer stdout
integer, parameter r8
Definition mod_kinds.F:28
integer, parameter dp
Definition mod_kinds.F:25
subroutine, public initialize_ocean(ng, tile, model)
Definition mod_ocean.F:1526
integer, dimension(:), allocatable first_tile
logical master
integer, dimension(:), allocatable last_tile
integer, parameter inlm
Definition mod_param.F:662
integer, dimension(:), allocatable nstr
Definition mod_param.F:646
integer, parameter iadm
Definition mod_param.F:665
integer ngrids
Definition mod_param.F:113
integer, parameter itlm
Definition mod_param.F:663
integer, dimension(:), allocatable nend
Definition mod_param.F:647
real(dp), parameter day2sec
integer, dimension(:), allocatable ntimes
integer, dimension(:), allocatable iic
real(dp), dimension(:), allocatable dt
logical lmultigst
integer, dimension(:), allocatable nconv
logical, dimension(:), allocatable synchro_flag
logical, dimension(:), allocatable predictor_2d_step
real(dp), dimension(:), allocatable tdays
logical, dimension(:), allocatable ldefadj
real(dp) dstart
real(dp), parameter sec2day
integer, dimension(:), allocatable ntend
integer exit_flag
logical, dimension(:), allocatable lwrtstate2d
integer, dimension(:), allocatable indx1
integer nintervals
integer, dimension(:), allocatable ntfirst
logical, dimension(:), allocatable lwrttlm
real(dp), dimension(:), allocatable time
logical, dimension(:), allocatable ldeftlm
integer, dimension(:), allocatable ntstart
integer nrun
integer, dimension(:), allocatable iif
integer noerror
integer, dimension(:), allocatable kstp
integer, dimension(:), allocatable knew
integer, dimension(:), allocatable nrhs
integer, dimension(:), allocatable nnew
integer, dimension(:), allocatable krhs
integer, dimension(:), allocatable nstp
type(t_storage), dimension(:), allocatable storage
Definition mod_storage.F:91
subroutine, public propagator_hso(runinterval, iter, state, ad_state)
subroutine, public set_depth(ng, tile, model)
Definition set_depth.F:34
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
subroutine ad_unpack(ng, tile, mstr, mend, state)
Definition packing.F:4712
subroutine tl_unpack(ng, tile, mstr, mend, state)
Definition packing.F:6561
subroutine tl_get_data(ng)
Definition tl_get_data.F:4
subroutine tl_get_idata(ng)
Definition tl_get_idata.F:4
subroutine tl_main2d
Definition tl_main2d.F:429
subroutine tl_main3d(runinterval)
Definition tl_main3d.F:4