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