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