ROMS
Loading...
Searching...
No Matches
ana_perturb.h
Go to the documentation of this file.
1!!
2 SUBROUTINE ana_perturb (ng, tile, model)
3!
4!! git $Id$
5!!======================================================================
6!! Copyright (c) 2002-2025 The ROMS Group !
7!! Licensed under a MIT/X style license !
8!! See License_ROMS.md !
9!=======================================================================
10! !
11! This routine perturbs initial conditions for momentum and tracers !
12! type variables using analytical expressions. !
13! !
14! It is also used to perturb the tangent linear and adjoint models !
15! at specified state variable and spatial (i,j,k) point to verify !
16! the correctness of these algorithms. This is activated with the !
17! SANITY_CHECK CPP switch. !
18! !
19! If each interior point is perturbed at one time, the resulting !
20! tangent linear (T) and adjoint (A) M-by-N matrices yield: !
21! !
22! T - tranpose(A) = 0 within round off !
23! !
24! That is, their inner product give a symmetric matrix. Here, M is !
25! the number of state points and N is the number of perturbations. !
26! In realistic applications, it is awkward to perturb all interior !
27! points for each state variable. Alternatively, random check at a !
28! specified points is inexpensive. The standard input "User" array !
29! is used to specify such point: !
30! !
31! INT(user(1)) => state variable to perturb !
32! INT(user(2)) => I-index to perturb !
33! INT(user(3)) => J-index to perturb !
34! INT(user(4)) => K-index to perturb (3D state fields) !
35! !
36!=======================================================================
37!
38 USE mod_param
39 USE mod_ncparam
40#ifdef ADJUST_BOUNDARY
41 USE mod_boundary
42#endif
43 USE mod_ocean
44#if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
45 USE mod_forces
46#endif
47 USE mod_stepping
48!
49! Imported variable declarations.
50!
51 integer, intent(in) :: ng, tile, model
52!
53! Local variable declarations.
54!
55 character (len=*), parameter :: MyFile = &
56 & __FILE__
57!
58#include "tile.h"
59!
60 CALL ana_perturb_tile (ng, tile, model, &
61 & lbi, ubi, lbj, ubj, lbij, ubij, &
62 & imins, imaxs, jmins, jmaxs, &
63 & kstp(ng), krhs(ng), knew(ng), &
64#ifdef SOLVE3D
65 & nstp(ng), nrhs(ng), nnew(ng), &
66#endif
67#ifdef ADJUST_BOUNDARY
68# ifdef SOLVE3D
69 & boundary(ng) % ad_t_obc, &
70 & boundary(ng) % ad_u_obc, &
71 & boundary(ng) % ad_v_obc, &
72# endif
73 & boundary(ng) % ad_ubar_obc, &
74 & boundary(ng) % ad_vbar_obc, &
75 & boundary(ng) % ad_zeta_obc, &
76#endif
77#ifdef ADJUST_WSTRESS
78 & forces(ng) % ad_ustr, &
79 & forces(ng) % ad_vstr, &
80#endif
81#if defined ADJUST_STFLUX && defined SOLVE3D
82 & forces(ng) % ad_tflux, &
83#endif
84#ifdef SOLVE3D
85 & ocean(ng) % ad_t, &
86 & ocean(ng) % ad_u, &
87 & ocean(ng) % ad_v, &
88#endif
89 & ocean(ng) % ad_ubar, &
90 & ocean(ng) % ad_vbar, &
91 & ocean(ng) % ad_zeta, &
92#ifdef ADJUST_BOUNDARY
93# ifdef SOLVE3D
94 & boundary(ng) % tl_t_obc, &
95 & boundary(ng) % tl_u_obc, &
96 & boundary(ng) % tl_v_obc, &
97# endif
98 & boundary(ng) % tl_ubar_obc, &
99 & boundary(ng) % tl_vbar_obc, &
100 & boundary(ng) % tl_zeta_obc, &
101#endif
102#ifdef ADJUST_WSTRESS
103 & forces(ng) % tl_ustr, &
104 & forces(ng) % tl_vstr, &
105#endif
106#if defined ADJUST_STFLUX && defined SOLVE3D
107 & forces(ng) % tl_tflux, &
108#endif
109#ifdef SOLVE3D
110 & ocean(ng) % tl_t, &
111 & ocean(ng) % tl_u, &
112 & ocean(ng) % tl_v, &
113#endif
114 & ocean(ng) % tl_ubar, &
115 & ocean(ng) % tl_vbar, &
116 & ocean(ng) % tl_zeta)
117!
118! Set analytical header file name used.
119!
120#ifdef DISTRIBUTE
121 IF (lanafile) THEN
122#else
123 IF (lanafile.and.(tile.eq.0)) THEN
124#endif
125 ananame(19)=myfile
126 END IF
127!
128 RETURN
129 END SUBROUTINE ana_perturb
130!
131!***********************************************************************
132 SUBROUTINE ana_perturb_tile (ng, tile, model, &
133 & LBi, UBi, LBj, UBj, LBij, UBij, &
134 & IminS, ImaxS, JminS, JmaxS, &
135 & kstp, krhs, knew, &
136#ifdef SOLVE3D
137 & nstp, nrhs, nnew, &
138#endif
139#ifdef ADJUST_BOUNDARY
140# ifdef SOLVE3D
141 & ad_t_obc, ad_u_obc, ad_v_obc, &
142# endif
143 & ad_ubar_obc, ad_vbar_obc, &
144 & ad_zeta_obc, &
145#endif
146#ifdef ADJUST_WSTRESS
147 & ad_ustr, ad_vstr, &
148#endif
149#if defined ADJUST_STFLUX && defined SOLVE3D
150 & ad_tflux, &
151#endif
152#ifdef SOLVE3D
153 & ad_t, ad_u, ad_v, &
154#endif
155 & ad_ubar, ad_vbar, ad_zeta, &
156
157
158#ifdef ADJUST_BOUNDARY
159# ifdef SOLVE3D
160 & tl_t_obc, tl_u_obc, tl_v_obc, &
161# endif
162 & tl_ubar_obc, tl_vbar_obc, &
163 & tl_zeta_obc, &
164#endif
165#ifdef ADJUST_WSTRESS
166 & tl_ustr, tl_vstr, &
167#endif
168#if defined ADJUST_STFLUX && defined SOLVE3D
169 & tl_tflux, &
170#endif
171#ifdef SOLVE3D
172 & tl_t, tl_u, tl_v, &
173#endif
174 & tl_ubar, tl_vbar, tl_zeta)
175!***********************************************************************
176!
177 USE mod_param
178 USE mod_parallel
179 USE mod_iounits
180 USE mod_ncparam
181 USE mod_scalars
182!
183! Imported variable declarations.
184!
185 integer, intent(in) :: ng, tile, model
186 integer, intent(in) :: LBi, UBi, LBj, UBj, LBij, UBij
187 integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
188 integer, intent(in) :: kstp, krhs, knew
189#ifdef SOLVE3D
190 integer, intent(in) :: nstp, nrhs, nnew
191#endif
192!
193#ifdef ASSUMED_SHAPE
194# ifdef ADJUST_BOUNDARY
195# ifdef SOLVE3D
196 real(r8), intent(inout) :: ad_t_obc(LBij:,:,:,:,:,:)
197 real(r8), intent(inout) :: ad_u_obc(LBij:,:,:,:,:)
198 real(r8), intent(inout) :: ad_v_obc(LBij:,:,:,:,:)
199# endif
200 real(r8), intent(inout) :: ad_ubar_obc(LBij:,:,:,:)
201 real(r8), intent(inout) :: ad_vbar_obc(LBij:,:,:,:)
202 real(r8), intent(inout) :: ad_zeta_obc(LBij:,:,:,:)
203# endif
204# ifdef ADJUST_WSTRESS
205 real(r8), intent(inout) :: ad_ustr(LBi:,LBj:,:,:)
206 real(r8), intent(inout) :: ad_vstr(LBi:,LBj:,:,:)
207# endif
208# if defined ADJUST_STFLUX && defined SOLVE3D
209 real(r8), intent(inout) :: ad_tflux(LBi:,LBj:,:,:,:)
210# endif
211# ifdef SOLVE3D
212 real(r8), intent(inout) :: ad_t(LBi:,LBj:,:,:,:)
213 real(r8), intent(inout) :: ad_u(LBi:,LBj:,:,:)
214 real(r8), intent(inout) :: ad_v(LBi:,LBj:,:,:)
215# endif
216 real(r8), intent(inout) :: ad_ubar(LBi:,LBj:,:)
217 real(r8), intent(inout) :: ad_vbar(LBi:,LBj:,:)
218 real(r8), intent(inout) :: ad_zeta(LBi:,LBj:,:)
219# ifdef ADJUST_BOUNDARY
220# ifdef SOLVE3D
221 real(r8), intent(inout) :: tl_t_obc(LBij:,:,:,:,:,:)
222 real(r8), intent(inout) :: tl_u_obc(LBij:,:,:,:,:)
223 real(r8), intent(inout) :: tl_v_obc(LBij:,:,:,:,:)
224# endif
225 real(r8), intent(inout) :: tl_ubar_obc(LBij:,:,:,:)
226 real(r8), intent(inout) :: tl_vbar_obc(LBij:,:,:,:)
227 real(r8), intent(inout) :: tl_zeta_obc(LBij:,:,:,:)
228# endif
229# ifdef ADJUST_WSTRESS
230 real(r8), intent(inout) :: tl_ustr(LBi:,LBj:,:,:)
231 real(r8), intent(inout) :: tl_vstr(LBi:,LBj:,:,:)
232# endif
233# if defined ADJUST_STFLUX && defined SOLVE3D
234 real(r8), intent(inout) :: tl_tflux(LBi:,LBj:,:,:,:)
235# endif
236# ifdef SOLVE3D
237 real(r8), intent(inout) :: tl_t(LBi:,LBj:,:,:,:)
238 real(r8), intent(inout) :: tl_u(LBi:,LBj:,:,:)
239 real(r8), intent(inout) :: tl_v(LBi:,LBj:,:,:)
240# endif
241 real(r8), intent(inout) :: tl_ubar(LBi:,LBj:,:)
242 real(r8), intent(inout) :: tl_vbar(LBi:,LBj:,:)
243 real(r8), intent(inout) :: tl_zeta(LBi:,LBj:,:)
244
245#else
246
247# ifdef ADJUST_BOUNDARY
248# ifdef SOLVE3D
249 real(r8), intent(inout) :: ad_t_obc(LBij:UBij,N(ng),4, &
250 & Nbrec(ng),2,NT(ng))
251 real(r8), intent(inout) :: ad_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
252 real(r8), intent(inout) :: ad_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
253# endif
254 real(r8), intent(inout) :: ad_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
255 real(r8), intent(inout) :: ad_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
256 real(r8), intent(inout) :: ad_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
257# endif
258# ifdef ADJUST_WSTRESS
259 real(r8), intent(inout) :: ad_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
260 real(r8), intent(inout) :: ad_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
261# endif
262# if defined ADJUST_STFLUX && defined SOLVE3D
263 real(r8), intent(inout) :: ad_tflux(LBi:UBi,LBj:UBj, &
264 & Nfrec(ng),2,NT(ng))
265# endif
266# ifdef SOLVE3D
267 real(r8), intent(inout) :: ad_t(LBi:UBI,LBj:UBj,N(ng),2,NT(ng))
268 real(r8), intent(inout) :: ad_u(LBi:UBi,LBj:UBj,N(ng),2)
269 real(r8), intent(inout) :: ad_v(LBi:UBi,LBj:UBj,N(ng),2)
270# endif
271 real(r8), intent(inout) :: ad_ubar(LBi:UBi,LBj:UBj,:)
272 real(r8), intent(inout) :: ad_vbar(LBi:UBi,LBj:UBj,:)
273 real(r8), intent(inout) :: ad_zeta(LBi:UBi,LBj:UBj,:)
274# ifdef ADJUST_BOUNDARY
275# ifdef SOLVE3D
276 real(r8), intent(inout) :: tl_t_obc(LBij:UBij,N(ng),4, &
277 & Nbrec(ng),2,NT(ng))
278 real(r8), intent(inout) :: tl_u_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
279 real(r8), intent(inout) :: tl_v_obc(LBij:UBij,N(ng),4,Nbrec(ng),2)
280# endif
281 real(r8), intent(inout) :: tl_ubar_obc(LBij:UBij,4,Nbrec(ng),2)
282 real(r8), intent(inout) :: tl_vbar_obc(LBij:UBij,4,Nbrec(ng),2)
283 real(r8), intent(inout) :: tl_zeta_obc(LBij:UBij,4,Nbrec(ng),2)
284# endif
285# ifdef ADJUST_WSTRESS
286 real(r8), intent(inout) :: tl_ustr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
287 real(r8), intent(inout) :: tl_vstr(LBi:UBi,LBj:UBj,Nfrec(ng),2)
288# endif
289# if defined ADJUST_STFLUX && defined SOLVE3D
290 real(r8), intent(inout) :: tl_tflux(LBi:UBi,LBj:UBj, &
291 & Nfrec(ng),2,NT(ng))
292# endif
293# ifdef SOLVE3D
294 real(r8), intent(inout) :: tl_u(LBi:UBi,LBj:UBj,N(ng),2)
295 real(r8), intent(inout) :: tl_v(LBi:UBi,LBj:UBj,N(ng),2)
296 real(r8), intent(inout) :: tl_t(LBi:UBi,LBj:UBj,N(ng),3,NT(ng))
297# endif
298 real(r8), intent(inout) :: tl_ubar(LBi:UBi,LBj:UBj,:)
299 real(r8), intent(inout) :: tl_vbar(LBi:UBi,LBj:UBj,:)
300 real(r8), intent(inout) :: tl_zeta(LBi:UBi,LBj:UBj,:)
301#endif
302!
303! Local variable declarations.
304!
305#ifdef ADJUST_BOUNDARY
306 logical :: Lperturb(4)
307!
308#endif
309 integer :: IperAD, JperAD, KperAD, ivarAD
310 integer :: IperTL, JperTL, KperTL, ivarTL
311 integer :: i, ib, ir, itrc, j, k
312
313#include "set_bounds.h"
314!
315!-----------------------------------------------------------------------
316! Set tangent and adjoint variable and random point to perturb.
317!-----------------------------------------------------------------------
318!
319 ivartl=int(user(1))
320 ivarad=int(user(2))
321 ipertl=int(user(3))
322 iperad=int(user(4))
323 jpertl=int(user(5))
324 jperad=int(user(6))
325#ifdef SOLVE3D
326 kpertl=int(user(7))
327 kperad=int(user(8))
328#endif
329 IF (master) THEN
330 IF (tlmodel) THEN
331 IF (ivartl.eq.isubar) THEN
332 WRITE (stdout,10) 'tl_ubar perturbed at (i,j) = ', &
333 & ipertl, jpertl
334 ELSE IF (ivartl.eq.isvbar) THEN
335 WRITE (stdout,10) 'tl_vbar perturbed at (i,j) = ', &
336 & ipertl, jpertl
337 ELSE IF (ivartl.eq.isfsur) THEN
338 WRITE (stdout,10) 'tl_zeta perturbed at (i,j) = ', &
339 & ipertl, jpertl
340#ifdef ADJUST_WSTRESS
341 ELSE IF (ivartl.eq.isustr) THEN
342 WRITE (stdout,10) 'tl_ustr perturbed at (i,j) = ', &
343 & ipertl, jpertl
344 ELSE IF (ivartl.eq.isvstr) THEN
345 WRITE (stdout,10) 'tl_vstr perturbed at (i,j) = ', &
346 & ipertl, jpertl
347#endif
348#ifdef SOLVE3D
349 ELSE IF (ivartl.eq.isuvel) THEN
350 WRITE (stdout,20) 'tl_u perturbed at (i,j,k) = ', &
351 & ipertl, jpertl, kpertl
352 ELSE IF (ivartl.eq.isvvel) THEN
353 WRITE (stdout,20) 'tl_v perturbed at (i,j,k) = ', &
354 & ipertl, jpertl, kpertl
355#endif
356 END IF
357#ifdef SOLVE3D
358 DO itrc=1,nt(ng)
359 IF (ivartl.eq.istvar(itrc)) THEN
360 WRITE (stdout,30) 'tl_t perturbed at (i,j,k,itrc) = ', &
361 & ipertl, jpertl, kpertl, itrc
362# ifdef ADJUST_STFLUX
363 ELSE IF (ivartl.eq.istsur(itrc)) THEN
364 WRITE (stdout,20) 'tl_tflux perturbed at (i,j,k,itrc) = ',&
365 & ipertl, jpertl, kpertl, itrc
366# endif
367 END IF
368 END DO
369#endif
370#ifdef ADJUST_BOUNDARY
371 IF (ivartl.eq.isubar) THEN
372 WRITE (stdout,10) 'tl_ubar_obc (S/N) perturbed at (i) = ', &
373 & ipertl
374 WRITE (stdout,10) 'tl_ubar_obc (E/W) perturbed at (j) = ', &
375 & jpertl
376 ELSE IF (ivartl.eq.isvbar) THEN
377 WRITE (stdout,10) 'tl_vbar_obc (S/N) perturbed at (i) = ', &
378 & ipertl
379 WRITE (stdout,10) 'tl_vbar_obc (E/W) perturbed at (j) = ', &
380 & jpertl
381 ELSE IF (ivartl.eq.isfsur) THEN
382 WRITE (stdout,10) 'tl_zeta_obc (S/N) perturbed at (i) = ', &
383 & ipertl
384 WRITE (stdout,10) 'tl_zeta_obc (E/W) perturbed at (j) = ', &
385 & jpertl
386# ifdef SOLVE3D
387 ELSE IF (ivartl.eq.isuvel) THEN
388 WRITE (stdout,10) 'tl_u_obc (S/N) perturbed at (i,k) = ', &
389 & ipertl, kpertl
390 WRITE (stdout,10) 'tl_u_obc (E/W) perturbed at (j,k) = ', &
391 & jpertl, kpertl
392 ELSE IF (ivartl.eq.isuvel) THEN
393 WRITE (stdout,10) 'tl_u_obc (S/N) perturbed at (i,k) = ', &
394 & ipertl, kpertl
395 WRITE (stdout,10) 'tl_u_obc (E/W) perturbed at (j,k) = ', &
396 & jpertl, kpertl
397# endif
398 END IF
399# ifdef SOLVE3D
400 DO itrc=1,nt(ng)
401 IF (ivartl.eq.istvar(itrc)) THEN
402 WRITE (stdout,20) 'tl_t_obc perturbed at (i,k,itrc) = ', &
403 & ipertl, kpertl, itrc
404 WRITE (stdout,20) 'tl_t_obc perturbed at (j,k,itrc) = ', &
405 & jpertl, kpertl, itrc
406 END IF
407 END DO
408# endif
409#endif
410 END IF
411 IF (admodel) THEN
412 IF (ivarad.eq.isubar) THEN
413 WRITE (stdout,40) 'ad_ubar perturbed at (i,j) = ', &
414 & iperad, jperad
415 ELSE IF (ivarad.eq.isvbar) THEN
416 WRITE (stdout,40) 'ad_vbar perturbed at (i,j) = ', &
417 & iperad, jperad
418 ELSE IF (ivarad.eq.isfsur) THEN
419 WRITE (stdout,40) 'ad_zeta perturbed at (i,j) = ', &
420 & iperad, jperad
421#ifdef ADJUST_WSTRESS
422 ELSE IF (ivarad.eq.isustr) THEN
423 WRITE (stdout,40) 'ad_ustr perturbed at (i,j,k) = ', &
424 & iperad, jperad
425 ELSE IF (ivarad.eq.isvstr) THEN
426 WRITE (stdout,40) 'ad_vstr perturbed at (i,j,k) = ', &
427 & iperad, jperad
428#endif
429#ifdef SOLVE3D
430 ELSE IF (ivarad.eq.isuvel) THEN
431 WRITE (stdout,50) 'ad_u perturbed at (i,j,k) = ', &
432 & iperad, jperad, kperad
433 ELSE IF (ivarad.eq.isvvel) THEN
434 WRITE (stdout,50) 'ad_v perturbed at (i,j,k) = ', &
435 & iperad, jperad, kperad
436#endif
437 END IF
438#ifdef SOLVE3D
439 DO itrc=1,nt(ng)
440 IF (ivarad.eq.istvar(itrc)) THEN
441 WRITE (stdout,60) 'ad_t perturbed at (i,j,k,itrc) = ', &
442 & iperad, jperad, kperad, itrc
443# ifdef ADJUST_STFLUX
444 ELSE IF (ivarad.eq.istsur(itrc)) THEN
445 WRITE (stdout,50) 'ad_tflux perturbed at (i,j,k,itrc) = ',&
446 & iperad, jperad, kperad, itrc
447# endif
448 END IF
449 END DO
450#endif
451#ifdef ADJUST_BOUNDARY
452 IF (ivarad.eq.isubar) THEN
453 WRITE (stdout,40) 'ad_ubar_obc (S/N) perturbed at (i) = ', &
454 & iperad
455 WRITE (stdout,40) 'ad_ubar_obc (E/W) perturbed at (j) = ', &
456 & jperad
457 ELSE IF (ivarad.eq.isvbar) THEN
458 WRITE (stdout,40) 'ad_vbar_obc (S/N) perturbed at (i) = ', &
459 & iperad
460 WRITE (stdout,40) 'ad_vbar_obc (E/W) perturbed at (j) = ', &
461 & jperad
462 ELSE IF (ivarad.eq.isfsur) THEN
463 WRITE (stdout,40) 'ad_zeta_obc (S/N) perturbed at (i) = ', &
464 & iperad
465 WRITE (stdout,40) 'ad_zeta_obc (E/W) perturbed at (j) = ', &
466 & jperad
467# ifdef SOLVE3D
468 ELSE IF (ivarad.eq.isuvel) THEN
469 WRITE (stdout,40) 'ad_u_obc (S/N) perturbed at (i,k) = ', &
470 & iperad, kperad
471 WRITE (stdout,40) 'ad_u_obc (E/W) perturbed at (j,k) = ', &
472 & jperad, kperad
473 ELSE IF (ivarad.eq.isuvel) THEN
474 WRITE (stdout,40) 'ad_u_obc (S/N) perturbed at (i,k) = ', &
475 & iperad, kperad
476 WRITE (stdout,40) 'ad_u_obc (E/W) perturbed at (j,k) = ', &
477 & jperad, kperad
478# endif
479 END IF
480# ifdef SOLVE3D
481 DO itrc=1,nt(ng)
482 IF (ivarad.eq.istvar(itrc)) THEN
483 WRITE (stdout,50) 'ad_t_obc perturbed at (i,k,itrc) = ', &
484 & iperad, kperad, itrc
485 WRITE (stdout,50) 'ad_t_obc perturbed at (j,k,itrc) = ', &
486 & jperad, kperad, itrc
487 END IF
488 END DO
489# endif
490#endif
491 END IF
492 END IF
493!
494!-----------------------------------------------------------------------
495! Peturb initial conditions for 2D momentum (m/s) components.
496!-----------------------------------------------------------------------
497!
498 IF (tlmodel) THEN
499 DO j=jstrt,jendt
500 DO i=istrp,iendt
501 IF ((ivartl.eq.isubar).and. &
502 & (i.eq.ipertl).and.(j.eq.jpertl)) THEN
503 tl_ubar(i,j,kstp)=1.0_r8
504 ELSE
505 tl_ubar(i,j,kstp)=0.0_r8
506 END IF
507 END DO
508 END DO
509 DO j=jstrp,jendt
510 DO i=istrt,iendt
511 IF ((ivartl.eq.isvbar).and. &
512 & (i.eq.ipertl).and.(j.eq.jpertl)) THEN
513 tl_vbar(i,j,kstp)=1.0_r8
514 ELSE
515 tl_vbar(i,j,kstp)=0.0_r8
516 END IF
517 END DO
518 END DO
519 END IF
520!
521 IF (admodel) THEN
522 DO j=jstrt,jendt
523 DO i=istrp,iendt
524 IF ((ivarad.eq.isubar).and. &
525 & (i.eq.iperad).and.(j.eq.jperad)) THEN
526 ad_ubar(i,j,knew)=1.0_r8
527 ELSE
528 ad_ubar(i,j,knew)=0.0_r8
529 END IF
530 END DO
531 END DO
532 DO j=jstrp,jendt
533 DO i=istrt,iendt
534 IF ((ivarad.eq.isvbar).and. &
535 & (i.eq.iperad).and.(j.eq.jperad)) THEN
536 ad_vbar(i,j,knew)=1.0_r8
537 ELSE
538 ad_vbar(i,j,knew)=0.0_r8
539 END IF
540 END DO
541 END DO
542 END IF
543#ifdef ADJUST_WSTRESS
544!
545!-----------------------------------------------------------------------
546! Peturb initial conditions for surface momentum stress (UNIT ????).
547!-----------------------------------------------------------------------
548!
549 IF (tlmodel) THEN
550 DO ir=1,nfrec(ng)
551 DO j=jstrt,jendt
552 DO i=istrp,iendt
553 IF ((ivartl.eq.isustr).and. &
554 & (i.eq.ipertl).and.(j.eq.jpertl).and. &
555 & (ir.eq.kpertl)) THEN
556 tl_ustr(i,j,ir,kstp)=1.0_r8
557 ELSE
558 tl_ustr(i,j,ir,kstp)=0.0_r8
559 END IF
560 END DO
561 END DO
562 DO j=jstrp,jendt
563 DO i=istrt,iendt
564 IF ((ivartl.eq.isvstr).and. &
565 & (i.eq.ipertl).and.(j.eq.jpertl).and. &
566 & (ir.eq.kpertl)) THEN
567 tl_vstr(i,j,ir,kstp)=1.0_r8
568 ELSE
569 tl_vstr(i,j,ir,kstp)=0.0_r8
570 END IF
571 END DO
572 END DO
573 END DO
574 END IF
575!
576 IF (admodel) THEN
577 DO ir=1,nfrec(ng)
578 DO j=jstrt,jendt
579 DO i=istrp,iendt
580 IF ((ivarad.eq.isustr).and. &
581 & (i.eq.iperad).and.(j.eq.jperad).and. &
582 & (ir.eq.kperad)) THEN
583 ad_ustr(i,j,ir,knew)=1.0_r8
584 ELSE
585 ad_ustr(i,j,ir,knew)=0.0_r8
586 END IF
587 END DO
588 END DO
589 DO j=jstrp,jendt
590 DO i=istrt,iendt
591 IF ((ivarad.eq.isvstr).and. &
592 & (i.eq.iperad).and.(j.eq.jperad).and. &
593 & (ir.eq.kperad)) THEN
594 ad_vstr(i,j,ir,knew)=1.0_r8
595 ELSE
596 ad_vstr(i,j,ir,knew)=0.0_r8
597 END IF
598 END DO
599 END DO
600 END DO
601 END IF
602#endif
603!
604!-----------------------------------------------------------------------
605! Perturb initial conditions for free-surface (m).
606!-----------------------------------------------------------------------
607!
608 IF (tlmodel) THEN
609 DO j=jstrt,jendt
610 DO i=istrt,iendt
611 IF ((ivartl.eq.isfsur).and. &
612 & (i.eq.ipertl).and.(j.eq.jpertl)) THEN
613 tl_zeta(i,j,kstp)=1.0_r8
614 ELSE
615 tl_zeta(i,j,kstp)=0.0_r8
616 END IF
617 END DO
618 END DO
619 END IF
620!
621 IF (admodel) THEN
622 DO j=jstrt,jendt
623 DO i=istrt,iendt
624 IF ((ivarad.eq.isfsur).and. &
625 & (i.eq.iperad).and.(j.eq.jperad)) THEN
626 ad_zeta(i,j,knew)=1.0_r8
627 ELSE
628 ad_zeta(i,j,knew)=0.0_r8
629 END IF
630 END DO
631 END DO
632 END IF
633
634#ifdef SOLVE3D
635!
636!-----------------------------------------------------------------------
637! Initial conditions for 3D momentum components (m/s).
638!-----------------------------------------------------------------------
639!
640 IF (tlmodel) THEN
641 DO k=1,n(ng)
642 DO j=jstrt,jendt
643 DO i=istrp,iendt
644 IF ((ivartl.eq.isuvel).and. &
645 & (i.eq.ipertl).and.(j.eq.jpertl).and. &
646 & (k.eq.kpertl)) THEN
647 tl_u(i,j,k,nstp)=1.0_r8
648 ELSE
649 tl_u(i,j,k,nstp)=0.0_r8
650 END IF
651 END DO
652 END DO
653 DO j=jstrp,jendt
654 DO i=istrt,iendt
655 IF ((ivartl.eq.isvvel).and. &
656 & (i.eq.ipertl).and.(j.eq.jpertl).and. &
657 & (k.eq.kpertl)) THEN
658 tl_v(i,j,k,nstp)=1.0_r8
659 ELSE
660 tl_v(i,j,k,nstp)=0.0_r8
661 END IF
662 END DO
663 END DO
664 END DO
665 END IF
666!
667 IF (admodel) THEN
668 DO k=1,n(ng)
669 DO j=jstrt,jendt
670 DO i=istrp,iendt
671 IF ((ivarad.eq.isuvel).and. &
672 & (i.eq.iperad).and.(j.eq.jperad).and. &
673 & (k.eq.kperad)) THEN
674 ad_u(i,j,k,nstp)=1.0_r8
675 ELSE
676 ad_u(i,j,k,nstp)=0.0_r8
677 END IF
678 END DO
679 END DO
680 DO j=jstrp,jendt
681 DO i=istrt,iendt
682 IF ((ivarad.eq.isvvel).and. &
683 & (i.eq.iperad).and.(j.eq.jperad).and. &
684 & (k.eq.kperad)) THEN
685 ad_v(i,j,k,nstp)=1.0_r8
686 ELSE
687 ad_v(i,j,k,nstp)=0.0_r8
688 END IF
689 END DO
690 END DO
691 END DO
692 END IF
693!
694!-----------------------------------------------------------------------
695! Perturb initial conditions for tracer type variables.
696!-----------------------------------------------------------------------
697!
698 IF (tlmodel) THEN
699 DO itrc=1,nt(ng)
700 DO k=1,n(ng)
701 DO j=jstrt,jendt
702 DO i=istrt,iendt
703 IF ((ivartl.eq.istvar(itrc)).and. &
704 & (i.eq.ipertl).and.(j.eq.jpertl).and. &
705 & (k.eq.kpertl)) THEN
706 tl_t(i,j,k,nstp,itrc)=1.0_r8
707 ELSE
708 tl_t(i,j,k,nstp,itrc)=0.0_r8
709 END IF
710 END DO
711 END DO
712 END DO
713 END DO
714 END IF
715!
716 IF (admodel) THEN
717 DO itrc=1,nt(ng)
718 DO k=1,n(ng)
719 DO j=jstrt,jendt
720 DO i=istrt,iendt
721 IF ((ivarad.eq.istvar(itrc)).and. &
722 & (i.eq.iperad).and.(j.eq.jperad).and. &
723 & (k.eq.kperad)) THEN
724 ad_t(i,j,k,nstp,itrc)=1.0_r8
725 ELSE
726 ad_t(i,j,k,nstp,itrc)=0.0_r8
727 END IF
728 END DO
729 END DO
730 END DO
731 END DO
732 END IF
733# ifdef ADJUST_STFLUX
734!
735!-----------------------------------------------------------------------
736! Perturb initial conditions for surface tracer flux.
737!-----------------------------------------------------------------------
738!
739 IF (tlmodel) THEN
740 DO itrc=1,nt(ng)
741 DO ir=1,nfrec(ng)
742 DO j=jstrt,jendt
743 DO i=istrt,iendt
744 IF ((ivartl.eq.istsur(itrc)).and. &
745 & (i.eq.ipertl).and.(j.eq.jpertl).and. &
746 & (ir.eq.kpertl)) THEN
747 tl_tflux(i,j,ir,nstp,itrc)=1.0_r8
748 ELSE
749 tl_tflux(i,j,ir,nstp,itrc)=0.0_r8
750 END IF
751 END DO
752 END DO
753 END DO
754 END DO
755 END IF
756!
757 IF (admodel) THEN
758 DO itrc=1,nt(ng)
759 DO ir=1,nfrec(ng)
760 DO j=jstrt,jendt
761 DO i=istrt,iendt
762 IF ((ivarad.eq.istsur(itrc)).and. &
763 & (i.eq.iperad).and.(j.eq.jperad).and. &
764 & (ir.eq.kperad)) THEN
765 ad_tflux(i,j,ir,nstp,itrc)=1.0_r8
766 ELSE
767 ad_tflux(i,j,ir,nstp,itrc)=0.0_r8
768 END IF
769 END DO
770 END DO
771 END DO
772 END DO
773 END IF
774# endif
775# ifdef ADJUST_BOUNDARY
776!
777!-----------------------------------------------------------------------
778! Perturb open boundary conditions.
779!-----------------------------------------------------------------------
780!
781 lperturb(iwest )=domain(ng)%Western_Edge (tile)
782 lperturb(ieast )=domain(ng)%Eastern_Edge (tile)
783 lperturb(isouth)=domain(ng)%Southern_Edge(tile)
784 lperturb(inorth)=domain(ng)%Northern_Edge(tile)
785
786 DO ir=1,nbrec(ng)
787 DO ib=1,4
788!
789! Perturb free-surface open boundaries.
790!
791 IF (lperturb(ib).and.lobc(ib,isfsur,ng)) THEN
792 IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
793 IF (tlmodel.and.(ivartl.eq.isfsur)) THEN
794 DO j=jstr,jend
795 IF (j.eq.jpertl) THEN
796 tl_zeta_obc(j,ib,ir,1)=1.0_r8
797 ELSE
798 tl_zeta_obc(j,ib,ir,1)=0.0_r8
799 END IF
800 END DO
801 ELSE IF (admodel.and.(ivarad.eq.isfsur)) THEN
802 DO j=jstr,jend
803 IF (j.eq.jperad) THEN
804 ad_zeta_obc(j,ib,ir,1)=1.0_r8
805 ELSE
806 ad_zeta_obc(j,ib,ir,1)=0.0_r8
807 END IF
808 END DO
809 END IF
810 ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
811 IF (tlmodel.and.(ivartl.eq.isfsur)) THEN
812 DO i=istr,iend
813 IF (i.eq.ipertl) THEN
814 tl_zeta_obc(i,ib,ir,1)=1.0_r8
815 ELSE
816 tl_zeta_obc(i,ib,ir,1)=0.0_r8
817 END IF
818 END DO
819 ELSE IF (admodel.and.(ivarad.eq.isfsur)) THEN
820 DO i=istr,iend
821 IF (i.eq.iperad) THEN
822 ad_zeta_obc(i,ib,ir,1)=1.0_r8
823 ELSE
824 ad_zeta_obc(i,ib,ir,1)=0.0_r8
825 END IF
826 END DO
827 END IF
828 END IF
829 END IF
830!
831! Perturb 2D U-momentum open boundaries.
832!
833 IF (lperturb(ib).and.lobc(ib,isubar,ng)) THEN
834 IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
835 IF (tlmodel.and.(ivartl.eq.isubar)) THEN
836 DO j=jstr,jend
837 IF (j.eq.jpertl) THEN
838 tl_ubar_obc(j,ib,ir,1)=1.0_r8
839 ELSE
840 tl_ubar_obc(j,ib,ir,1)=0.0_r8
841 END IF
842 END DO
843 ELSE IF (admodel.and.(ivarad.eq.isubar)) THEN
844 DO j=jstr,jend
845 IF (j.eq.jperad) THEN
846 ad_ubar_obc(j,ib,ir,1)=1.0_r8
847 ELSE
848 ad_ubar_obc(j,ib,ir,1)=0.0_r8
849 END IF
850 END DO
851 END IF
852 ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
853 IF (tlmodel.and.(ivartl.eq.isubar)) THEN
854 DO i=istru,iend
855 IF (i.eq.ipertl) THEN
856 tl_ubar_obc(i,ib,ir,1)=1.0_r8
857 ELSE
858 tl_ubar_obc(i,ib,ir,1)=0.0_r8
859 END IF
860 END DO
861 ELSE IF (admodel.and.(ivarad.eq.isubar)) THEN
862 DO i=istru,iend
863 IF (i.eq.iperad) THEN
864 ad_ubar_obc(i,ib,ir,1)=1.0_r8
865 ELSE
866 ad_ubar_obc(i,ib,ir,1)=0.0_r8
867 END IF
868 END DO
869 END IF
870 END IF
871 END IF
872!
873! Perturb 2D V-momentum open boundaries.
874!
875 IF (lperturb(ib).and.lobc(ib,isvbar,ng)) THEN
876 IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
877 IF (tlmodel.and.(ivartl.eq.isvbar)) THEN
878 DO j=jstrv,jend
879 IF (j.eq.jpertl) THEN
880 tl_vbar_obc(j,ib,ir,1)=1.0_r8
881 ELSE
882 tl_vbar_obc(j,ib,ir,1)=0.0_r8
883 END IF
884 END DO
885 ELSE IF (admodel.and.(ivarad.eq.isvbar)) THEN
886 DO j=jstrv,jend
887 IF (j.eq.jperad) THEN
888 ad_vbar_obc(j,ib,ir,1)=1.0_r8
889 ELSE
890 ad_vbar_obc(j,ib,ir,1)=0.0_r8
891 END IF
892 END DO
893 END IF
894 ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
895 IF (tlmodel.and.(ivartl.eq.isvbar)) THEN
896 DO i=istr,iend
897 IF (i.eq.ipertl) THEN
898 tl_vbar_obc(i,ib,ir,1)=1.0_r8
899 ELSE
900 tl_vbar_obc(i,ib,ir,1)=0.0_r8
901 END IF
902 END DO
903 ELSE IF (admodel.and.(ivarad.eq.isvbar)) THEN
904 DO i=istr,iend
905 IF (i.eq.iperad) THEN
906 ad_vbar_obc(i,ib,ir,1)=1.0_r8
907 ELSE
908 ad_vbar_obc(i,ib,ir,1)=0.0_r8
909 END IF
910 END DO
911 END IF
912 END IF
913 END IF
914
915# ifdef SOLVE3D
916!
917! Perturb 3D U-momentum open boundaries.
918!
919 IF (lperturb(ib).and.lobc(ib,isuvel,ng)) THEN
920 IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
921 IF (tlmodel.and.(ivartl.eq.isuvel)) THEN
922 DO k=1,n(ng)
923 DO j=jstr,jend
924 IF ((j.eq.jpertl).and.(k.eq.kpertl)) THEN
925 tl_u_obc(j,k,ib,ir,1)=1.0_r8
926 ELSE
927 tl_u_obc(j,k,ib,ir,1)=0.0_r8
928 END IF
929 END DO
930 END DO
931 ELSE IF (admodel.and.(ivarad.eq.isuvel)) THEN
932 DO k=1,n(ng)
933 DO j=jstr,jend
934 IF ((j.eq.jperad).and.(k.eq.kperad)) THEN
935 ad_u_obc(j,k,ib,ir,1)=1.0_r8
936 ELSE
937 ad_u_obc(j,k,ib,ir,1)=0.0_r8
938 END IF
939 END DO
940 END DO
941 END IF
942 ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
943 IF (tlmodel.and.(ivartl.eq.isuvel)) THEN
944 DO k=1,n(ng)
945 DO i=istru,iend
946 IF ((i.eq.ipertl).and.(k.eq.kpertl)) THEN
947 tl_u_obc(i,k,ib,ir,1)=1.0_r8
948 ELSE
949 tl_u_obc(i,k,ib,ir,1)=0.0_r8
950 END IF
951 END DO
952 END DO
953 ELSE IF (admodel.and.(ivarad.eq.isuvel)) THEN
954 DO k=1,n(ng)
955 DO i=istru,iend
956 IF ((i.eq.iperad).and.(k.eq.kperad)) THEN
957 ad_u_obc(i,k,ib,ir,1)=1.0_r8
958 ELSE
959 ad_u_obc(i,k,ib,ir,1)=0.0_r8
960 END IF
961 END DO
962 END DO
963 END IF
964 END IF
965 END IF
966!
967! Perturb 3D V-momentum open boundaries.
968!
969 IF (lperturb(ib).and.lobc(ib,isvvel,ng)) THEN
970 IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
971 IF (tlmodel.and.(ivartl.eq.isvvel)) THEN
972 DO k=1,n(ng)
973 DO j=jstrv,jend
974 IF ((j.eq.jpertl).and.(k.eq.kpertl)) THEN
975 tl_v_obc(j,k,ib,ir,1)=1.0_r8
976 ELSE
977 tl_v_obc(j,k,ib,ir,1)=0.0_r8
978 END IF
979 END DO
980 END DO
981 ELSE IF (admodel.and.(ivarad.eq.isvvel)) THEN
982 DO k=1,n(ng)
983 DO j=jstrv,jend
984 IF ((j.eq.jperad).and.(k.eq.kperad)) THEN
985 ad_v_obc(j,k,ib,ir,1)=1.0_r8
986 ELSE
987 ad_v_obc(j,k,ib,ir,1)=0.0_r8
988 END IF
989 END DO
990 END DO
991 END IF
992 ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
993 IF (tlmodel.and.(ivartl.eq.isvvel)) THEN
994 DO k=1,n(ng)
995 DO i=istr,iend
996 IF ((i.eq.ipertl).and.(k.eq.kpertl)) THEN
997 tl_v_obc(i,k,ib,ir,1)=1.0_r8
998 ELSE
999 tl_v_obc(i,k,ib,ir,1)=0.0_r8
1000 END IF
1001 END DO
1002 END DO
1003 ELSE IF (admodel.and.(ivarad.eq.isvvel)) THEN
1004 DO k=1,n(ng)
1005 DO i=istr,iend
1006 IF ((i.eq.iperad).and.(k.eq.kperad)) THEN
1007 ad_v_obc(i,k,ib,ir,1)=1.0_r8
1008 ELSE
1009 ad_v_obc(i,k,ib,ir,1)=0.0_r8
1010 END IF
1011 END DO
1012 END DO
1013 END IF
1014 END IF
1015 END IF
1016!
1017! Perturb tracers open boundaries.
1018!
1019 DO itrc=1,nt(ng)
1020 IF (lperturb(ib).and.lobc(ib,istvar(itrc),ng)) THEN
1021 IF ((ib.eq.iwest).or.(ib.eq.ieast)) THEN
1022 IF (tlmodel.and.(ivartl.eq.istvar(itrc))) THEN
1023 DO k=1,n(ng)
1024 DO j=jstr,jend
1025 IF ((j.eq.jpertl).and.(k.eq.kpertl)) THEN
1026 tl_t_obc(j,k,ib,ir,1,itrc)=1.0_r8
1027 ELSE
1028 tl_t_obc(j,k,ib,ir,1,itrc)=0.0_r8
1029 END IF
1030 END DO
1031 END DO
1032 ELSE IF (admodel.and.(ivarad.eq.istvar(itrc))) THEN
1033 DO k=1,n(ng)
1034 DO j=jstr,jend
1035 IF ((j.eq.jperad).and.(k.eq.kperad)) THEN
1036 ad_t_obc(j,k,ib,ir,1,itrc)=1.0_r8
1037 ELSE
1038 ad_t_obc(j,k,ib,ir,1,itrc)=0.0_r8
1039 END IF
1040 END DO
1041 END DO
1042 END IF
1043 ELSE IF ((ib.eq.isouth).or.(ib.eq.inorth)) THEN
1044 IF (tlmodel.and.(ivartl.eq.istvar(itrc))) THEN
1045 DO k=1,n(ng)
1046 DO i=istr,iend
1047 IF ((i.eq.ipertl).and.(k.eq.kpertl)) THEN
1048 tl_t_obc(i,k,ib,ir,1,itrc)=1.0_r8
1049 ELSE
1050 tl_t_obc(i,k,ib,ir,1,itrc)=0.0_r8
1051 END IF
1052 END DO
1053 END DO
1054 ELSE IF (admodel.and.(ivarad.eq.istvar(itrc))) THEN
1055 DO k=1,n(ng)
1056 DO i=istr,iend
1057 IF ((i.eq.iperad).and.(k.eq.kperad)) THEN
1058 ad_t_obc(i,k,ib,ir,1,itrc)=1.0_r8
1059 ELSE
1060 ad_t_obc(i,k,ib,ir,1,itrc)=0.0_r8
1061 END IF
1062 END DO
1063 END DO
1064 END IF
1065 END IF
1066 END IF
1067 END DO
1068# endif
1069 END DO
1070 END DO
1071# endif
1072#endif
1073!
1074 10 FORMAT (/,' ANA_PERTURB - Tangent ', a, 2i4,/)
1075#ifdef SOLVE3D
1076 20 FORMAT (/,' ANA_PERTURB - Tangent ', a, 3i4,/)
1077 30 FORMAT (/,' ANA_PERTURB - Tangent ', a, 4i4,/)
1078#endif
1079 40 FORMAT (/,' ANA_PERTURB - Adjoint ', a, 2i4,/)
1080#ifdef SOLVE3D
1081 50 FORMAT (/,' ANA_PERTURB - Adjoint ', a, 3i4,/)
1082 60 FORMAT (/,' ANA_PERTURB - Adjoint ', a, 4i4,/)
1083#endif
1084!
1085 RETURN
1086 END SUBROUTINE ana_perturb_tile
subroutine ana_perturb(ng, tile, model)
Definition ana_perturb.h:3
subroutine ana_perturb_tile(ng, tile, model, lbi, ubi, lbj, ubj, lbij, ubij, imins, imaxs, jmins, jmaxs, kstp, krhs, knew, nstp, nrhs, nnew, ad_t_obc, ad_u_obc, ad_v_obc, ad_ubar_obc, ad_vbar_obc, ad_zeta_obc, ad_ustr, ad_vstr, ad_tflux, ad_t, ad_u, ad_v, ad_ubar, ad_vbar, ad_zeta, tl_t_obc, tl_u_obc, tl_v_obc, tl_ubar_obc, tl_vbar_obc, tl_zeta_obc, tl_ustr, tl_vstr, tl_tflux, tl_t, tl_u, tl_v, tl_ubar, tl_vbar, tl_zeta)
type(t_boundary), dimension(:), allocatable boundary
type(t_forces), dimension(:), allocatable forces
Definition mod_forces.F:554
integer stdout
integer isvvel
integer isvbar
integer isvstr
integer, dimension(:), allocatable istvar
integer isuvel
integer isfsur
logical lanafile
integer isustr
integer isubar
integer, dimension(:), allocatable istsur
character(len=256), dimension(39) ananame
type(t_ocean), dimension(:), allocatable ocean
Definition mod_ocean.F:351
logical master
type(t_domain), dimension(:), allocatable domain
Definition mod_param.F:329
logical tlmodel
logical, dimension(:,:,:), allocatable lobc
integer, parameter iwest
real(r8), dimension(:), allocatable user
logical admodel
integer, parameter isouth
integer, parameter ieast
integer, parameter inorth
integer, dimension(:), allocatable kstp
integer, dimension(:), allocatable knew
integer, dimension(:), allocatable nrhs
integer, dimension(:), allocatable nnew
integer, dimension(:), allocatable krhs
integer, dimension(:), allocatable nstp