New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
step_tam.F90 in branches/2010_and_older/TAM_V3_2_2/NEMOTAM/OPATAM_SRC – NEMO

source: branches/2010_and_older/TAM_V3_2_2/NEMOTAM/OPATAM_SRC/step_tam.F90 @ 5226

Last change on this file since 5226 was 2578, checked in by rblod, 13 years ago

first import of NEMOTAM 3.2.2

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