1 | MODULE step_tam |
---|
2 | #ifdef key_tam |
---|
3 | !!====================================================================== |
---|
4 | !! *** MODULE step_tam *** |
---|
5 | !! Time-stepping : manager of the adjoint ocean time stepping |
---|
6 | !! Tangent and Adjoint module |
---|
7 | !!====================================================================== |
---|
8 | !! History of |
---|
9 | !! the direct: ! 91-03 () Original code |
---|
10 | !! ! 91-11 (G. Madec) |
---|
11 | !! ! 92-06 (M. Imbard) add a first output record |
---|
12 | !! ! 96-04 (G. Madec) introduction of dynspg |
---|
13 | !! ! 96-04 (M.A. Foujols) introduction of passive tracer |
---|
14 | !! 8.0 ! 97-06 (G. Madec) new architecture of call |
---|
15 | !! 8.2 ! 97-06 (G. Madec, M. Imbard, G. Roullet) free surface |
---|
16 | !! 8.2 ! 99-02 (G. Madec, N. Grima) hpg implicit |
---|
17 | !! 8.2 ! 00-07 (J-M Molines, M. Imbard) Open Bondary Conditions |
---|
18 | !! 9.0 ! 02-06 (G. Madec) free form, suppress macro-tasking |
---|
19 | !! " " ! 04-08 (C. Talandier) New trends organization |
---|
20 | !! " " ! 05-01 (C. Ethe) Add the KPP closure scheme |
---|
21 | !! " " ! 05-11 (V. Garnier) Surface pressure gradient organization |
---|
22 | !! " " ! 05-11 (G. Madec) Reorganisation of tra and dyn calls |
---|
23 | !! " " ! 06-01 (L. Debreu, C. Mazauric) Agrif implementation |
---|
24 | !! " " ! 06-07 (S. Masson) restart using iom |
---|
25 | !! " " " ! 07-04 (K. Mogensen, A. Weaver, M. Martin) Assimilation interface |
---|
26 | !! History of the TAM |
---|
27 | !! ! 08-06 (A. Vidard and A. Weaver) Tangent and Adjoint version of 9.0 |
---|
28 | !! ! 08-11 (A. Vidard) Nemo v3 update |
---|
29 | !! ! 06-09 (F. Vigilant) Modified to split NEMOVAR / NEMOTAM |
---|
30 | !!---------------------------------------------------------------------- |
---|
31 | !! stp_tam : OPA system time-stepping (tangent linear) |
---|
32 | !! stp_adj : OPA system time-stepping (adjoint) |
---|
33 | !!---------------------------------------------------------------------- |
---|
34 | |
---|
35 | USE par_kind , ONLY: & ! Precision variables |
---|
36 | & wp |
---|
37 | USE par_oce , ONLY: & ! Ocean space and time domain variables |
---|
38 | & jpi, & |
---|
39 | & jpj, & |
---|
40 | & jpk, & |
---|
41 | & jpiglo |
---|
42 | USE oce , ONLY: & ! ocean dynamics and tracers variables |
---|
43 | & tb, sb, tn, sn, ta, & |
---|
44 | & un, vn, sshn, sshb, & |
---|
45 | & sa, ub, vb, & |
---|
46 | & ln_dynhpg_imp |
---|
47 | USE zdfkpp , ONLY: & |
---|
48 | & lk_zdfkpp ! KPP vertical mixing flag |
---|
49 | USE dynspg_oce , ONLY: & |
---|
50 | & lk_dynspg_rl ! Rigid-lid flag |
---|
51 | USE dom_oce , ONLY: & ! ocean space and time domain variables |
---|
52 | & n_cla, e1u, e2u, & |
---|
53 | & e1v, e2v, e1t, e2t, & |
---|
54 | # if defined key_vvl |
---|
55 | & e3t_1, & |
---|
56 | # else |
---|
57 | # if defined key_zco |
---|
58 | & e3t_0, e3w_0, & |
---|
59 | # else |
---|
60 | & e3t, e3u, e3v, e3w, & |
---|
61 | # endif |
---|
62 | # endif |
---|
63 | & nldi, nldj, nlei, & |
---|
64 | & nlej, lk_vvl, ln_zps,& |
---|
65 | & narea, mig, mjg, & |
---|
66 | & umask, vmask, tmask |
---|
67 | USE in_out_manager, ONLY: & ! I/O manager |
---|
68 | & numout, lwp, & |
---|
69 | & nit000, nitend, & |
---|
70 | & nstop, ctl_stop |
---|
71 | USE daymod , ONLY: & |
---|
72 | & adatrj |
---|
73 | USE trabbc , ONLY: &! bottom boundary condition |
---|
74 | & lk_trabbc |
---|
75 | USE traqsr , ONLY: & |
---|
76 | & ln_traqsr |
---|
77 | ! solar radiation penetration flag |
---|
78 | USE oce_tam , ONLY: & ! Tangent linear and adjoint variables |
---|
79 | & oce_tam_init, & |
---|
80 | & un_tl, vn_tl, tn_tl, & |
---|
81 | & sn_tl, ua_tl, va_tl, & |
---|
82 | & ta_tl, sa_tl, ub_tl, & |
---|
83 | & vb_tl, tb_tl, sb_tl, & |
---|
84 | & un_ad, vn_ad, tn_ad, & |
---|
85 | & sn_ad, ua_ad, va_ad, & |
---|
86 | & ta_ad, sa_ad, ub_ad, & |
---|
87 | & vb_ad, tb_ad, sb_ad, & |
---|
88 | & sshn_tl, sshn_ad, & |
---|
89 | #if defined key_obc |
---|
90 | & sshb_ad, & |
---|
91 | #endif |
---|
92 | & rn2_tl, rn2_ad, & |
---|
93 | & rhd_tl, rhop_tl, & |
---|
94 | & rhd_ad, rhop_ad, & |
---|
95 | & gtu_tl, gsu_tl, & |
---|
96 | & gru_tl, gtv_tl, & |
---|
97 | & gsv_tl, grv_tl, & |
---|
98 | & gtu_ad, gsu_ad, & |
---|
99 | & gru_ad, gtv_ad, & |
---|
100 | & gsv_ad, grv_ad, & |
---|
101 | & ssha_tl, sshb_tl, & |
---|
102 | & rotb_tl, hdivb_tl, & |
---|
103 | & rotn_tl, hdivn_tl, & |
---|
104 | & wn_tl |
---|
105 | USE lbclnk_tam |
---|
106 | USE daymod_tam ! calendar (adjoint of day routine) |
---|
107 | USE sbc_oce_tam |
---|
108 | USE sbcmod_tam |
---|
109 | USE traqsr_tam ! solar radiation penetration (adjoint of tra_qsr routine) |
---|
110 | USE trasbc_tam ! surface boundary condition (adjoint of tra_sbc routine) |
---|
111 | USE trabbc_tam ! bottom boundary condition (adjoint of tra_bbc routine) |
---|
112 | USE tradmp_tam ! internal damping (adjoint of tra_dmp routine) |
---|
113 | USE traadv_tam ! advection scheme control (adjoint of tra_adv_ctl routine) |
---|
114 | USE traldf_tam ! lateral mixing (adjoint of tra_ldf routine) |
---|
115 | USE cla_tam ! cross land advection (adjoint of tra_cla routine) |
---|
116 | ! zdfkpp ! KPP non-local tracer fluxes (adjoint of tra_kpp routine) |
---|
117 | USE trazdf_tam ! vertical mixing (adjoint of tra_zdf routine) |
---|
118 | USE tranxt_tam ! time-stepping (adjoint of tra_nxt routine) |
---|
119 | USE eosbn2_tam ! equation of state (adjoint of eos_bn2 routine) |
---|
120 | |
---|
121 | USE dynadv_tam ! advection (adjoint of dyn_adv routine) |
---|
122 | USE dynvor_tam ! vorticity term (adjoint of dyn_vor routine) |
---|
123 | USE dynhpg_tam ! hydrostatic pressure grad. (adjoint of dyn_hpg routine) |
---|
124 | USE dynldf_tam ! lateral momentum diffusion (adjoint of dyn_ldf routine) |
---|
125 | USE dynzdf_tam ! vertical diffusion (adjoint of dyn_zdf routine) |
---|
126 | USE dynspg_tam ! surface pressure gradient (adjoint of dyn_spg routine) |
---|
127 | USE dynnxt_tam ! time-stepping (adjoint of dyn_nxt routine) |
---|
128 | |
---|
129 | ! USE bdy_par_tam |
---|
130 | ! USE bdydta_tam |
---|
131 | |
---|
132 | USE divcur_tam ! hor. divergence and curl (adjoint of div & cur routines) |
---|
133 | USE cla_div_tam ! cross land: hor. divergence (adjoint of div_cla routine) |
---|
134 | USE wzvmod_tam ! vertical velocity (adjoint of wzv routine) |
---|
135 | |
---|
136 | USE zdfkpp_tam ! KPP vertical mixing |
---|
137 | |
---|
138 | USE zpshde_tam ! partial step: hor. derivative (adjoint of zps_hde routine) |
---|
139 | |
---|
140 | !! USE diaobs_tam ! obs-minus-model (adjoint of assimilation) (adjoint of dia_obs routine) |
---|
141 | |
---|
142 | USE trj_tam |
---|
143 | |
---|
144 | USE stpctl_tam ! time stepping control (adjoint of stp_ctl routine) |
---|
145 | |
---|
146 | USE gridrandom, ONLY: & |
---|
147 | ! Random Gaussian noise on grids |
---|
148 | & grid_random, & |
---|
149 | & grid_rd_sd |
---|
150 | USE dotprodfld, ONLY: & ! Computes dot product for 3D and 2D fields |
---|
151 | & dot_product |
---|
152 | USE tstool_tam, ONLY: & |
---|
153 | & prntst_adj, & |
---|
154 | & prntst_tlm, & |
---|
155 | & stdemp, & |
---|
156 | & stdu,stdv, & |
---|
157 | & stdt, & |
---|
158 | & stds, stdssh, & |
---|
159 | & stdr |
---|
160 | |
---|
161 | USE paresp, ONLY: & |
---|
162 | ! Weights for an energy-type scalar product |
---|
163 | & wesp_t, & |
---|
164 | & wesp_s, & |
---|
165 | & wesp_u, & |
---|
166 | & wesp_ssh |
---|
167 | |
---|
168 | USE istate_tam !: Initial state setting (istate_init routine) |
---|
169 | USE sol_oce_tam |
---|
170 | USE trc_oce_tam |
---|
171 | |
---|
172 | USE step |
---|
173 | USE iom ! for dumping an abort state |
---|
174 | #if defined key_obc |
---|
175 | USE obc_par, ONLY : lk_obc |
---|
176 | USE obcdta ! open boundary condition data (obc_dta routine) |
---|
177 | #endif |
---|
178 | |
---|
179 | #if defined key_agrif |
---|
180 | #error 'agrif not yet implemented in nemotam' |
---|
181 | #endif |
---|
182 | |
---|
183 | IMPLICIT NONE |
---|
184 | PRIVATE |
---|
185 | |
---|
186 | PUBLIC stp_tan, & |
---|
187 | & stp_adj, & ! called by simvar.F90 |
---|
188 | & stp_tlm_tst, & |
---|
189 | & stp_adj_tst |
---|
190 | |
---|
191 | !! * Substitutions |
---|
192 | # include "domzgr_substitute.h90" |
---|
193 | # include "zdfddm_substitute.h90" |
---|
194 | |
---|
195 | CONTAINS |
---|
196 | #if defined key_agrif |
---|
197 | #error 'agrif not yet implemented in nemotam' |
---|
198 | ! SUBROUTINE stp_tan( ) |
---|
199 | #else |
---|
200 | SUBROUTINE stp_tan( kstp ) |
---|
201 | #endif |
---|
202 | !!---------------------------------------------------------------------- |
---|
203 | !! *** ROUTINE stp_tan *** |
---|
204 | !! |
---|
205 | !! ** Purpose of the direct routine: |
---|
206 | !! - Time stepping of OPA (momentum and active tracer eqs.) |
---|
207 | !! - Time stepping of LIM (dynamic and thermodynamic eqs.) |
---|
208 | !! - Tme stepping of TRC (passive tracer eqs.) |
---|
209 | !! |
---|
210 | !! ** Method of the direct routine: |
---|
211 | !! -1- Update forcings and data |
---|
212 | !! -2- Update ocean physics |
---|
213 | !! -3- Compute the t and s trends |
---|
214 | !! -4- Update t and s |
---|
215 | !! -5- Compute the momentum trends |
---|
216 | !! -6- Update the horizontal velocity |
---|
217 | !! -7- Compute the diagnostics variables (rd,N2, div,cur,w) |
---|
218 | !! -8- Outputs and diagnostics |
---|
219 | !!---------------------------------------------------------------------- |
---|
220 | !! * Arguments |
---|
221 | #if defined key_agrif |
---|
222 | INTEGER :: kstp ! ocean time-step index |
---|
223 | #else |
---|
224 | INTEGER, INTENT( in ) :: kstp ! ocean time-step index |
---|
225 | #endif |
---|
226 | |
---|
227 | !! * local declarations |
---|
228 | INTEGER :: indic ! error indicator if < 0 |
---|
229 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zta_tmp, zsa_tmp |
---|
230 | !! --------------------------------------------------------------------- |
---|
231 | |
---|
232 | #if defined key_agrif |
---|
233 | kstp = nit000 + Agrif_Nb_Step() |
---|
234 | ! IF ( Agrif_Root() .and. lwp) Write(*,*) '---' |
---|
235 | ! IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp |
---|
236 | #endif |
---|
237 | indic = 1 ! reset to no error condition |
---|
238 | |
---|
239 | CALL day_tam( kstp, 0 ) ! Calendar |
---|
240 | |
---|
241 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
242 | ! Update data, open boundaries, surface boundary condition (including sea-ice) |
---|
243 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
244 | |
---|
245 | ! update 3D temperature data ... not needed in tangent |
---|
246 | |
---|
247 | ! update 3D salinity data ... not needed in tangent (to be investigated, see sbc_ssr) |
---|
248 | |
---|
249 | CALL sbc_tan ( kstp ) ! Sea boundary condition (including sea-ice) |
---|
250 | |
---|
251 | ! IF( lk_obc ) CALL obc_dta( kstp ) ! ??? ! update dynamic and tracer data at open boundaries |
---|
252 | ! update dynamic and tracer data at open boundaries ... not needed for the time being, to investigate whenever we do obc. |
---|
253 | |
---|
254 | ! compute phase velocities at open boundaries ... not needed for the time being, to investigate whenever we do obc. |
---|
255 | ! IF( lk_bdy ) CALL bdy_dta_tan( kstp ) ! update dynamic and tracer data at unstructured open boundary |
---|
256 | |
---|
257 | ! Output the initial state and forcings ... not needed in tangent |
---|
258 | |
---|
259 | ! saving direct variables ua,va, ta, sa before entering in tracer |
---|
260 | zta_tmp (:,:,:) = ta (:,:,:) |
---|
261 | zsa_tmp (:,:,:) = sa (:,:,:) |
---|
262 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
263 | ! Ocean physics update |
---|
264 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
265 | !----------------------------------------------------------------------- |
---|
266 | ! VERTICAL PHYSICS |
---|
267 | !----------------------------------------------------------------------- |
---|
268 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
269 | !----------------------------------------------------------------------- |
---|
270 | |
---|
271 | CALL bn2_tan( tb, sb, tb_tl, sb_tl, rn2_tl ) ! before Brunt-Vaisala frequency |
---|
272 | |
---|
273 | ! ! Vertical eddy viscosity and diffusivity coefficients |
---|
274 | ! Richardson number dependent Kz ... not available |
---|
275 | ! TKE closure scheme for Kz ... not available |
---|
276 | ! KPP closure scheme for Kz ... not available |
---|
277 | |
---|
278 | ! Constant Kz (reset avt, avm[uv] to the background value)... not available |
---|
279 | |
---|
280 | ! lk_zdfcst_tan: Constant Kz read from the reference trajectory |
---|
281 | |
---|
282 | |
---|
283 | ! ! increase diffusivity at rivers mouths... not needed in tangent |
---|
284 | |
---|
285 | ! enhanced vertical eddy diffusivity ... not needed in tangent with lk_zdfcst_tan |
---|
286 | |
---|
287 | ! double diffusive mixing ... not needed in tangent with lk_zdfcst_tan |
---|
288 | |
---|
289 | ! bottom friction... not needed in tangent with lk_zdfcst_tan |
---|
290 | |
---|
291 | ! mixed layer depth... not needed in tangent with lk_zdfcst_tan |
---|
292 | |
---|
293 | |
---|
294 | !----------------------------------------------------------------------- |
---|
295 | ! LATERAL PHYSICS |
---|
296 | !----------------------------------------------------------------------- |
---|
297 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
298 | !----------------------------------------------------------------------- |
---|
299 | |
---|
300 | ! before slope of the lateral mixing... not needed in tangent with lk_zdfcst_tan |
---|
301 | |
---|
302 | #if defined key_traldf_c2d |
---|
303 | ! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_tan |
---|
304 | #endif |
---|
305 | |
---|
306 | |
---|
307 | #if defined key_top |
---|
308 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
309 | ! Passive Tracer Model |
---|
310 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
311 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
312 | !----------------------------------------------------------------------- |
---|
313 | |
---|
314 | ! time-stepping... not needed in tangent for the time being |
---|
315 | |
---|
316 | #endif |
---|
317 | |
---|
318 | ta (:,:,:) = zta_tmp (:,:,:) |
---|
319 | sa (:,:,:) = zsa_tmp (:,:,:) |
---|
320 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
321 | ! Active tracers |
---|
322 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
323 | ! N.B. ua, va arrays are used as workspace in this section |
---|
324 | !----------------------------------------------------------------------- |
---|
325 | |
---|
326 | ta_tl(:,:,:) = 0.e0 ! set tracer trends to zero |
---|
327 | sa_tl(:,:,:) = 0.e0 |
---|
328 | |
---|
329 | ! apply tracer assimilation increment ... not needed in tangent |
---|
330 | ! apply bias directly to tn, sn, tb, sb ... not needed in tangent |
---|
331 | CALL tra_sbc_tan( kstp ) ! surface boundary condition |
---|
332 | |
---|
333 | IF( ln_traqsr ) CALL tra_qsr_tan( kstp ) ! penetrative solar radiation qsr |
---|
334 | |
---|
335 | IF( lk_trabbc ) CALL tra_bbc_tan( kstp ) ! bottom heat flux |
---|
336 | |
---|
337 | ! diffusive bottom boundary layer scheme ... currently not available |
---|
338 | ! advective (and/or diffusive) bottom boundary layer scheme ... currently not available |
---|
339 | |
---|
340 | IF( lk_tradmp ) CALL tra_dmp_tan( kstp ) ! internal damping trends |
---|
341 | |
---|
342 | CALL tra_adv_tan( kstp ) ! horizontal & vertical advection |
---|
343 | |
---|
344 | IF( n_cla == 1 ) CALL tra_cla_tan( kstp ) ! Cross Land Advection (Update Hor. advection) |
---|
345 | |
---|
346 | ! IF( lk_zdfkpp ) CALL tra_kpp_tan( kstp ) ! KPP non-local tracer fluxes |
---|
347 | |
---|
348 | CALL tra_ldf_tan( kstp ) ! lateral mixing |
---|
349 | |
---|
350 | #if defined key_agrif |
---|
351 | ! tracers sponge ... not available |
---|
352 | #endif |
---|
353 | CALL tra_zdf_tan( kstp ) ! vertical mixing |
---|
354 | |
---|
355 | ! update the new (t,s) fields by non |
---|
356 | ! penetrative convective adjustment ... not available |
---|
357 | |
---|
358 | CALL tra_nxt_tan( kstp ) ! tracer fields at next time step |
---|
359 | |
---|
360 | IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg |
---|
361 | CALL eos_tan( ta, sa, ta_tl, sa_tl, rhd_tl, rhop_tl ) ! Time-filtered in situ density used in dynhpg module |
---|
362 | IF( ln_zps ) CALL zps_hde_tan( kstp, ta, sa, ta_tl, sa_tl, rhd_tl,& ! Partial steps: time filtered hor. gradient |
---|
363 | & gtu_tl, gsu_tl, gru_tl, & ! of t, s, rd at the bottom ocean level |
---|
364 | & gtv_tl, gsv_tl, grv_tl ) |
---|
365 | ELSE ! centered hpg (default case) |
---|
366 | CALL eos_tan( tb, sb, tb_tl, sb_tl, rhd_tl, rhop_tl ) ! now (swap=before) in situ density for dynhpg module |
---|
367 | IF( ln_zps ) CALL zps_hde_tan( kstp, tb, sb, tb_tl, sb_tl, rhd_tl,& ! Partial steps: now horizontal gradient |
---|
368 | & gtu_tl, gsu_tl, gru_tl, & ! of t, s, rd at the bottom ocean level |
---|
369 | & gtv_tl, gsv_tl, grv_tl ) |
---|
370 | ENDIF |
---|
371 | |
---|
372 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
373 | ! Dynamics |
---|
374 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
375 | ! N.B. ta, sa arrays are used as workspace in this section |
---|
376 | !----------------------------------------------------------------------- |
---|
377 | |
---|
378 | ua_tl(:,:,:) = 0.e0 ! set dynamics trends to zero |
---|
379 | va_tl(:,:,:) = 0.e0 |
---|
380 | |
---|
381 | ! apply dynamics assimilation increment ... not needed in tangent |
---|
382 | |
---|
383 | CALL dyn_adv_tan( kstp ) ! advection (vector or flux form) |
---|
384 | CALL dyn_vor_tan( kstp ) ! vorticity term including Coriolis |
---|
385 | CALL dyn_ldf_tan( kstp ) ! lateral mixing |
---|
386 | #if defined key_agrif |
---|
387 | ! momemtum sponge ... not available |
---|
388 | #endif |
---|
389 | |
---|
390 | CALL dyn_hpg_tan( kstp ) ! horizontal gradient of Hydrostatic pressure |
---|
391 | CALL dyn_zdf_tan( kstp ) ! vertical diffusion |
---|
392 | |
---|
393 | IF( lk_dynspg_rl ) THEN |
---|
394 | ! surface pressure gradient at open boundaries ... not available |
---|
395 | ENDIF |
---|
396 | indic=0 |
---|
397 | !i |
---|
398 | CALL dyn_spg_tan( kstp, indic ) ! surface pressure gradient |
---|
399 | CALL dyn_nxt_tan( kstp ) ! lateral velocity at next time step |
---|
400 | |
---|
401 | ! vertical mesh at next time step ... not available |
---|
402 | |
---|
403 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
404 | ! Computation of diagnostic variables |
---|
405 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
406 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
407 | !----------------------------------------------------------------------- |
---|
408 | |
---|
409 | CALL div_cur_tan( kstp ) ! Horizontal divergence & Relative vorticity |
---|
410 | |
---|
411 | IF( n_cla == 1 ) CALL div_cla_tan( kstp ) ! Cross Land Advection (Update Hor. divergence) |
---|
412 | |
---|
413 | CALL wzv_tan( kstp ) ! Vertical velocity |
---|
414 | |
---|
415 | CALL trj_rea( kstp, 1) ! ... Read basic state trajectory at end of current step |
---|
416 | |
---|
417 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
418 | ! Control, and restarts |
---|
419 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
420 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
421 | !----------------------------------------------------------------------- |
---|
422 | ! ! Time loop: control and print |
---|
423 | CALL stp_ctl_tan( kstp, indic, 0 ) |
---|
424 | IF( indic < 0 ) CALL ctl_stop( 'step_tan: indic < 0' ) |
---|
425 | |
---|
426 | ! close input ocean restart file ... not needed in tangent |
---|
427 | ! write output ocean restart file... not needed in tangent |
---|
428 | ! write open boundary restart file... not needed in tangent |
---|
429 | |
---|
430 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
431 | ! diagnostics and outputs |
---|
432 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
433 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
434 | !----------------------------------------------------------------------- |
---|
435 | |
---|
436 | IF ( nstop == 0 ) THEN ! Diagnostics |
---|
437 | ! drifting Floats... not needed in tangent |
---|
438 | ! trends: dynamics ... not needed in tangent |
---|
439 | ! trends: active tracers... not needed in tangent |
---|
440 | ! trends: Mixed-layer ... not needed in tangent |
---|
441 | ! trends: vorticity budget... not needed in tangent |
---|
442 | ! Surface pressure diagnostics... not needed in tangent |
---|
443 | ! Thermocline depth (20 degres isotherm depth)... not needed in tangent |
---|
444 | ! basin averaged diagnostics... not needed in tangent |
---|
445 | ! dynamical heigh diagnostics... not needed in tangent |
---|
446 | ! Fresh water budget diagnostics... not needed in tangent |
---|
447 | !!! IF( lk_diaobs ) CALL dia_obs_tan( kstp ) ! obs-minus-model (assimilation) diagnostics NOTE: may be better off outside this module |
---|
448 | ! Poleward TRansports diagnostics... not needed in tangent |
---|
449 | |
---|
450 | |
---|
451 | ! ! Outputs |
---|
452 | ! ocean model: outputs... not needed in tangent |
---|
453 | |
---|
454 | ENDIF |
---|
455 | |
---|
456 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
457 | ! Coupled mode |
---|
458 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
459 | |
---|
460 | ! coupled mode : field exchanges ... not available for the next 20 years |
---|
461 | ! |
---|
462 | END SUBROUTINE stp_tan |
---|
463 | |
---|
464 | #if defined key_agrif |
---|
465 | #error 'agrif not yet implemented in nemotam' |
---|
466 | ! SUBROUTINE stp_adj( ) |
---|
467 | #else |
---|
468 | SUBROUTINE stp_adj( kstp ) |
---|
469 | #endif |
---|
470 | !!---------------------------------------------------------------------- |
---|
471 | !! *** ROUTINE stp_adj *** |
---|
472 | !! |
---|
473 | !! ** Purpose of the direct routine: |
---|
474 | !! - Time stepping of OPA (momentum and active tracer eqs.) |
---|
475 | !! - Time stepping of LIM (dynamic and thermodynamic eqs.) |
---|
476 | !! - Tme stepping of TRC (passive tracer eqs.) |
---|
477 | !! |
---|
478 | !! ** Method of the direct routine: |
---|
479 | !! -1- Update forcings and data |
---|
480 | !! -2- Update ocean physics |
---|
481 | !! -3- Compute the t and s trends |
---|
482 | !! -4- Update t and s |
---|
483 | !! -5- Compute the momentum trends |
---|
484 | !! -6- Update the horizontal velocity |
---|
485 | !! -7- Compute the diagnostics variables (rd,N2, div,cur,w) |
---|
486 | !! -8- Outputs and diagnostics |
---|
487 | !!---------------------------------------------------------------------- |
---|
488 | !! * Arguments |
---|
489 | #if defined key_agrif |
---|
490 | INTEGER :: kstp ! ocean time-step index |
---|
491 | #else |
---|
492 | INTEGER, INTENT( in ) :: kstp ! ocean time-step index |
---|
493 | #endif |
---|
494 | |
---|
495 | REAL(wp), DIMENSION(jpi,jpj,jpk) :: zta_tmp, zsa_tmp |
---|
496 | |
---|
497 | !! * local declarations |
---|
498 | INTEGER :: indic ! error indicator if < 0 |
---|
499 | !! --------------------------------------------------------------------- |
---|
500 | |
---|
501 | #if defined key_agrif |
---|
502 | kstp = nit000 + Agrif_Nb_Step() |
---|
503 | ! IF ( Agrif_Root() .and. lwp) Write(*,*) '---' |
---|
504 | ! IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp |
---|
505 | #endif |
---|
506 | indic = 1 ! reset to no error condition |
---|
507 | |
---|
508 | CALL day_tam( kstp, 0 ) ! Calendar |
---|
509 | |
---|
510 | ! ... Read basic state trajectory at end of previous step |
---|
511 | CALL trj_rea( ( kstp - 1 ), -1 ) |
---|
512 | |
---|
513 | zta_tmp (:,:,:) = ta (:,:,:) |
---|
514 | zsa_tmp (:,:,:) = sa (:,:,:) |
---|
515 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
516 | ! Computation of diagnostic variables |
---|
517 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
518 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
519 | !----------------------------------------------------------------------- |
---|
520 | |
---|
521 | ! ocean surface freezing temperature ... not needed in adjoint |
---|
522 | |
---|
523 | CALL wzv_adj( kstp ) ! Vertical velocity |
---|
524 | IF( n_cla == 1 ) CALL div_cla_adj( kstp ) ! Cross Land Advection (Update Hor. divergence) |
---|
525 | |
---|
526 | |
---|
527 | CALL div_cur_adj( kstp ) ! Horizontal divergence & Relative vorticity |
---|
528 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
529 | ! Dynamics |
---|
530 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
531 | ! N.B. ta, sa arrays are used as workspace in this section |
---|
532 | !----------------------------------------------------------------------- |
---|
533 | ! apply dynamics assimilation increment ... not needed in adjoint |
---|
534 | indic=0 |
---|
535 | !i bug lbc sur emp |
---|
536 | |
---|
537 | ! vertical mesh at next time step ... not available |
---|
538 | CALL dyn_nxt_adj( kstp ) ! lateral velocity at next time step |
---|
539 | CALL dyn_spg_adj( kstp, indic ) ! surface pressure gradient |
---|
540 | !i |
---|
541 | IF( lk_dynspg_rl ) THEN |
---|
542 | ! surface pressure gradient at open boundaries ... not available |
---|
543 | ENDIF |
---|
544 | CALL dyn_zdf_adj( kstp ) ! vertical diffusion |
---|
545 | CALL dyn_hpg_adj( kstp ) ! horizontal gradient of Hydrostatic pressure |
---|
546 | #if defined key_agrif |
---|
547 | ! momemtum sponge ... not available |
---|
548 | #endif |
---|
549 | CALL dyn_ldf_adj( kstp ) ! lateral mixing |
---|
550 | CALL dyn_vor_adj( kstp ) ! vorticity term including Coriolis |
---|
551 | CALL dyn_adv_adj( kstp ) ! advection (vector or flux form) |
---|
552 | |
---|
553 | ta (:,:,:) = zta_tmp (:,:,:) |
---|
554 | sa (:,:,:) = zsa_tmp (:,:,:) |
---|
555 | |
---|
556 | ua_ad(:,:,:) = 0.e0 ! set tracer trends to zero |
---|
557 | va_ad(:,:,:) = 0.e0 |
---|
558 | |
---|
559 | IF( ln_dynhpg_imp ) THEN ! semi-implicit hpg |
---|
560 | IF( ln_zps ) CALL zps_hde_adj( kstp, ta, sa, ta_ad, sa_ad, rhd_ad,& ! Partial steps: time filtered hor. gradient |
---|
561 | & gtu_ad, gsu_ad, gru_ad, & ! of t, s, rd at the bottom ocean level |
---|
562 | & gtv_ad, gsv_ad, grv_ad ) |
---|
563 | CALL eos_adj( ta, sa, ta_ad, sa_ad, rhd_ad, rhop_ad ) ! Time-filtered in situ density used in dynhpg module |
---|
564 | ELSE ! centered hpg (default case) |
---|
565 | IF( ln_zps ) CALL zps_hde_adj( kstp, tb, sb, & |
---|
566 | & tb_ad, sb_ad,rhd_ad,& ! Partial steps: now horizontal gradient |
---|
567 | & gtu_ad, gsu_ad, gru_ad, & ! of t, s, rd at the bottom ocean level |
---|
568 | & gtv_ad, gsv_ad, grv_ad ) |
---|
569 | CALL eos_adj( tb, sb, tb_ad, sb_ad, rhd_ad, rhop_ad ) ! now (swap=before) in situ density for dynhpg module |
---|
570 | ENDIF |
---|
571 | |
---|
572 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
573 | ! Active tracers |
---|
574 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
575 | ! N.B. ua, va arrays are used as workspace in this section |
---|
576 | !----------------------------------------------------------------------- |
---|
577 | |
---|
578 | CALL tra_nxt_adj( kstp ) ! tracer fields at next time step |
---|
579 | ! update the new (t,s) fields by non |
---|
580 | ! penetrative convective adjustment ... not available |
---|
581 | CALL tra_zdf_adj( kstp ) ! vertical mixing |
---|
582 | |
---|
583 | #if defined key_agrif |
---|
584 | ! tracers sponge ... not available |
---|
585 | #endif |
---|
586 | |
---|
587 | CALL tra_ldf_adj( kstp ) ! lateral mixing |
---|
588 | |
---|
589 | ! IF( lk_zdfkpp ) CALL tra_kpp_adj( kstp ) ! KPP non-local tracer fluxes |
---|
590 | |
---|
591 | IF( n_cla == 1 ) CALL tra_cla_adj( kstp ) ! Cross Land Advection (Update Hor. advection) |
---|
592 | |
---|
593 | CALL tra_adv_adj( kstp ) ! horizontal & vertical advection |
---|
594 | |
---|
595 | IF( lk_tradmp ) CALL tra_dmp_adj( kstp ) ! internal damping trends |
---|
596 | ! advective (and/or diffusive) bottom boundary layer scheme ... currently not available |
---|
597 | ! diffusive bottom boundary layer scheme ... currently not available |
---|
598 | |
---|
599 | IF( lk_trabbc ) CALL tra_bbc_adj( kstp ) ! bottom heat flux |
---|
600 | |
---|
601 | IF( ln_traqsr ) CALL tra_qsr_adj( kstp ) ! penetrative solar radiation qsr |
---|
602 | |
---|
603 | CALL tra_sbc_adj( kstp ) ! surface boundary condition |
---|
604 | |
---|
605 | ta_ad(:,:,:) = 0.e0 ! set tracer trends to zero |
---|
606 | sa_ad(:,:,:) = 0.e0 |
---|
607 | ! apply tracer assimilation increment ... not needed in adjoint |
---|
608 | #if defined key_top |
---|
609 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
610 | ! Passive Tracer Model |
---|
611 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
612 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
613 | !----------------------------------------------------------------------- |
---|
614 | |
---|
615 | ! time-stepping... not needed in adjoint for the time being |
---|
616 | |
---|
617 | #endif |
---|
618 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
619 | ! Ocean physics update |
---|
620 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
621 | !----------------------------------------------------------------------- |
---|
622 | ! LATERAL PHYSICS |
---|
623 | !----------------------------------------------------------------------- |
---|
624 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
625 | !----------------------------------------------------------------------- |
---|
626 | #if defined key_traldf_c2d |
---|
627 | ! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_adj |
---|
628 | #endif |
---|
629 | ! before slope of the lateral mixing... not needed in adjoint with lk_zdfcst_adj |
---|
630 | !----------------------------------------------------------------------- |
---|
631 | ! VERTICAL PHYSICS |
---|
632 | !----------------------------------------------------------------------- |
---|
633 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
634 | !----------------------------------------------------------------------- |
---|
635 | ! mixed layer depth... not needed in adjoint with lk_zdfcst_adj |
---|
636 | ! bottom friction... not needed in adjoint with lk_zdfcst_adj |
---|
637 | ! double diffusive mixing ... not needed in adjoint with lk_zdfcst_adj |
---|
638 | ! enhanced vertical eddy diffusivity ... not needed in adjoint with lk_zdfcst_adj |
---|
639 | ! ! increase diffusivity at rivers mouths... not needed in tangent |
---|
640 | ! lk_zdfcst_adj: Constant Kz read from the reference trajectory |
---|
641 | ! Constant Kz (reset avt, avm[uv] to the background value)... not available |
---|
642 | ! KPP closure scheme for Kz ... not available |
---|
643 | ! TKE closure scheme for Kz ... not available |
---|
644 | ! Richardson number dependent Kz ... not available |
---|
645 | ! ! Vertical eddy viscosity and diffusivity coefficients |
---|
646 | CALL bn2_adj( tb, sb, tb_ad, sb_ad, rn2_ad ) ! before Brunt-Vaisala frequency |
---|
647 | |
---|
648 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
649 | ! Update data, open boundaries and Forcings |
---|
650 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
651 | ! Output the initial state and forcings ... not needed in adjoint |
---|
652 | ! IF( lk_bdy ) CALL bdy_dta_adj( kstp ) ! update dynamic and tracer data at unstructured open boundary |
---|
653 | ! compute phase velocities at open boundaries ... not needed for the time being, to investigate whenever we do obc. |
---|
654 | ! update dynamic and tracer data at open boundaries ... not needed for the time being, to investigate whenever we do obc. |
---|
655 | |
---|
656 | CALL sbc_adj ( kstp ) ! Sea boundary condition (including sea-ice) |
---|
657 | ! update 3D salinity data ... not needed in tangent |
---|
658 | ! update 3D temperature data ... not needed in adjoint |
---|
659 | |
---|
660 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
661 | ! Control, and restarts |
---|
662 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
663 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
664 | !----------------------------------------------------------------------- |
---|
665 | ! ! Time loop: control and print |
---|
666 | CALL stp_ctl_adj( kstp, indic, 0 ) |
---|
667 | IF( indic < 0 ) CALL ctl_stop( 'step_adj: indic < 0' ) |
---|
668 | |
---|
669 | ! close input ocean restart file ... not needed in adjoint |
---|
670 | ! write output ocean restart file... not needed in adjoint |
---|
671 | ! write open boundary restart file... not needed in adjoint |
---|
672 | |
---|
673 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
674 | ! diagnostics and outputs |
---|
675 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
676 | ! N.B. ua, va, ta, sa arrays are used as workspace in this section |
---|
677 | !----------------------------------------------------------------------- |
---|
678 | |
---|
679 | IF ( nstop == 0 ) THEN ! Diagnostics |
---|
680 | ! drifting Floats... not needed in adjoint |
---|
681 | ! trends: dynamics ... not needed in adjoint |
---|
682 | ! trends: active tracers... not needed in adjoint |
---|
683 | ! trends: Mixed-layer ... not needed in adjoint |
---|
684 | ! trends: vorticity budget... not needed in adjoint |
---|
685 | ! Surface pressure diagnostics... not needed in adjoint |
---|
686 | ! Thermocline depth (20 degres isotherm depth)... not needed in adjoint |
---|
687 | ! basin averaged diagnostics... not needed in adjoint |
---|
688 | ! dynamical heigh diagnostics... not needed in adjoint |
---|
689 | ! Fresh water budget diagnostics... not needed in adjoint |
---|
690 | !!! IF( lk_diaobs ) CALL dia_obs_adj( kstp ) ! obs-minus-model (assimilation) diagnostics NOTE: may be better off outside this module |
---|
691 | ! Poleward TRansports diagnostics... not needed in adjoint |
---|
692 | |
---|
693 | |
---|
694 | ! ! Outputs |
---|
695 | ! ocean model: outputs... not needed in adjoint |
---|
696 | |
---|
697 | ENDIF |
---|
698 | |
---|
699 | !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> |
---|
700 | ! Coupled mode |
---|
701 | !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< |
---|
702 | |
---|
703 | ! coupled mode : field exchanges ... not available for the next 20 years |
---|
704 | ! |
---|
705 | END SUBROUTINE stp_adj |
---|
706 | |
---|
707 | SUBROUTINE stp_adj_tst( kumadt ) |
---|
708 | !!----------------------------------------------------------------------- |
---|
709 | !! |
---|
710 | !! *** ROUTINE example_adj_tst *** |
---|
711 | !! |
---|
712 | !! ** Purpose : Test the adjoint routine. |
---|
713 | !! |
---|
714 | !! ** Method : Verify the scalar product |
---|
715 | !! |
---|
716 | !! ( L dx )^T W dy = dx^T L^T W dy |
---|
717 | !! |
---|
718 | !! where L = tangent routine |
---|
719 | !! L^T = adjoint routine |
---|
720 | !! W = diagonal matrix of scale factors |
---|
721 | !! dx = input perturbation (random field) |
---|
722 | !! dy = L dx |
---|
723 | !! |
---|
724 | !! |
---|
725 | !! History : |
---|
726 | !! ! 09-03 (F. Vigilant) |
---|
727 | !!----------------------------------------------------------------------- |
---|
728 | !! * Modules used |
---|
729 | USE sbcssr_tam, ONLY: & |
---|
730 | & qrp_tl, qrp_ad, & |
---|
731 | & erp_tl, erp_ad |
---|
732 | USE sbcfwb_tam, ONLY: & |
---|
733 | & a_fwb_tl, & ! for before year. |
---|
734 | & a_fwb_ad ! for before year. |
---|
735 | |
---|
736 | !! * Arguments |
---|
737 | INTEGER, INTENT(IN) :: & |
---|
738 | & kumadt ! Output unit |
---|
739 | |
---|
740 | !! * Local declarations |
---|
741 | INTEGER :: & |
---|
742 | & ji, & ! dummy loop indices |
---|
743 | & jj, & |
---|
744 | & jk, & |
---|
745 | & istp |
---|
746 | INTEGER, DIMENSION(jpi,jpj) :: & |
---|
747 | & iseed_2d ! 2D seed for the random number generator |
---|
748 | REAL(KIND=wp) :: & |
---|
749 | & zsp1 , & ! scalar product involving the tangent routine |
---|
750 | & zsp1_U , & |
---|
751 | & zsp1_V , & |
---|
752 | & zsp1_T , & |
---|
753 | & zsp1_S , & |
---|
754 | & zsp1_SSH, & |
---|
755 | & zsp2 , & ! scalar product involving the adjoint routine |
---|
756 | & zsp2_U , & |
---|
757 | & zsp2_V , & |
---|
758 | & zsp2_T , & |
---|
759 | & zsp2_S , & |
---|
760 | & zsp2_SSH |
---|
761 | REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
762 | & zun_tlin , & ! Tangent input |
---|
763 | & zvn_tlin , & ! Tangent input |
---|
764 | & ztn_tlin , & ! Tangent input |
---|
765 | & zsn_tlin , & ! Tangent input |
---|
766 | & zun_tlout, & ! Tangent output |
---|
767 | & zvn_tlout, & ! Tangent output |
---|
768 | & ztn_tlout, & ! Tangent output |
---|
769 | & zsn_tlout, & ! Tangent outpu |
---|
770 | & zun_adin , & ! Adjoint input |
---|
771 | & zvn_adin , & ! Adjoint input |
---|
772 | & ztn_adin , & ! Adjoint input |
---|
773 | & zsn_adin , & ! Adjoint input |
---|
774 | #if defined key_obc |
---|
775 | & zub_adin , & ! Adjoint input |
---|
776 | & zvb_adin , & ! Adjoint input |
---|
777 | & ztb_adin , & ! Adjoint input |
---|
778 | & zsb_adin , & ! Adjoint input |
---|
779 | & zub_tlout , & ! Adjoint input |
---|
780 | & zvb_tlout , & ! Adjoint input |
---|
781 | & ztb_tlout , & ! Adjoint input |
---|
782 | & zsb_tlout , & ! Adjoint input |
---|
783 | #endif |
---|
784 | & zun_adout, & ! Adjoint output |
---|
785 | & zvn_adout, & ! Adjoint output |
---|
786 | & ztn_adout, & ! Adjoint output |
---|
787 | & zsn_adout, & ! Adjoint output |
---|
788 | & z3r ! 3D random field |
---|
789 | REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
790 | & zsshn_tlin , & ! Tangent input |
---|
791 | & zsshn_tlout , & ! Tangent input |
---|
792 | & zsshn_adin , & ! Adjoint output |
---|
793 | & zsshn_adout , & ! Adjoint output |
---|
794 | #if defined key_obc |
---|
795 | & zsshb_tlout , & ! Tangent input |
---|
796 | & zsshb_adin , & ! Tangent input |
---|
797 | #endif |
---|
798 | & z2r ! 2D random field |
---|
799 | |
---|
800 | CHARACTER(LEN=14) :: & |
---|
801 | & cl_name |
---|
802 | |
---|
803 | INTEGER :: & |
---|
804 | & jpert |
---|
805 | INTEGER, PARAMETER :: jpertmax = 6 |
---|
806 | |
---|
807 | ! Allocate memory |
---|
808 | |
---|
809 | ALLOCATE( & |
---|
810 | & zun_tlin( jpi,jpj,jpk), zvn_tlin( jpi,jpj,jpk), ztn_tlin( jpi,jpj,jpk), & |
---|
811 | & zsn_tlin( jpi,jpj,jpk), zsshn_tlin( jpi,jpj ), zun_tlout( jpi,jpj,jpk), & |
---|
812 | & zvn_tlout( jpi,jpj,jpk), ztn_tlout( jpi,jpj,jpk), zsn_tlout( jpi,jpj,jpk), & |
---|
813 | & zsshn_tlout(jpi,jpj ), zun_adin( jpi,jpj,jpk), zvn_adin( jpi,jpj,jpk), & |
---|
814 | & ztn_adin( jpi,jpj,jpk), zsn_adin( jpi,jpj,jpk), zsshn_adin(jpi,jpj ), & |
---|
815 | & zun_adout( jpi,jpj,jpk), zvn_adout( jpi,jpj,jpk), ztn_adout( jpi,jpj,jpk), & |
---|
816 | & zsn_adout( jpi,jpj,jpk), zsshn_adout(jpi,jpj ), z2r( jpi,jpj ), & |
---|
817 | & z3r( jpi,jpj,jpk) & |
---|
818 | & ) |
---|
819 | |
---|
820 | #if defined key_obc |
---|
821 | ALLOCATE( zub_adin (jpi,jpj,jpk), zvb_adin (jpi,jpj,jpk) , & |
---|
822 | & ztb_adin (jpi,jpj,jpk), zsb_adin (jpi,jpj,jpk) , & |
---|
823 | & zub_tlout (jpi,jpj,jpk), zvb_tlout (jpi,jpj,jpk) , & |
---|
824 | & ztb_tlout (jpi,jpj,jpk), zsb_tlout (jpi,jpj,jpk) , & |
---|
825 | & zsshb_tlout(jpi,jpj) , zsshb_adin(jpi,jpj) ) |
---|
826 | #endif |
---|
827 | !================================================================== |
---|
828 | ! 1) dx = ( un_tl, vn_tl, tn_tl, sn_tl, sshn_tl ) and |
---|
829 | ! dy = ( ..... ) |
---|
830 | !================================================================== |
---|
831 | |
---|
832 | |
---|
833 | DO jpert = 1, jpertmax |
---|
834 | |
---|
835 | !-------------------------------------------------------------------- |
---|
836 | ! Reset the tangent and adjoint variables |
---|
837 | !-------------------------------------------------------------------- |
---|
838 | zun_tlin (:,:,:) = 0.0_wp |
---|
839 | zvn_tlin (:,:,:) = 0.0_wp |
---|
840 | ztn_tlin (:,:,:) = 0.0_wp |
---|
841 | zsn_tlin (:,:,:) = 0.0_wp |
---|
842 | zsshn_tlin ( :,:) = 0.0_wp |
---|
843 | zun_tlout (:,:,:) = 0.0_wp |
---|
844 | zvn_tlout (:,:,:) = 0.0_wp |
---|
845 | ztn_tlout (:,:,:) = 0.0_wp |
---|
846 | zsn_tlout (:,:,:) = 0.0_wp |
---|
847 | zsshn_tlout( :,:) = 0.0_wp |
---|
848 | zun_adin (:,:,:) = 0.0_wp |
---|
849 | zvn_adin (:,:,:) = 0.0_wp |
---|
850 | ztn_adin (:,:,:) = 0.0_wp |
---|
851 | zsn_adin (:,:,:) = 0.0_wp |
---|
852 | zsshn_adin ( :,:) = 0.0_wp |
---|
853 | zun_adout (:,:,:) = 0.0_wp |
---|
854 | zvn_adout (:,:,:) = 0.0_wp |
---|
855 | ztn_adout (:,:,:) = 0.0_wp |
---|
856 | zsn_adout (:,:,:) = 0.0_wp |
---|
857 | zsshn_adout( :,:) = 0.0_wp |
---|
858 | z3r (:,:,:) = 0.0_wp |
---|
859 | z2r ( :,:) = 0.0_wp |
---|
860 | |
---|
861 | !-------------------------------------------------------------------- |
---|
862 | ! Initialize the tangent input with random noise: dx |
---|
863 | !-------------------------------------------------------------------- |
---|
864 | |
---|
865 | IF ( (jpert == 1) .OR. (jpert == jpertmax) ) THEN |
---|
866 | DO jj = 1, jpj |
---|
867 | DO ji = 1, jpi |
---|
868 | iseed_2d(ji,jj) = - ( 596035 + & |
---|
869 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
870 | END DO |
---|
871 | END DO |
---|
872 | CALL grid_random( iseed_2d, z3r, 'U', 0.0_wp, stdu ) |
---|
873 | DO jk = 1, jpk |
---|
874 | DO jj = nldj, nlej |
---|
875 | DO ji = nldi, nlei |
---|
876 | zun_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
877 | END DO |
---|
878 | END DO |
---|
879 | END DO |
---|
880 | ENDIF |
---|
881 | IF ( (jpert == 2) .OR. (jpert == jpertmax) ) THEN |
---|
882 | DO jj = 1, jpj |
---|
883 | DO ji = 1, jpi |
---|
884 | iseed_2d(ji,jj) = - ( 371836 + & |
---|
885 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
886 | END DO |
---|
887 | END DO |
---|
888 | CALL grid_random( iseed_2d, z3r, 'V', 0.0_wp, stdv ) |
---|
889 | DO jk = 1, jpk |
---|
890 | DO jj = nldj, nlej |
---|
891 | DO ji = nldi, nlei |
---|
892 | zvn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
893 | END DO |
---|
894 | END DO |
---|
895 | END DO |
---|
896 | ENDIF |
---|
897 | IF ( (jpert == 3) .OR. (jpert == jpertmax) ) THEN |
---|
898 | DO jj = 1, jpj |
---|
899 | DO ji = 1, jpi |
---|
900 | iseed_2d(ji,jj) = - ( 284035 + & |
---|
901 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
902 | END DO |
---|
903 | END DO |
---|
904 | CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stdt ) |
---|
905 | DO jk = 1, jpk |
---|
906 | DO jj = nldj, nlej |
---|
907 | DO ji = nldi, nlei |
---|
908 | ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
909 | END DO |
---|
910 | END DO |
---|
911 | END DO |
---|
912 | ENDIF |
---|
913 | IF ( (jpert == 4) .OR. (jpert == jpertmax) ) THEN |
---|
914 | DO jj = 1, jpj |
---|
915 | DO ji = 1, jpi |
---|
916 | iseed_2d(ji,jj) = - ( 471426 + & |
---|
917 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
918 | END DO |
---|
919 | END DO |
---|
920 | CALL grid_random( iseed_2d, z3r, 'T', 0.0_wp, stds ) |
---|
921 | DO jk = 1, jpk |
---|
922 | DO jj = nldj, nlej |
---|
923 | DO ji = nldi, nlei |
---|
924 | zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
925 | END DO |
---|
926 | END DO |
---|
927 | END DO |
---|
928 | ENDIF |
---|
929 | IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN |
---|
930 | DO jj = 1, jpj |
---|
931 | DO ji = 1, jpi |
---|
932 | iseed_2d(ji,jj) = - ( 712035 + & |
---|
933 | & mig(ji) + ( mjg(jj) - 1 ) * jpiglo ) |
---|
934 | END DO |
---|
935 | END DO |
---|
936 | CALL grid_random( iseed_2d, z2r, 'T', 0.0_wp, stdssh ) |
---|
937 | DO jj = nldj, nlej |
---|
938 | DO ji = nldi, nlei |
---|
939 | zsshn_tlin(ji,jj) = z2r(ji,jj) |
---|
940 | END DO |
---|
941 | END DO |
---|
942 | ENDIF |
---|
943 | IF (jpert>0) THEN |
---|
944 | CALL oce_tam_init( 0 ) ! kindic = 0 allocate/initialize tl and adj variables |
---|
945 | CALL sbc_oce_tam_init( 0 ) ! kindic = 0 allocate/initialize tl and adj variables |
---|
946 | CALL sol_oce_tam_init( 1 ) |
---|
947 | CALL sol_oce_tam_init( 2 ) |
---|
948 | CALL trc_oce_tam_init( 0 ) |
---|
949 | |
---|
950 | qrp_tl = 0.0_wp |
---|
951 | qrp_ad = 0.0_wp |
---|
952 | erp_tl = 0.0_wp |
---|
953 | erp_ad = 0.0_wp |
---|
954 | |
---|
955 | emp_tl = 0.e0 |
---|
956 | emp_ad = 0.e0 |
---|
957 | |
---|
958 | a_fwb_tl = 0.0_wp |
---|
959 | a_fwb_ad = 0.0_wp |
---|
960 | |
---|
961 | istp = nit000 - 1 |
---|
962 | CALL trj_rea( istp, 1 ) |
---|
963 | !-------------------------------------------------------------------- |
---|
964 | ! Initialize the tangent variables: dy^* = W dy |
---|
965 | !-------------------------------------------------------------------- |
---|
966 | |
---|
967 | un_tl (:,:,:) = zun_tlin (:,:,:) |
---|
968 | vn_tl (:,:,:) = zvn_tlin (:,:,:) |
---|
969 | tn_tl (:,:,:) = ztn_tlin (:,:,:) |
---|
970 | sn_tl (:,:,:) = zsn_tlin (:,:,:) |
---|
971 | sshn_tl( :,:) = zsshn_tlin ( :,:) |
---|
972 | |
---|
973 | #if defined key_pomme_r025 |
---|
974 | IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN |
---|
975 | !DO ji = 1, jpi |
---|
976 | ! DO jj = 1, jpj |
---|
977 | ! sshn_tl(ji,jj) = cos( (2.*rpi)*(FLOAT(ji)-0.5)/FLOAT(jpi) + rpi/2. ) & |
---|
978 | ! * sin( (2.*rpi)*(FLOAT(jj)-0.5)/FLOAT(jpj) ) / 100. |
---|
979 | ! END DO |
---|
980 | !END DO |
---|
981 | |
---|
982 | DO ji = 1, jpi |
---|
983 | DO jj = 1, jpj |
---|
984 | zsshn_tlin(ji,jj) = exp( -((float(ji)-float(jpi)/2.)/(float(jpi)/5.))**2 & |
---|
985 | -((float(jj)-float(jpj)/2.)/(float(jpj)/5.))**2 ) / 100. & |
---|
986 | * tmask(ji,jj,1) |
---|
987 | END DO |
---|
988 | END DO |
---|
989 | sshn_tl(:,:) = zsshn_tlin(:,:) |
---|
990 | ENDIF |
---|
991 | #endif |
---|
992 | |
---|
993 | !----------------------------------------------------------------------- |
---|
994 | ! Initialization of the dynamics and tracer fields for the tangent |
---|
995 | !----------------------------------------------------------------------- |
---|
996 | |
---|
997 | CALL istate_init_tan |
---|
998 | |
---|
999 | DO istp = nit000, nitend, 1 |
---|
1000 | CALL stp_tan( istp ) |
---|
1001 | END DO |
---|
1002 | |
---|
1003 | |
---|
1004 | zun_tlout ( :,:,:) = un_tl (:,:,:) |
---|
1005 | zvn_tlout ( :,:,:) = vn_tl (:,:,:) |
---|
1006 | ztn_tlout ( :,:,:) = tn_tl (:,:,:) |
---|
1007 | zsn_tlout ( :,:,:) = sn_tl (:,:,:) |
---|
1008 | zsshn_tlout( :,:) = sshn_tl ( :,:) |
---|
1009 | |
---|
1010 | #if defined key_obc |
---|
1011 | zub_tlout (:,:,:) = ub_tl (:,:,:) |
---|
1012 | zvb_tlout (:,:,:) = vb_tl (:,:,:) |
---|
1013 | ztb_tlout (:,:,:) = tb_tl (:,:,:) |
---|
1014 | zsb_tlout (:,:,:) = sb_tl (:,:,:) |
---|
1015 | zsshb_tlout(:,:) = sshb_tl(:,:) |
---|
1016 | #endif |
---|
1017 | |
---|
1018 | !-------------------------------------------------------------------- |
---|
1019 | ! Initialize the adjoint variables: dy^* = W dy |
---|
1020 | !-------------------------------------------------------------------- |
---|
1021 | |
---|
1022 | DO jk = 1, jpk |
---|
1023 | DO jj = nldj, nlej |
---|
1024 | DO ji = nldi, nlei |
---|
1025 | zun_adin(ji,jj,jk) = zun_tlout(ji,jj,jk) & |
---|
1026 | & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & |
---|
1027 | & * umask(ji,jj,jk) * wesp_u |
---|
1028 | #if defined key_obc |
---|
1029 | zub_adin(ji,jj,jk) = zub_tlout(ji,jj,jk) & |
---|
1030 | & * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) & |
---|
1031 | & * umask(ji,jj,jk) * wesp_u |
---|
1032 | #endif |
---|
1033 | END DO |
---|
1034 | END DO |
---|
1035 | END DO |
---|
1036 | |
---|
1037 | DO jk = 1, jpk |
---|
1038 | DO jj = nldj, nlej |
---|
1039 | DO ji = nldi, nlei |
---|
1040 | zvn_adin(ji,jj,jk) = zvn_tlout(ji,jj,jk) & |
---|
1041 | & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) & |
---|
1042 | & * vmask(ji,jj,jk) * wesp_u |
---|
1043 | #if defined key_obc |
---|
1044 | zvb_adin(ji,jj,jk) = zvb_tlout(ji,jj,jk) & |
---|
1045 | & * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) & |
---|
1046 | & * vmask(ji,jj,jk) * wesp_u |
---|
1047 | #endif |
---|
1048 | END DO |
---|
1049 | END DO |
---|
1050 | END DO |
---|
1051 | |
---|
1052 | DO jk = 1, jpk |
---|
1053 | DO jj = nldj, nlej |
---|
1054 | DO ji = nldi, nlei |
---|
1055 | ztn_adin(ji,jj,jk) = ztn_tlout(ji,jj,jk) & |
---|
1056 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) & |
---|
1057 | & * tmask(ji,jj,jk) * wesp_t(jk) |
---|
1058 | #if defined key_obc |
---|
1059 | ztb_adin(ji,jj,jk) = ztb_tlout(ji,jj,jk) & |
---|
1060 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) & |
---|
1061 | & * tmask(ji,jj,jk) * wesp_t(jk) |
---|
1062 | #endif |
---|
1063 | END DO |
---|
1064 | END DO |
---|
1065 | END DO |
---|
1066 | |
---|
1067 | DO jk = 1, jpk |
---|
1068 | DO jj = nldj, nlej |
---|
1069 | DO ji = nldi, nlei |
---|
1070 | zsn_adin(ji,jj,jk) = zsn_tlout(ji,jj,jk) & |
---|
1071 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) & |
---|
1072 | & * tmask(ji,jj,jk) * wesp_s(jk) |
---|
1073 | #if defined key_obc |
---|
1074 | zsb_adin(ji,jj,jk) = zsb_tlout(ji,jj,jk) & |
---|
1075 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) & |
---|
1076 | & * tmask(ji,jj,jk) * wesp_s(jk) |
---|
1077 | #endif |
---|
1078 | END DO |
---|
1079 | END DO |
---|
1080 | END DO |
---|
1081 | |
---|
1082 | DO jj = nldj, nlej |
---|
1083 | DO ji = nldi, nlei |
---|
1084 | zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) & |
---|
1085 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) & |
---|
1086 | & * tmask(ji,jj,1) * wesp_ssh |
---|
1087 | #if defined key_obc |
---|
1088 | zsshb_adin(ji,jj) = zsshb_tlout(ji,jj) & |
---|
1089 | & * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) & |
---|
1090 | & * tmask(ji,jj,1) * wesp_ssh |
---|
1091 | #endif |
---|
1092 | END DO |
---|
1093 | END DO |
---|
1094 | |
---|
1095 | !-------------------------------------------------------------------- |
---|
1096 | ! Compute the scalar product: ( L dx )^T W dy |
---|
1097 | !-------------------------------------------------------------------- |
---|
1098 | |
---|
1099 | |
---|
1100 | zsp1_U = DOT_PRODUCT( zun_tlout , zun_adin ) |
---|
1101 | zsp1_V = DOT_PRODUCT( zvn_tlout , zvn_adin ) |
---|
1102 | zsp1_T = DOT_PRODUCT( ztn_tlout , ztn_adin ) |
---|
1103 | zsp1_S = DOT_PRODUCT( zsn_tlout , zsn_adin ) |
---|
1104 | zsp1_SSH = DOT_PRODUCT( zsshn_tlout, zsshn_adin ) |
---|
1105 | |
---|
1106 | zsp1 = zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH |
---|
1107 | |
---|
1108 | #if defined key_obc |
---|
1109 | zsp1_U = DOT_PRODUCT( zub_tlout , zub_adin ) |
---|
1110 | zsp1_V = DOT_PRODUCT( zvb_tlout , zvb_adin ) |
---|
1111 | zsp1_T = DOT_PRODUCT( ztb_tlout , ztb_adin ) |
---|
1112 | zsp1_S = DOT_PRODUCT( zsb_tlout , zsb_adin ) |
---|
1113 | zsp1_SSH = DOT_PRODUCT( zsshb_tlout, zsshb_adin ) |
---|
1114 | |
---|
1115 | zsp1 = zsp1 + ( zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH ) |
---|
1116 | #endif |
---|
1117 | !-------------------------------------------------------------------- |
---|
1118 | ! Call the adjoint routine: dx^* = L^T dy^* |
---|
1119 | !-------------------------------------------------------------------- |
---|
1120 | |
---|
1121 | un_ad (:,:,:) = zun_adin (:,:,:) |
---|
1122 | vn_ad (:,:,:) = zvn_adin (:,:,:) |
---|
1123 | tn_ad (:,:,:) = ztn_adin (:,:,:) |
---|
1124 | sn_ad (:,:,:) = zsn_adin (:,:,:) |
---|
1125 | sshn_ad( :,:) = zsshn_adin ( :,:) |
---|
1126 | |
---|
1127 | #if defined key_obc |
---|
1128 | ub_ad (:,:,:) = zub_adin (:,:,:) |
---|
1129 | vb_ad (:,:,:) = zvb_adin (:,:,:) |
---|
1130 | tb_ad (:,:,:) = ztb_adin (:,:,:) |
---|
1131 | sb_ad (:,:,:) = zsb_adin (:,:,:) |
---|
1132 | sshb_ad(:,:) = zsshb_adin (:,:) |
---|
1133 | #endif |
---|
1134 | |
---|
1135 | DO istp = nitend, nit000, -1 |
---|
1136 | CALL stp_adj(istp) |
---|
1137 | END DO |
---|
1138 | |
---|
1139 | CALL istate_init_adj |
---|
1140 | |
---|
1141 | zun_adout ( :,:,:) = un_ad (:,:,:) |
---|
1142 | zvn_adout ( :,:,:) = vn_ad (:,:,:) |
---|
1143 | ztn_adout ( :,:,:) = tn_ad (:,:,:) |
---|
1144 | zsn_adout ( :,:,:) = sn_ad (:,:,:) |
---|
1145 | zsshn_adout( :,:) = sshn_ad ( :,:) |
---|
1146 | |
---|
1147 | zsp2_U = DOT_PRODUCT( zun_tlin , zun_adout ) |
---|
1148 | zsp2_V = DOT_PRODUCT( zvn_tlin , zvn_adout ) |
---|
1149 | zsp2_T = DOT_PRODUCT( ztn_tlin , ztn_adout ) |
---|
1150 | zsp2_S = DOT_PRODUCT( zsn_tlin , zsn_adout ) |
---|
1151 | zsp2_SSH = DOT_PRODUCT( zsshn_tlin, zsshn_adout ) |
---|
1152 | |
---|
1153 | zsp2 = zsp2_U + zsp2_V + zsp2_T + zsp2_S + zsp2_SSH |
---|
1154 | |
---|
1155 | ! 14 char:'12345678901234' |
---|
1156 | SELECT CASE (jpert) |
---|
1157 | CASE(1) |
---|
1158 | cl_name = 'step_adj U' |
---|
1159 | CASE(2) |
---|
1160 | cl_name = 'step_adj V' |
---|
1161 | CASE(3) |
---|
1162 | cl_name = 'step_adj T' |
---|
1163 | CASE(4) |
---|
1164 | cl_name = 'step_adj S' |
---|
1165 | CASE(5) |
---|
1166 | cl_name = 'step_adj SSH' |
---|
1167 | CASE(jpertmax) |
---|
1168 | cl_name = 'step_adj ' |
---|
1169 | END SELECT |
---|
1170 | CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 ) |
---|
1171 | endif |
---|
1172 | END DO |
---|
1173 | |
---|
1174 | DEALLOCATE( & |
---|
1175 | & zun_tlin , zvn_tlin , ztn_tlin , & |
---|
1176 | & zsn_tlin , zsshn_tlin , zun_tlout , & |
---|
1177 | & zvn_tlout , ztn_tlout , zsn_tlout , & |
---|
1178 | & zsshn_tlout, zun_adin , zvn_adin , & |
---|
1179 | & ztn_adin , zsn_adin , zsshn_adin , & |
---|
1180 | & zun_adout , zvn_adout , ztn_adout , & |
---|
1181 | & zsn_adout , zsshn_adout, z2r , & |
---|
1182 | & z3r & |
---|
1183 | & ) |
---|
1184 | |
---|
1185 | END SUBROUTINE stp_adj_tst |
---|
1186 | |
---|
1187 | SUBROUTINE stp_tlm_tst( kumadt ) |
---|
1188 | !!----------------------------------------------------------------------- |
---|
1189 | !! |
---|
1190 | !! *** ROUTINE example_tl_tst *** |
---|
1191 | !! |
---|
1192 | !! ** Purpose : Test the tangent linear routine. |
---|
1193 | !! |
---|
1194 | !! ** Method : Verify the tangent with Taylor expansion |
---|
1195 | !! |
---|
1196 | !! M(x+hdx) = M(x) + L(hdx) + O(h^2) |
---|
1197 | !! |
---|
1198 | !! where L = tangent routine |
---|
1199 | !! M = direct routine |
---|
1200 | !! dx = input perturbation (random field) |
---|
1201 | !! h = ration on perturbation |
---|
1202 | !! |
---|
1203 | !! In the tangent test we verify that: |
---|
1204 | !! M(x+h*dx) - M(x) |
---|
1205 | !! g(h) = ------------------ ---> 1 as h ---> 0 |
---|
1206 | !! L(h*dx) |
---|
1207 | !! and |
---|
1208 | !! g(h) - 1 |
---|
1209 | !! f(h) = ---------- ---> k (costant) as h ---> 0 |
---|
1210 | !! p |
---|
1211 | !! |
---|
1212 | !! History : |
---|
1213 | !! ! 09-06 (F. Vigilant) |
---|
1214 | !!----------------------------------------------------------------------- |
---|
1215 | !! * Modules used |
---|
1216 | USE opatam_tst_ini, ONLY: & |
---|
1217 | & tlm_namrd |
---|
1218 | USE sbcssr_tam, ONLY: & |
---|
1219 | & qrp_tl,erp_tl |
---|
1220 | USE sbcfwb_tam, ONLY: & |
---|
1221 | & a_fwb_tl ! for before year. |
---|
1222 | USE step_c1d ! Time stepping loop for the 1D configuration |
---|
1223 | USE tamtrj ! writing out state trajectory |
---|
1224 | USE lib_mpp, ONLY: & ! distributed memory computing |
---|
1225 | & lk_mpp, & |
---|
1226 | & mpp_max |
---|
1227 | USE c1d, ONLY: & ! 1D configuration |
---|
1228 | & lk_c1d |
---|
1229 | USE par_tlm, ONLY: & |
---|
1230 | & cur_loop, & |
---|
1231 | & h_ratio |
---|
1232 | USE istate_mod |
---|
1233 | USE wzvmod ! vertical velocity |
---|
1234 | USE sbc_oce , ONLY: & ! ocean dynamics and tracers variables |
---|
1235 | & emp, & |
---|
1236 | & emps |
---|
1237 | USE tamctl, ONLY: & ! Control parameters |
---|
1238 | & numtan, numtan_sc |
---|
1239 | USE step |
---|
1240 | !! * Arguments |
---|
1241 | INTEGER, INTENT(IN) :: & |
---|
1242 | & kumadt ! Output unit |
---|
1243 | |
---|
1244 | !! * Local declarations |
---|
1245 | INTEGER :: & |
---|
1246 | & ji, & ! dummy loop indices |
---|
1247 | & jj, & |
---|
1248 | & jk, & |
---|
1249 | & istp, & |
---|
1250 | & pgamma, & ! control the number of pertubation loops |
---|
1251 | & igamma, & |
---|
1252 | & nit_loop |
---|
1253 | INTEGER, DIMENSION(jpi,jpj) :: & |
---|
1254 | & iseed_2d ! 2D seed for the random number generator |
---|
1255 | REAL(KIND=wp) :: & |
---|
1256 | & zsp1, zsp2, zsp3, zzsp, & ! scalar product |
---|
1257 | & zsp1_U, zsp1_V, zsp1_T, zsp1_S, zsp1_SSH, & |
---|
1258 | & zsp2_U, zsp2_V, zsp2_T, zsp2_S, zsp2_SSH, & |
---|
1259 | & zsp3_U, zsp3_V, zsp3_T, zsp3_S, zsp3_SSH, & |
---|
1260 | & zzsp_U, zzsp_V, zzsp_T, zzsp_S, zzsp_SSH, & |
---|
1261 | & gamma, & |
---|
1262 | & zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, & |
---|
1263 | & zgsp6, zgsp7 |
---|
1264 | REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: & |
---|
1265 | & ta_tmp , & |
---|
1266 | & sa_tmp , & |
---|
1267 | & zun_tlin , & ! Pertubation input |
---|
1268 | & zvn_tlin , & ! Pertubation input |
---|
1269 | & ztn_tlin , & ! Pertubation input |
---|
1270 | & zsn_tlin , & ! Pertubation input |
---|
1271 | & zun_out , & ! Direct output |
---|
1272 | & zvn_out , & ! Direct output |
---|
1273 | & ztn_out , & ! Direct output |
---|
1274 | & zsn_out , & ! Direct output |
---|
1275 | & zun_wop , & ! Direct without pertubation output |
---|
1276 | & zvn_wop , & ! Direct without pertubation output |
---|
1277 | & ztn_wop , & ! Direct without pertubation output |
---|
1278 | & zsn_wop , & ! Direct without pertubation output |
---|
1279 | & z3r ! 3D random field |
---|
1280 | REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: & |
---|
1281 | & emp_tmp, & |
---|
1282 | & emps_tmp, & |
---|
1283 | & sshn_tmp , & |
---|
1284 | & zsshn_tlin , & ! Pertubation input |
---|
1285 | & zsshn_out , & ! Direct output |
---|
1286 | & zsshn_wop , & ! Direct without pertubation output |
---|
1287 | & z2r ! 2D random field |
---|
1288 | |
---|
1289 | CHARACTER(LEN=14) :: cl_name |
---|
1290 | CHARACTER (LEN=128) :: file_out, file_wop, file_wop2 |
---|
1291 | CHARACTER (LEN=90) :: FMT |
---|
1292 | |
---|
1293 | REAL(KIND=wp), DIMENSION(100):: & |
---|
1294 | & zscu , zscv , zsct , zscs , zscssh, & |
---|
1295 | & zscerru, zscerrv, zscerrt, zscerrs, zscerrssh |
---|
1296 | INTEGER, DIMENSION(100):: & |
---|
1297 | & iiposu, iiposv, iipost, iiposs, iiposssh, & |
---|
1298 | & ijposu, ijposv, ijpost, ijposs, ijposssh, & |
---|
1299 | & ikposu, ikposv, ikpost, ikposs |
---|
1300 | INTEGER:: & |
---|
1301 | & ii, & |
---|
1302 | & isamp=40, & |
---|
1303 | & jsamp=40, & |
---|
1304 | & ksamp=10, & |
---|
1305 | & numsctlm |
---|
1306 | REAL(KIND=wp), DIMENSION(jpi,jpj,jpk) :: & |
---|
1307 | & zerru, zerrv, zerrt, zerrs ! Pertubation input |
---|
1308 | REAL(KIND=wp), DIMENSION(jpi,jpj) :: & |
---|
1309 | & zerrssh ! Pertubation input |
---|
1310 | INTEGER :: & |
---|
1311 | & jdir, jdirmax=6 ! Perturbation direction (U, V, T, S, SSH, all) |
---|
1312 | |
---|
1313 | ! Allocate memory |
---|
1314 | |
---|
1315 | ALLOCATE( & |
---|
1316 | & zun_tlin( jpi,jpj,jpk), & |
---|
1317 | & zvn_tlin( jpi,jpj,jpk), & |
---|
1318 | & ztn_tlin( jpi,jpj,jpk), & |
---|
1319 | & zsn_tlin( jpi,jpj,jpk), & |
---|
1320 | & zsshn_tlin( jpi,jpj ), & |
---|
1321 | & zun_out( jpi,jpj,jpk), & |
---|
1322 | & zvn_out( jpi,jpj,jpk), & |
---|
1323 | & ztn_out( jpi,jpj,jpk), & |
---|
1324 | & zsn_out( jpi,jpj,jpk), & |
---|
1325 | & zsshn_out( jpi,jpj ), & |
---|
1326 | & zun_wop( jpi,jpj,jpk), & |
---|
1327 | & zvn_wop( jpi,jpj,jpk), & |
---|
1328 | & ztn_wop( jpi,jpj,jpk), & |
---|
1329 | & zsn_wop( jpi,jpj,jpk), & |
---|
1330 | & zsshn_wop( jpi,jpj ), & |
---|
1331 | & z2r( jpi,jpj ), & |
---|
1332 | & z3r( jpi,jpj,jpk), & |
---|
1333 | & emp_tmp( jpi,jpj ), & |
---|
1334 | & emps_tmp( jpi,jpj ), & |
---|
1335 | & ta_tmp( jpi,jpj,jpk), & |
---|
1336 | & sa_tmp( jpi,jpj,jpk), & |
---|
1337 | & sshn_tmp( jpi,jpj ) & |
---|
1338 | & ) |
---|
1339 | |
---|
1340 | !-------------------------------------------------------------------- |
---|
1341 | ! Reset the scalar |
---|
1342 | !-------------------------------------------------------------------- |
---|
1343 | zscu (:) = 0.0_wp |
---|
1344 | zscv (:) = 0.0_wp |
---|
1345 | zsct (:) = 0.0_wp |
---|
1346 | zscs (:) = 0.0_wp |
---|
1347 | zscssh (:) = 0.0_wp |
---|
1348 | |
---|
1349 | iiposu (:) = 0 |
---|
1350 | iiposv (:) = 0 |
---|
1351 | iipost (:) = 0 |
---|
1352 | iiposs (:) = 0 |
---|
1353 | iiposssh (:) = 0 |
---|
1354 | ijposu (:) = 0 |
---|
1355 | ijposv (:) = 0 |
---|
1356 | ijpost (:) = 0 |
---|
1357 | ijposs (:) = 0 |
---|
1358 | ijposssh (:) = 0 |
---|
1359 | ikposu (:) = 0 |
---|
1360 | ikposv (:) = 0 |
---|
1361 | ikpost (:) = 0 |
---|
1362 | ikposs (:) = 0 |
---|
1363 | |
---|
1364 | zscerru (:) = 0.0_wp |
---|
1365 | zscerrv (:) = 0.0_wp |
---|
1366 | zscerrt (:) = 0.0_wp |
---|
1367 | zscerrs (:) = 0.0_wp |
---|
1368 | zscerrssh (:) = 0.0_wp |
---|
1369 | |
---|
1370 | zerru (:,:,:) = 0.0_wp |
---|
1371 | zerrv (:,:,:) = 0.0_wp |
---|
1372 | zerrt (:,:,:) = 0.0_wp |
---|
1373 | zerrs (:,:,:) = 0.0_wp |
---|
1374 | zerrssh (:,:) = 0.0_wp |
---|
1375 | |
---|
1376 | !-------------------------------------------------------------------- |
---|
1377 | ! Reset the tangent variables |
---|
1378 | !-------------------------------------------------------------------- |
---|
1379 | zun_tlin ( :,:,:) = 0.0_wp |
---|
1380 | zvn_tlin ( :,:,:) = 0.0_wp |
---|
1381 | ztn_tlin ( :,:,:) = 0.0_wp |
---|
1382 | zsn_tlin ( :,:,:) = 0.0_wp |
---|
1383 | zsshn_tlin ( :,:) = 0.0_wp |
---|
1384 | z3r ( :,:,:) = 0.0_wp |
---|
1385 | z2r ( :,:) = 0.0_wp |
---|
1386 | |
---|
1387 | !-------------------------------------------------------------------- |
---|
1388 | ! Output filename Xn=F(X0) |
---|
1389 | !-------------------------------------------------------------------- |
---|
1390 | file_wop='trj_wop_step' |
---|
1391 | CALL tlm_namrd |
---|
1392 | gamma = h_ratio |
---|
1393 | !-------------------------------------------------------------------- |
---|
1394 | ! Initialize the tangent input with random noise: dx |
---|
1395 | !-------------------------------------------------------------------- |
---|
1396 | jdir = jdirmax |
---|
1397 | IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN |
---|
1398 | IF (( jdir .EQ. 1 ) .OR. ( jdir .EQ. jdirmax ) ) THEN |
---|
1399 | CALL grid_rd_sd( 596035, z3r, 'U', 0.0_wp, stdu) |
---|
1400 | DO jk = 1, jpk |
---|
1401 | DO jj = nldj, nlej |
---|
1402 | DO ji = nldi, nlei |
---|
1403 | zun_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
1404 | END DO |
---|
1405 | END DO |
---|
1406 | END DO |
---|
1407 | |
---|
1408 | ENDIF |
---|
1409 | IF (( jdir .EQ. 2 ) .OR. ( jdir .EQ. jdirmax ) ) THEN |
---|
1410 | CALL grid_rd_sd( 371836, z3r, 'V', 0.0_wp, stdv) |
---|
1411 | DO jk = 1, jpk |
---|
1412 | DO jj = nldj, nlej |
---|
1413 | DO ji = nldi, nlei |
---|
1414 | zvn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
1415 | END DO |
---|
1416 | END DO |
---|
1417 | END DO |
---|
1418 | |
---|
1419 | ENDIF |
---|
1420 | IF (( jdir .EQ. 3 ) .OR. ( jdir .EQ. jdirmax ) ) THEN |
---|
1421 | CALL grid_rd_sd( 284035, z3r, 'T', 0.0_wp, stdt) |
---|
1422 | DO jk = 1, jpk |
---|
1423 | DO jj = nldj, nlej |
---|
1424 | DO ji = nldi, nlei |
---|
1425 | ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
1426 | END DO |
---|
1427 | END DO |
---|
1428 | END DO |
---|
1429 | |
---|
1430 | ENDIF |
---|
1431 | IF (( jdir .EQ. 4 ) .OR. ( jdir .EQ. jdirmax ) ) THEN |
---|
1432 | CALL grid_rd_sd( 471426, z3r, 'T', 0.0_wp, stds) |
---|
1433 | DO jk = 1, jpk |
---|
1434 | DO jj = nldj, nlej |
---|
1435 | DO ji = nldi, nlei |
---|
1436 | zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk) |
---|
1437 | END DO |
---|
1438 | END DO |
---|
1439 | END DO |
---|
1440 | |
---|
1441 | ENDIF |
---|
1442 | IF (( jdir .EQ. 5 ) .OR. ( jdir .EQ. jdirmax ) ) THEN |
---|
1443 | CALL grid_rd_sd( 712035, z2r,'T', 0.0_wp, stdssh) |
---|
1444 | DO jj = nldj, nlej |
---|
1445 | DO ji = nldi, nlei |
---|
1446 | zsshn_tlin(ji,jj) = z2r(ji,jj) |
---|
1447 | END DO |
---|
1448 | END DO |
---|
1449 | |
---|
1450 | ENDIF |
---|
1451 | ENDIF |
---|
1452 | CALL istate_p |
---|
1453 | |
---|
1454 | !-------------------------------------------------------------------- |
---|
1455 | ! Open File to save scalar result |
---|
1456 | !-------------------------------------------------------------------- |
---|
1457 | |
---|
1458 | PRINT*,'IN TST_STP h_ratio, cur_loop, gamma', h_ratio, ' ',cur_loop,' ', gamma |
---|
1459 | Call flush(numout) |
---|
1460 | |
---|
1461 | ! check that all process are still there... If some process have an error, |
---|
1462 | ! they will never enter in step and other processes will wait until the end of the cpu time! |
---|
1463 | IF( lk_mpp ) CALL mpp_max( nstop ) |
---|
1464 | |
---|
1465 | istp = nit000 |
---|
1466 | IF( lk_c1d ) THEN ! 1D configuration (no AGRIF zoom) |
---|
1467 | ! |
---|
1468 | DO WHILE ( istp <= nitend .AND. nstop == 0 ) |
---|
1469 | CALL stp_c1d( istp ) |
---|
1470 | istp = istp + 1 |
---|
1471 | END DO |
---|
1472 | ELSE |
---|
1473 | istp = nit000 - 1 |
---|
1474 | IF( ln_trjwri ) CALL tam_trj_wri( istp ) ! Output trajectory fields |
---|
1475 | ENDIF |
---|
1476 | |
---|
1477 | IF (( cur_loop .NE. 0) .OR. ( gamma .NE. 0.0_wp) )THEN |
---|
1478 | zun_tlin(:,:,:) = gamma * zun_tlin(:,:,:) |
---|
1479 | un(:,:,:) = un(:,:,:) + zun_tlin(:,:,:) |
---|
1480 | |
---|
1481 | zvn_tlin(:,:,:) = gamma * zvn_tlin(:,:,:) |
---|
1482 | vn(:,:,:) = vn(:,:,:) + zvn_tlin(:,:,:) |
---|
1483 | |
---|
1484 | ztn_tlin(:,:,:) = gamma * ztn_tlin(:,:,:) |
---|
1485 | tn(:,:,:) = tn(:,:,:) + ztn_tlin(:,:,:) |
---|
1486 | |
---|
1487 | zsn_tlin(:,:,:) = gamma * zsn_tlin(:,:,:) |
---|
1488 | sn(:,:,:) = sn(:,:,:) + zsn_tlin(:,:,:) |
---|
1489 | |
---|
1490 | zsshn_tlin(:,:) = gamma * zsshn_tlin(:,:) |
---|
1491 | sshn (:,:) = sshn (:,:) + zsshn_tlin(:,:) |
---|
1492 | ENDIF |
---|
1493 | |
---|
1494 | IF( .NOT. lk_vvl ) CALL wzv(nit000) |
---|
1495 | |
---|
1496 | !-------------------------------------------------------------------- |
---|
1497 | ! Compute the direct model F(X0,t=n) = Xn |
---|
1498 | !-------------------------------------------------------------------- |
---|
1499 | DO istp = nit000, nitend, 1 |
---|
1500 | CALL stp( istp ) |
---|
1501 | END DO |
---|
1502 | IF ( cur_loop .EQ. 0) CALL trj_wri_spl(file_wop) |
---|
1503 | |
---|
1504 | !-------------------------------------------------------------------- |
---|
1505 | ! Compute the Tangent |
---|
1506 | !-------------------------------------------------------------------- |
---|
1507 | IF ( cur_loop .NE. 0) THEN |
---|
1508 | |
---|
1509 | !-------------------------------------------------------------------- |
---|
1510 | ! Storing data |
---|
1511 | !-------------------------------------------------------------------- |
---|
1512 | zun_out (:,:,:) = un (:,:,:) |
---|
1513 | zvn_out (:,:,:) = vn (:,:,:) |
---|
1514 | ztn_out (:,:,:) = tn (:,:,:) |
---|
1515 | zsn_out (:,:,:) = sn (:,:,:) |
---|
1516 | zsshn_out (:,: ) = sshn (:,: ) |
---|
1517 | |
---|
1518 | !-------------------------------------------------------------------- |
---|
1519 | ! Initialize the tangent variables: dy^* = W dy |
---|
1520 | !-------------------------------------------------------------------- |
---|
1521 | qrp_tl = 0.0_wp |
---|
1522 | |
---|
1523 | a_fwb_tl = 0.0_wp |
---|
1524 | |
---|
1525 | CALL trj_rea( nit000-1, 1 ) |
---|
1526 | |
---|
1527 | un_tl (:,:,:) = zun_tlin (:,:,:) |
---|
1528 | vn_tl (:,:,:) = zvn_tlin (:,:,:) |
---|
1529 | tn_tl (:,:,:) = ztn_tlin (:,:,:) |
---|
1530 | sn_tl (:,:,:) = zsn_tlin (:,:,:) |
---|
1531 | sshn_tl( :,:) = zsshn_tlin ( :,:) |
---|
1532 | |
---|
1533 | |
---|
1534 | !----------------------------------------------------------------------- |
---|
1535 | ! Initialization of the dynamics and tracer fields for the tangent |
---|
1536 | !----------------------------------------------------------------------- |
---|
1537 | CALL istate_init_tan |
---|
1538 | DO istp = nit000, nitend, 1 |
---|
1539 | CALL stp_tan( istp ) |
---|
1540 | END DO |
---|
1541 | |
---|
1542 | !-------------------------------------------------------------------- |
---|
1543 | ! Compute the scalar product: ( L(t0,tn) gamma dx0 ) ) |
---|
1544 | !-------------------------------------------------------------------- |
---|
1545 | |
---|
1546 | zsp2_U = DOT_PRODUCT( un_tl , un_tl ) |
---|
1547 | zsp2_V = DOT_PRODUCT( vn_tl , vn_tl ) |
---|
1548 | zsp2_T = DOT_PRODUCT( tn_tl , tn_tl ) |
---|
1549 | zsp2_S = DOT_PRODUCT( sn_tl , sn_tl ) |
---|
1550 | zsp2_SSH = DOT_PRODUCT( sshn_tl, sshn_tl ) |
---|
1551 | |
---|
1552 | zsp2 = zsp2_U + zsp2_V + zsp2_T + zsp2_S + zsp2_SSH |
---|
1553 | |
---|
1554 | !-------------------------------------------------------------------- |
---|
1555 | ! Storing data |
---|
1556 | !-------------------------------------------------------------------- |
---|
1557 | CALL trj_rd_spl(file_wop) |
---|
1558 | |
---|
1559 | zun_wop (:,:,:) = un (:,:,:) |
---|
1560 | zvn_wop (:,:,:) = vn (:,:,:) |
---|
1561 | ztn_wop (:,:,:) = tn (:,:,:) |
---|
1562 | zsn_wop (:,:,:) = sn (:,:,:) |
---|
1563 | zsshn_wop (:,: ) = sshn (:,: ) |
---|
1564 | !-------------------------------------------------------------------- |
---|
1565 | ! Compute the Linearization Error |
---|
1566 | ! Nn = M( X0+gamma.dX0, t0,tn) - M(X0, t0,tn) |
---|
1567 | ! and |
---|
1568 | ! Compute the Linearization Error |
---|
1569 | ! En = Nn -TL(gamma.dX0, t0,tn) |
---|
1570 | !-------------------------------------------------------------------- |
---|
1571 | ! Warning: Here we re-use local variables z()_out and z()_wop |
---|
1572 | ii=0 |
---|
1573 | DO jk = 1, jpk |
---|
1574 | DO jj = 1, jpj |
---|
1575 | DO ji = 1, jpi |
---|
1576 | zun_out (ji,jj,jk) = zun_out (ji,jj,jk) - zun_wop (ji,jj,jk) |
---|
1577 | zun_wop (ji,jj,jk) = zun_out (ji,jj,jk) - un_tl (ji,jj,jk) |
---|
1578 | IF ( un_tl(ji,jj,jk) .NE. 0.0_wp ) & |
---|
1579 | & zerru(ji,jj,jk) = zun_out(ji,jj,jk)/un_tl(ji,jj,jk) |
---|
1580 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1581 | & (MOD(jj, jsamp) .EQ. 0) .AND. & |
---|
1582 | & (MOD(jk, ksamp) .EQ. 0) ) THEN |
---|
1583 | ii = ii+1 |
---|
1584 | iiposu(ii) = ji |
---|
1585 | ijposu(ii) = jj |
---|
1586 | ikposu(ii) = jk |
---|
1587 | IF ( INT(umask(ji,jj,jk)) .NE. 0) THEN |
---|
1588 | zscu (ii) = zun_out(ji,jj,jk) |
---|
1589 | zscerru (ii) = ( zerru(ji,jj,jk) - 1.0_wp ) / gamma |
---|
1590 | ENDIF |
---|
1591 | ENDIF |
---|
1592 | END DO |
---|
1593 | END DO |
---|
1594 | END DO |
---|
1595 | ii=0 |
---|
1596 | DO jk = 1, jpk |
---|
1597 | DO jj = 1, jpj |
---|
1598 | DO ji = 1, jpi |
---|
1599 | zvn_out (ji,jj,jk) = zvn_out (ji,jj,jk) - zvn_wop (ji,jj,jk) |
---|
1600 | zvn_wop (ji,jj,jk) = zvn_out (ji,jj,jk) - vn_tl (ji,jj,jk) |
---|
1601 | IF ( vn_tl(ji,jj,jk) .NE. 0.0_wp ) & |
---|
1602 | & zerrv(ji,jj,jk) = zvn_out(ji,jj,jk)/vn_tl(ji,jj,jk) |
---|
1603 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1604 | & (MOD(jj, jsamp) .EQ. 0) .AND. & |
---|
1605 | & (MOD(jk, ksamp) .EQ. 0) ) THEN |
---|
1606 | ii = ii+1 |
---|
1607 | iiposv(ii) = ji |
---|
1608 | ijposv(ii) = jj |
---|
1609 | ikposv(ii) = jk |
---|
1610 | IF ( INT(vmask(ji,jj,jk)) .NE. 0) THEN |
---|
1611 | zscv (ii) = zvn_out (ji,jj,jk) |
---|
1612 | zscerrv (ii) = ( zerrv(ji,jj,jk) - 1.0_wp ) / gamma |
---|
1613 | ENDIF |
---|
1614 | ENDIF |
---|
1615 | END DO |
---|
1616 | END DO |
---|
1617 | END DO |
---|
1618 | ii=0 |
---|
1619 | DO jk = 1, jpk |
---|
1620 | DO jj = 1, jpj |
---|
1621 | DO ji = 1, jpi |
---|
1622 | ztn_out (ji,jj,jk) = ztn_out (ji,jj,jk) - ztn_wop (ji,jj,jk) |
---|
1623 | ztn_wop (ji,jj,jk) = ztn_out (ji,jj,jk) - tn_tl (ji,jj,jk) |
---|
1624 | IF ( tn_tl(ji,jj,jk) .NE. 0.0_wp ) & |
---|
1625 | & zerrt(ji,jj,jk) = ztn_out(ji,jj,jk)/tn_tl(ji,jj,jk) |
---|
1626 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1627 | & (MOD(jj, jsamp) .EQ. 0) .AND. & |
---|
1628 | & (MOD(jk, ksamp) .EQ. 0) ) THEN |
---|
1629 | ii = ii+1 |
---|
1630 | iipost(ii) = ji |
---|
1631 | ijpost(ii) = jj |
---|
1632 | ikpost(ii) = jk |
---|
1633 | IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN |
---|
1634 | zsct (ii) = ztn_out (ji,jj,jk) |
---|
1635 | zscerrt (ii) = ( zerrt(ji,jj,jk) - 1.0_wp ) / gamma |
---|
1636 | ENDIF |
---|
1637 | ENDIF |
---|
1638 | END DO |
---|
1639 | END DO |
---|
1640 | END DO |
---|
1641 | ii=0 |
---|
1642 | DO jk = 1, jpk |
---|
1643 | DO jj = 1, jpj |
---|
1644 | DO ji = 1, jpi |
---|
1645 | zsn_out (ji,jj,jk) = zsn_out (ji,jj,jk) - zsn_wop (ji,jj,jk) |
---|
1646 | zsn_wop (ji,jj,jk) = zsn_out (ji,jj,jk) - sn_tl (ji,jj,jk) |
---|
1647 | IF ( sn_tl(ji,jj,jk) .NE. 0.0_wp ) & |
---|
1648 | & zerrs(ji,jj,jk) = zsn_out(ji,jj,jk)/sn_tl(ji,jj,jk) |
---|
1649 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1650 | & (MOD(jj, jsamp) .EQ. 0) .AND. & |
---|
1651 | & (MOD(jk, ksamp) .EQ. 0) ) THEN |
---|
1652 | ii = ii+1 |
---|
1653 | iiposs(ii) = ji |
---|
1654 | ijposs(ii) = jj |
---|
1655 | ikposs(ii) = jk |
---|
1656 | IF ( INT(tmask(ji,jj,jk)) .NE. 0) THEN |
---|
1657 | zscs (ii) = zsn_out (ji,jj,jk) |
---|
1658 | zscerrs (ii) = ( zerrs(ji,jj,jk) - 1.0_wp ) / gamma |
---|
1659 | ENDIF |
---|
1660 | ENDIF |
---|
1661 | END DO |
---|
1662 | END DO |
---|
1663 | END DO |
---|
1664 | ii=0 |
---|
1665 | jk=1 |
---|
1666 | DO jj = 1, jpj |
---|
1667 | DO ji = 1, jpi |
---|
1668 | zsshn_out (ji,jj ) = zsshn_out (ji,jj ) - zsshn_wop (ji,jj ) |
---|
1669 | zsshn_wop (ji,jj ) = zsshn_out (ji,jj ) - sshn_tl (ji,jj ) |
---|
1670 | IF ( sshn_tl(ji,jj) .NE. 0.0_wp ) & |
---|
1671 | & zerrssh(ji,jj) = zsshn_out(ji,jj)/sshn_tl(ji,jj) |
---|
1672 | IF( (MOD(ji, isamp) .EQ. 0) .AND. & |
---|
1673 | & (MOD(jj, jsamp) .EQ. 0) ) THEN |
---|
1674 | ii = ii+1 |
---|
1675 | iiposssh(ii) = ji |
---|
1676 | ijposssh(ii) = jj |
---|
1677 | IF ( INT(tmask(ji,jj,jk)) .NE. 0 ) THEN |
---|
1678 | zscssh (ii) = zsshn_out (ji,jj) |
---|
1679 | zscerrssh (ii) = ( zerrssh(ji,jj) - 1.0_wp ) / gamma |
---|
1680 | ENDIF |
---|
1681 | ENDIF |
---|
1682 | END DO |
---|
1683 | END DO |
---|
1684 | |
---|
1685 | zsp1_U = DOT_PRODUCT( zun_out , zun_out ) |
---|
1686 | zsp1_V = DOT_PRODUCT( zvn_out , zvn_out ) |
---|
1687 | zsp1_T = DOT_PRODUCT( ztn_out , ztn_out ) |
---|
1688 | zsp1_S = DOT_PRODUCT( zsn_out , zsn_out ) |
---|
1689 | zsp1_SSH = DOT_PRODUCT( zsshn_out, zsshn_out ) |
---|
1690 | |
---|
1691 | zsp1 = zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH |
---|
1692 | |
---|
1693 | zsp3_U = DOT_PRODUCT( zun_wop , zun_wop ) |
---|
1694 | zsp3_V = DOT_PRODUCT( zvn_wop , zvn_wop ) |
---|
1695 | zsp3_T = DOT_PRODUCT( ztn_wop , ztn_wop ) |
---|
1696 | zsp3_S = DOT_PRODUCT( zsn_wop , zsn_wop ) |
---|
1697 | zsp3_SSH = DOT_PRODUCT( zsshn_wop, zsshn_wop ) |
---|
1698 | |
---|
1699 | zsp3 = zsp3_U + zsp3_V + zsp3_T + zsp3_S + zsp3_SSH |
---|
1700 | |
---|
1701 | !-------------------------------------------------------------------- |
---|
1702 | ! Print the linearization error En - norme 2 |
---|
1703 | !-------------------------------------------------------------------- |
---|
1704 | ! 14 char:'12345678901234' |
---|
1705 | cl_name = 'step_tam:En ' |
---|
1706 | zzsp = SQRT(zsp3) |
---|
1707 | zzsp_U = SQRT(zsp3_U) |
---|
1708 | zzsp_V = SQRT(zsp3_V) |
---|
1709 | zzsp_T = SQRT(zsp3_T) |
---|
1710 | zzsp_S = SQRT(zsp3_S) |
---|
1711 | zzsp_SSH = SQRT(zsp3_SSH) |
---|
1712 | zgsp5 = zzsp |
---|
1713 | CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) |
---|
1714 | |
---|
1715 | CALL prntst_tlm( 'U ', kumadt, zzsp_U , h_ratio ) |
---|
1716 | CALL prntst_tlm( 'V ', kumadt, zzsp_V , h_ratio ) |
---|
1717 | CALL prntst_tlm( 'T ', kumadt, zzsp_T , h_ratio ) |
---|
1718 | CALL prntst_tlm( 'S ', kumadt, zzsp_S , h_ratio ) |
---|
1719 | CALL prntst_tlm( 'SSH ', kumadt, zzsp_SSH, h_ratio ) |
---|
1720 | |
---|
1721 | !-------------------------------------------------------------------- |
---|
1722 | ! Compute TLM norm2 |
---|
1723 | !-------------------------------------------------------------------- |
---|
1724 | zzsp = SQRT(zsp2) |
---|
1725 | zzsp_U = SQRT(zsp2_U) |
---|
1726 | zzsp_V = SQRT(zsp2_V) |
---|
1727 | zzsp_T = SQRT(zsp2_T) |
---|
1728 | zzsp_S = SQRT(zsp2_S) |
---|
1729 | zzsp_SSH = SQRT(zsp2_SSH) |
---|
1730 | zgsp4 = zzsp |
---|
1731 | cl_name = 'step_tam:Ln2 ' |
---|
1732 | CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) |
---|
1733 | |
---|
1734 | CALL prntst_tlm( 'U ', kumadt, zzsp_U , h_ratio ) |
---|
1735 | CALL prntst_tlm( 'V ', kumadt, zzsp_V , h_ratio ) |
---|
1736 | CALL prntst_tlm( 'T ', kumadt, zzsp_T , h_ratio ) |
---|
1737 | CALL prntst_tlm( 'S ', kumadt, zzsp_S , h_ratio ) |
---|
1738 | CALL prntst_tlm( 'SSH ', kumadt, zzsp_SSH, h_ratio ) |
---|
1739 | |
---|
1740 | !-------------------------------------------------------------------- |
---|
1741 | ! Print the linearization error Nn - norme 2 |
---|
1742 | !-------------------------------------------------------------------- |
---|
1743 | ! 14 char:'12345678901234' |
---|
1744 | cl_name = 'step_tam:Nn ' |
---|
1745 | zzsp = SQRT(zsp1) |
---|
1746 | zzsp_U = SQRT(zsp1_U) |
---|
1747 | zzsp_V = SQRT(zsp1_V) |
---|
1748 | zzsp_T = SQRT(zsp1_T) |
---|
1749 | zzsp_S = SQRT(zsp1_S) |
---|
1750 | zzsp_SSH = SQRT(zsp1_SSH) |
---|
1751 | ! zzsp = zzsp_U + zzsp_V + zzsp_T + zzsp_S + zzsp_SSH |
---|
1752 | CALL prntst_tlm( cl_name, kumadt, zzsp, h_ratio ) |
---|
1753 | |
---|
1754 | CALL prntst_tlm( 'U ', kumadt, zzsp_U , h_ratio ) |
---|
1755 | CALL prntst_tlm( 'V ', kumadt, zzsp_V , h_ratio ) |
---|
1756 | CALL prntst_tlm( 'T ', kumadt, zzsp_T , h_ratio ) |
---|
1757 | CALL prntst_tlm( 'S ', kumadt, zzsp_S , h_ratio ) |
---|
1758 | CALL prntst_tlm( 'SSH ', kumadt, zzsp_SSH, h_ratio ) |
---|
1759 | |
---|
1760 | zgsp3 = SQRT( zsp3/zsp2 ) |
---|
1761 | zgsp7 = zgsp3/gamma |
---|
1762 | zgsp1 = zzsp |
---|
1763 | zgsp2 = zgsp1 / zgsp4 |
---|
1764 | zgsp6 = (zgsp2 - 1.0_wp)/gamma |
---|
1765 | |
---|
1766 | FMT = "(A8,2X,I4.4,2X,E6.1,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13,2X,E20.13)" |
---|
1767 | WRITE(numtan,FMT) 'step ', cur_loop, h_ratio, zgsp1, zgsp2, zgsp3, zgsp4, zgsp5, zgsp6, zgsp7 |
---|
1768 | |
---|
1769 | !-------------------------------------------------------------------- |
---|
1770 | ! Unitary calculus |
---|
1771 | !-------------------------------------------------------------------- |
---|
1772 | FMT = "(A8,2X,A8,2X,I4.4,2X,E6.1,2X,I4.4,2X,I4.4,2X,I4.4,2X,E20.13,1X)" |
---|
1773 | cl_name = 'step ' |
---|
1774 | IF (lwp) THEN |
---|
1775 | DO ii=1, 100, 1 |
---|
1776 | IF ( zscu(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscu ', & |
---|
1777 | & cur_loop, h_ratio, iiposu(ii), ijposu(ii), ikposu(ii), zscu(ii) |
---|
1778 | ENDDO |
---|
1779 | DO ii=1, 100, 1 |
---|
1780 | IF ( zscv(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscv ', & |
---|
1781 | & cur_loop, h_ratio, iiposv(ii), ijposv(ii), ikposv(ii), zscv(ii) |
---|
1782 | ENDDO |
---|
1783 | DO ii=1, 100, 1 |
---|
1784 | IF ( zsct(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zsct ', & |
---|
1785 | & cur_loop, h_ratio, iipost(ii), ijpost(ii), ikpost(ii), zsct(ii) |
---|
1786 | ENDDO |
---|
1787 | DO ii=1, 100, 1 |
---|
1788 | IF ( zscs(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscs ', & |
---|
1789 | & cur_loop, h_ratio, iiposs(ii), ijposs(ii), ikposs(ii), zscs(ii) |
---|
1790 | ENDDO |
---|
1791 | DO ii=1, 100, 1 |
---|
1792 | IF ( zscssh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscssh ', & |
---|
1793 | & cur_loop, h_ratio, iiposssh(ii), ijposssh(ii), ijposssh(ii), zscssh(ii) |
---|
1794 | ENDDO |
---|
1795 | |
---|
1796 | DO ii=1, 100, 1 |
---|
1797 | IF ( zscerru(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerru ', & |
---|
1798 | & cur_loop, h_ratio, iiposu(ii), ijposu(ii), ikposu(ii), zscerru(ii) |
---|
1799 | ENDDO |
---|
1800 | DO ii=1, 100, 1 |
---|
1801 | IF ( zscerrv(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrv ', & |
---|
1802 | & cur_loop, h_ratio, iiposv(ii), ijposv(ii), ikposv(ii), zscerrv(ii) |
---|
1803 | ENDDO |
---|
1804 | DO ii=1, 100, 1 |
---|
1805 | IF ( zscerrt(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrt ', & |
---|
1806 | & cur_loop, h_ratio, iipost(ii), ijpost(ii), ikpost(ii), zscerrt(ii) |
---|
1807 | ENDDO |
---|
1808 | DO ii=1, 100, 1 |
---|
1809 | IF ( zscerrs(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerrs ', & |
---|
1810 | & cur_loop, h_ratio, iiposs(ii), ijposs(ii), ikposs(ii), zscerrs(ii) |
---|
1811 | ENDDO |
---|
1812 | DO ii=1, 100, 1 |
---|
1813 | IF ( zscerrssh(ii) .NE. 0.0_wp ) WRITE(numtan_sc,FMT) cl_name, 'zscerssh', & |
---|
1814 | & cur_loop, h_ratio, iiposssh(ii), ijposssh(ii), ijposssh(ii), zscerrssh(ii) |
---|
1815 | ENDDO |
---|
1816 | ! write separator |
---|
1817 | WRITE(numtan_sc,"(A4)") '====' |
---|
1818 | ENDIF |
---|
1819 | |
---|
1820 | ENDIF |
---|
1821 | |
---|
1822 | DEALLOCATE( & |
---|
1823 | & zun_tlin , & |
---|
1824 | & zvn_tlin , & |
---|
1825 | & ztn_tlin , & |
---|
1826 | & zsn_tlin , & |
---|
1827 | & zsshn_tlin , & |
---|
1828 | & zun_out , & |
---|
1829 | & zvn_out , & |
---|
1830 | & ztn_out , & |
---|
1831 | & zsn_out , & |
---|
1832 | & zsshn_out , & |
---|
1833 | & zun_wop , & |
---|
1834 | & zvn_wop , & |
---|
1835 | & ztn_wop , & |
---|
1836 | & zsn_wop , & |
---|
1837 | & zsshn_wop , & |
---|
1838 | & z2r , & |
---|
1839 | & z3r & |
---|
1840 | & ) |
---|
1841 | |
---|
1842 | END SUBROUTINE stp_tlm_tst |
---|
1843 | |
---|
1844 | !!---------------------------------------------------------------------- |
---|
1845 | |
---|
1846 | SUBROUTINE asm_trj_wop_tasa_wri |
---|
1847 | !!----------------------------------------------------------------------- |
---|
1848 | !! |
---|
1849 | !! *** ROUTINE asm_trj_wop_tasa_wri *** |
---|
1850 | !! |
---|
1851 | !! ** Purpose : Write to file the model state trajectory for use with |
---|
1852 | !! 4D-Var. |
---|
1853 | !! |
---|
1854 | !! ** Method : |
---|
1855 | !! |
---|
1856 | !! ** Action : |
---|
1857 | !! |
---|
1858 | !! History : |
---|
1859 | !! ! 09-07 (F. Vigilant) |
---|
1860 | !!----------------------------------------------------------------------- |
---|
1861 | !! *Module udes |
---|
1862 | USE iom |
---|
1863 | USE oce , ONLY: & ! ocean dynamics and tracers variables |
---|
1864 | & tn, sn |
---|
1865 | !! * Arguments |
---|
1866 | !! * Local declarations |
---|
1867 | INTEGER :: & |
---|
1868 | & inum ! File unit number |
---|
1869 | CHARACTER (LEN=50) :: & |
---|
1870 | & filename |
---|
1871 | ! Define the output file |
---|
1872 | filename='trj_wop_tasa' |
---|
1873 | WRITE(filename, FMT='(A,A)' ) TRIM( filename ), '.nc' |
---|
1874 | filename = TRIM( filename ) |
---|
1875 | CALL iom_open( filename, inum, ldwrt = .TRUE., kiolib = jprstlib) |
---|
1876 | |
---|
1877 | ! Output trajectory fields |
---|
1878 | CALL iom_rstput( 1, 1, inum, 'tn' , tn ) |
---|
1879 | CALL iom_rstput( 1, 1, inum, 'sn' , sn ) |
---|
1880 | CALL iom_close( inum ) |
---|
1881 | END SUBROUTINE asm_trj_wop_tasa_wri |
---|
1882 | |
---|
1883 | SUBROUTINE asm_trj_wop_tasa_rd |
---|
1884 | !!----------------------------------------------------------------------- |
---|
1885 | !! |
---|
1886 | !! *** ROUTINE asm_trj_wop_tasa_rd *** |
---|
1887 | !! |
---|
1888 | !! ** Purpose : Read file the model state trajectory for use with |
---|
1889 | !! 4D-Var. |
---|
1890 | !! |
---|
1891 | !! ** Method : |
---|
1892 | !! |
---|
1893 | !! ** Action : |
---|
1894 | !! |
---|
1895 | !! History : |
---|
1896 | !! ! 09-07 (F. Vigilant) |
---|
1897 | !!----------------------------------------------------------------------- |
---|
1898 | !! *Module udes |
---|
1899 | USE iom ! I/O module |
---|
1900 | USE oce , ONLY: & ! ocean dynamics and tracers variables |
---|
1901 | & tn, sn |
---|
1902 | !! * Arguments |
---|
1903 | !! * Local declarations |
---|
1904 | INTEGER :: & |
---|
1905 | & inum ! File unit number |
---|
1906 | CHARACTER (LEN=50) :: & |
---|
1907 | & filename |
---|
1908 | ! Define the output file |
---|
1909 | filename='trj_wop_tasa' |
---|
1910 | WRITE(filename, FMT='(A,A)' ) TRIM( filename ), '.nc' |
---|
1911 | filename = TRIM( filename ) |
---|
1912 | CALL iom_open( filename, inum) |
---|
1913 | |
---|
1914 | ! Output trajectory fields |
---|
1915 | CALL iom_get( inum, jpdom_data, 'tn' , tn, 1 ) |
---|
1916 | CALL iom_get( inum, jpdom_data, 'sn' , sn, 1 ) |
---|
1917 | CALL iom_close( inum ) |
---|
1918 | END SUBROUTINE asm_trj_wop_tasa_rd |
---|
1919 | |
---|
1920 | |
---|
1921 | !!====================================================================== |
---|
1922 | #endif |
---|
1923 | END MODULE step_tam |
---|