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