ROMS
Loading...
Searching...
No Matches
propagator_op.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! Optimal Perturbation (singular vectors) Propagator: !
11! !
12! This routine solves the singular vectors of the propagator R(0,t) !
13! which measure the fastest growing of all possible perturbations !
14! over a given time interval. It involves a single integration of a !
15! perturbation forward in time with the tangent linear model over !
16! [0,t], followed by an integration of the result backward in time !
17! with the adjoint model over [t,0]. !
18! !
19! Reference: !
20! !
21! Moore, A.M. et al., 2004: A comprehensive ocean prediction and !
22! analysis system based on the tangent linear and adjoint of a !
23! regional ocean model, Ocean Modelling, 7, 227-258. !
24! !
25!=======================================================================
26!
27 USE mod_kinds
28!
29 implicit none
30!
31 PRIVATE
32 PUBLIC :: propagator_op
33!
34 CONTAINS
35!
36!***********************************************************************
37 SUBROUTINE propagator_op (RunInterval, state, ad_state)
38!***********************************************************************
39!
40 USE mod_param
41 USE mod_parallel
42#ifdef SOLVE3D
43 USE mod_coupling
44#endif
45 USE mod_iounits
46 USE mod_ocean
47 USE mod_scalars
48 USE mod_stepping
49!
50 USE close_io_mod, ONLY : close_inp
53 USE packing_mod, ONLY : tl_unpack, ad_pack
54#ifdef SOLVE3D
55 USE set_depth_mod, ONLY : set_depth
56#endif
57 USE strings_mod, ONLY : founderror
58!
59! Imported variable declarations.
60!
61 real(dp), intent(in) :: runinterval
62!
63 TYPE (t_gst), intent(in) :: state(ngrids)
64 TYPE (t_gst), intent(inout) :: ad_state(ngrids)
65!
66! Local variable declarations.
67!
68 integer :: ng, tile
69 integer :: ktmp, ntmp
70!
71 real(r8) :: statenorm(ngrids)
72!
73 character (len=*), parameter :: myfile = &
74 & __FILE__
75!
76!=======================================================================
77! Forward integration of the tangent linear model.
78!=======================================================================
79!
80 nrun=nrun+1
81 IF (master) THEN
82 DO ng=1,ngrids
83 WRITE (stdout,10) ' PROPAGATOR - Grid: ', ng, &
84 & ', Iteration: ', nrun, &
85 & ', number converged RITZ values: ', &
86 & nconv(ng)
87 END DO
88 END IF
89!
90! Initialize time stepping indices and counters.
91!
92 DO ng=1,ngrids
93 iif(ng)=1
94 indx1(ng)=1
95 kstp(ng)=1
96 krhs(ng)=1
97 knew(ng)=1
98 predictor_2d_step(ng)=.false.
99!
100 iic(ng)=0
101 nstp(ng)=1
102 nrhs(ng)=1
103 nnew(ng)=1
104!
105 synchro_flag(ng)=.true.
106 tdays(ng)=dstart
107 time(ng)=tdays(ng)*day2sec
108 ntstart(ng)=int((time(ng)-dstart*day2sec)/dt(ng))+1
109 ntend(ng)=ntimes(ng)
110 ntfirst(ng)=ntstart(ng)
111 END DO
112!
113!-----------------------------------------------------------------------
114! Clear tangent linear state variables. There is not need to clean
115! the basic state arrays since they were zeroth out at initialization
116! and bottom of previous iteration.
117!-----------------------------------------------------------------------
118!
119 DO ng=1,ngrids
120 DO tile=first_tile(ng),last_tile(ng),+1
121 CALL initialize_ocean (ng, tile, itlm)
122 END DO
123 END DO
124
125#ifdef SOLVE3D
126!
127!-----------------------------------------------------------------------
128! Compute basic state initial level thicknesses used for state norm
129! scaling. It uses zero time-averaged free-surface (rest state).
130! Therefore, the norm scaling is time invariant.
131!-----------------------------------------------------------------------
132!
133 DO ng=1,ngrids
134 DO tile=last_tile(ng),first_tile(ng),-1
135 CALL set_depth (ng, tile, itlm)
136 END DO
137 END DO
138#endif
139!
140!-----------------------------------------------------------------------
141! Unpack tangent linear initial conditions from state vector.
142!-----------------------------------------------------------------------
143!
144 DO ng=1,ngrids
145 DO tile=first_tile(ng),last_tile(ng),+1
146 CALL tl_unpack (ng, tile, nstr(ng), nend(ng), &
147 & state(ng)%vector)
148 END DO
149 END DO
150!
151!-----------------------------------------------------------------------
152! Compute initial tangent linear state dot product norm.
153!-----------------------------------------------------------------------
154!
155 DO ng=1,ngrids
156 DO tile=last_tile(ng),first_tile(ng),-1
157 CALL tl_statenorm (ng, tile, kstp(ng), nstp(ng), &
158 & statenorm(ng))
159 END DO
160 IF (master) THEN
161 WRITE (stdout,20) ' PROPAGATOR - Grid: ', ng, &
162 & ', Tangent Initial Norm: ', statenorm(ng)
163 END IF
164 END DO
165!
166!-----------------------------------------------------------------------
167! Read in initial forcing, climatology and assimilation data from
168! input NetCDF files. It loads the first relevant data record for
169! the time-interpolation between snapshots.
170!-----------------------------------------------------------------------
171!
172 DO ng=1,ngrids
173 CALL close_inp (ng, itlm)
174 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
175
176 CALL tl_get_idata (ng)
177 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
178
179 CALL tl_get_data (ng)
180 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
181 END DO
182!
183!-----------------------------------------------------------------------
184! Time-step the tangent linear model.
185!-----------------------------------------------------------------------
186!
187 DO ng=1,ngrids
188 IF (master) THEN
189 WRITE (stdout,30) 'TL', ng, ntstart(ng), ntend(ng)
190 END IF
191 time(ng)=time(ng)-dt(ng)
192 iic(ng)=ntstart(ng)-1
193 END DO
194
195#ifdef SOLVE3D
196 CALL tl_main3d (runinterval)
197#else
198 CALL tl_main2d (runinterval)
199#endif
200 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
201!
202!-----------------------------------------------------------------------
203! Clear nonlinear (basic state) and adjoint state variables.
204!-----------------------------------------------------------------------
205!
206 DO ng=1,ngrids
207 DO tile=first_tile(ng),last_tile(ng),+1
208 CALL initialize_ocean (ng, tile, inlm)
209 CALL initialize_ocean (ng, tile, iadm)
210#ifdef SOLVE3D
211 CALL initialize_coupling (ng, tile, 0)
212#endif
213 END DO
214 END DO
215
216#ifdef SOLVE3D
217!
218!-----------------------------------------------------------------------
219! Compute basic state final level thicknesses used for state norm
220! scaling. It uses zero time-averaged free-surface (rest state).
221! Therefore, the norm scaling is time invariant.
222!-----------------------------------------------------------------------
223!
224 DO ng=1,ngrids
225 DO tile=last_tile(ng),first_tile(ng),-1
226 CALL set_depth (ng, tile, itlm)
227 END DO
228 END DO
229#endif
230!
231!-----------------------------------------------------------------------
232! Compute final tangent linear state dot product norm.
233!-----------------------------------------------------------------------
234!
235 DO ng=1,ngrids
236 DO tile=first_tile(ng),last_tile(ng),+1
237 CALL tl_statenorm (ng, tile, knew(ng), nstp(ng), &
238 & statenorm(ng))
239 END DO
240 IF (master) THEN
241 WRITE (stdout,20) ' PROPAGATOR - Grid: ', ng, &
242 & ', Tangent Final Norm: ', statenorm(ng)
243 END IF
244 END DO
245!
246!=======================================================================
247! Backward integration with the adjoint model.
248!=======================================================================
249!
250! Initialize time stepping indices and counters.
251!
252 DO ng=1,ngrids
253 iif(ng)=1
254 indx1(ng)=1
255 ktmp=knew(ng)
256 kstp(ng)=1
257 krhs(ng)=3
258 knew(ng)=2
259 predictor_2d_step(ng)=.false.
260!
261 iic(ng)=0
262 ntmp=nstp(ng)
263 nstp(ng)=1
264 nrhs(ng)=1
265 nnew(ng)=2
266!
267 synchro_flag(ng)=.true.
268 tdays(ng)=dstart+dt(ng)*real(ntimes(ng),r8)*sec2day
269 time(ng)=tdays(ng)*day2sec
270 ntstart(ng)=ntimes(ng)+1
271 ntend(ng)=1
272 ntfirst(ng)=ntend(ng)
273 END DO
274!
275!-----------------------------------------------------------------------
276! Initialize adjoint model with the final tangent linear solution
277! scaled by the energy norm.
278!-----------------------------------------------------------------------
279!
280 DO ng=1,ngrids
281 DO tile=last_tile(ng),first_tile(ng),-1
282 CALL ad_ini_perturb (ng, tile, &
283 & ktmp, knew(ng), ntmp, nstp(ng))
284 END DO
285 END DO
286!
287!-----------------------------------------------------------------------
288! Read in initial forcing, climatology and assimilation data from
289! input NetCDF files. It loads the first relevant data record for
290! the time-interpolation between snapshots.
291!-----------------------------------------------------------------------
292!
293 DO ng=1,ngrids
294 CALL close_inp (ng, iadm)
295 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
296
297 CALL ad_get_idata (ng)
298 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
299
300 CALL ad_get_data (ng)
301 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
302 END DO
303!
304!-----------------------------------------------------------------------
305! Time-step the adjoint model backwards.
306!-----------------------------------------------------------------------
307!
308 DO ng=1,ngrids
309 IF (master) THEN
310 WRITE (stdout,30) 'AD', ng, ntstart(ng), ntend(ng)
311 END IF
312 time(ng)=time(ng)+dt(ng)
313 iic(ng)=ntstart(ng)+1
314 END DO
315
316#ifdef SOLVE3D
317 CALL ad_main3d (runinterval)
318#else
319 CALL ad_main2d (runinterval)
320#endif
321 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
322!
323!-----------------------------------------------------------------------
324! Clear nonlinear state (basic state) variables for next iteration
325! and to insure a rest state time averaged free-surface before adjoint
326! state norm scaling.
327!-----------------------------------------------------------------------
328!
329 DO ng=1,ngrids
330 DO tile=first_tile(ng),last_tile(ng),+1
331 CALL initialize_ocean (ng, tile, inlm)
332#ifdef SOLVE3D
333 CALL initialize_coupling (ng, tile, 0)
334#endif
335 END DO
336 END DO
337
338#ifdef SOLVE3D
339!
340!-----------------------------------------------------------------------
341! Compute basic state initial level thicknesses used for state norm
342! scaling. It uses zero free-surface (rest state). Therefore, the
343! norm scaling is time invariant.
344!-----------------------------------------------------------------------
345!
346 DO ng=1,ngrids
347 DO tile=last_tile(ng),first_tile(ng),-1
348 CALL set_depth (ng, tile, iadm)
349 END DO
350 END DO
351#endif
352!
353!-----------------------------------------------------------------------
354! Pack final adjoint solution into adjoint state vector.
355!-----------------------------------------------------------------------
356!
357 DO ng=1,ngrids
358 DO tile=first_tile(ng),last_tile(ng),+1
359 CALL ad_pack (ng, tile, nstr(ng), nend(ng), &
360 & ad_state(ng)%vector)
361 END DO
362 END DO
363!
364 10 FORMAT (/,a,i2.2,a,i3.3,a,i3.3/)
365 20 FORMAT (/,a,i2.2,a,1p,e15.6,/)
366 30 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
367 & ' (Grid: ',i2.2,' TimeSteps: ',i8.8,' - ',i8.8,')')
368!
369 RETURN
370 END SUBROUTINE propagator_op
371
372 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)
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
integer, dimension(:), allocatable nconv
logical, dimension(:), allocatable synchro_flag
logical, dimension(:), allocatable predictor_2d_step
real(dp), dimension(:), allocatable tdays
real(dp) dstart
real(dp), parameter sec2day
integer, dimension(:), allocatable ntend
integer exit_flag
integer, dimension(:), allocatable indx1
integer, dimension(:), allocatable ntfirst
real(dp), dimension(:), allocatable time
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_op(runinterval, 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 ad_pack(ng, tile, mstr, mend, ad_state)
Definition packing.F:341
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