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_cpd.F90 in branches/TAM_V3_0/NEMOTAM/OPATAM_SRC – NEMO

source: branches/TAM_V3_0/NEMOTAM/OPATAM_SRC/step_tam_cpd.F90 @ 2587

Last change on this file since 2587 was 2587, checked in by vidard, 13 years ago

refer to ticket #798

File size: 19.8 KB
Line 
1MODULE step_tam_cpd
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE step_tam  ***
5   !! Time-stepping    : manager of the adjoint ocean time stepping
6   !!                    Tangent and Adjoint module
7   !!======================================================================
8   !! History of
9   !! the direct:      !  91-03  ()  Original code
10   !!                  !  91-11  (G. Madec)
11   !!                  !  92-06  (M. Imbard)  add a first output record
12   !!                  !  96-04  (G. Madec)  introduction of dynspg
13   !!                  !  96-04  (M.A. Foujols)  introduction of passive tracer
14   !!             8.0  !  97-06  (G. Madec)  new architecture of call
15   !!             8.2  !  97-06  (G. Madec, M. Imbard, G. Roullet)  free surface
16   !!             8.2  !  99-02  (G. Madec, N. Grima)  hpg implicit
17   !!             8.2  !  00-07  (J-M Molines, M. Imbard)  Open Bondary Conditions
18   !!             9.0  !  02-06  (G. Madec)  free form, suppress macro-tasking
19   !!             " "  !  04-08  (C. Talandier) New trends organization
20   !!             " "  !  05-01  (C. Ethe) Add the KPP closure scheme
21   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization
22   !!             " "  !  05-11  (G. Madec)  Reorganisation of tra and dyn calls
23   !!             " "  !  06-01  (L. Debreu, C. Mazauric)  Agrif implementation
24   !!             " "  !  06-07  (S. Masson)  restart using iom
25   !!    "        " "  !  07-04  (K. Mogensen, A. Weaver, M. Martin) Assimilation interface
26   !! History of the TAM
27   !!                  !  08-06  (A. Vidard and A. Weaver) Tangent and Adjoint version of 9.0
28   !!                  !  08-11  (A. Vidard) Nemo v3 update
29   !!                  !  06-09  (F. Vigilant)  Modified to split NEMOVAR / NEMOTAM
30   !!----------------------------------------------------------------------
31   !!   stp_tam        : OPA system time-stepping (tangent linear)
32   !!   stp_adj        : OPA system time-stepping (adjoint)
33   !!----------------------------------------------------------------------
34
35   USE par_kind      , ONLY: & ! Precision variables
36      & wp
37   USE par_oce       , ONLY: & ! Ocean space and time domain variables
38      & jpi,                 &
39      & jpj,                 & 
40      & jpk,                 &
41      & jpiglo
42   USE oce           , ONLY: & ! ocean dynamics and tracers variables
43      & tb, sb, tn, sn, ta,  &
44      & un, vn, sshn, sshb,  &
45      & sa, ub, vb,          &
46      & ln_dynhpg_imp
47   USE zdfkpp        , ONLY: &
48      & lk_zdfkpp      ! KPP vertical mixing flag
49   USE dynspg_oce    , ONLY: &
50      &  lk_dynspg_rl  ! Rigid-lid flag
51   USE dom_oce       , ONLY: & ! ocean space and time domain variables
52      & n_cla, e1u, e2u,     &
53      & e1v, e2v, e1t, e2t,  &
54# if defined key_vvl
55      & e3t_1,               &
56# else
57#  if defined key_zco
58      & e3t_0, e3w_0,        &
59#  else
60      & e3t, e3u, e3v, e3w,  &
61#  endif
62# endif
63      & nldi, nldj, nlei,    &
64      & nlej, lk_vvl, ln_zps,&
65      & narea, mig, mjg,     &
66      & umask, vmask,  tmask
67   USE in_out_manager, ONLY: & ! I/O manager
68      & numout, lwp,         &
69      & nit000, nitend,      &
70      & nstop, ctl_stop
71   USE daymod        , ONLY: &
72      & adatrj
73
74   USE trabbc        , ONLY: &! bottom boundary condition
75      & lk_trabbc
76   USE traqsr        , ONLY: &
77      &  ln_traqsr
78          ! solar radiation penetration flag
79   USE oce_tam       , ONLY: & ! Tangent linear and adjoint variables
80      & oce_tam_init,        &
81      & un_tl, vn_tl, tn_tl, &
82      & sn_tl, ua_tl, va_tl, &
83      & ta_tl, sa_tl, ub_tl, &
84      & vb_tl, tb_tl, sb_tl, & 
85      & un_ad, vn_ad, tn_ad, &
86      & sn_ad, ua_ad, va_ad, &
87      & ta_ad, sa_ad, ub_ad, &
88      & vb_ad, tb_ad, sb_ad, & 
89      & sshn_tl, sshn_ad,    &
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   USE lbclnk_tam
101   USE daymod_tam      ! calendar                         (adjoint of day     routine)
102   USE sbc_oce_tam
103   USE sbcmod_tam
104   USE traqsr_tam      ! solar radiation penetration      (adjoint of tra_qsr routine)
105   USE trasbc_tam      ! surface boundary condition       (adjoint of tra_sbc routine)
106   USE trabbc_tam      ! bottom boundary condition        (adjoint of tra_bbc routine)
107   USE tradmp_tam      ! internal damping                 (adjoint of tra_dmp routine)
108   USE traadv_tam      ! advection scheme control     (adjoint of tra_adv_ctl routine)
109   USE traldf_tam      ! lateral mixing                   (adjoint of tra_ldf routine)
110   USE cla_tam         ! cross land advection             (adjoint of tra_cla routine)
111   !   zdfkpp          ! KPP non-local tracer fluxes      (adjoint of tra_kpp routine)
112   USE trazdf_tam      ! vertical mixing                  (adjoint of tra_zdf routine)
113   USE tranxt_tam      ! time-stepping                    (adjoint of tra_nxt routine)
114   USE eosbn2_tam      ! equation of state                (adjoint of eos_bn2 routine)
115
116   USE dynadv_tam      ! advection                        (adjoint of dyn_adv routine)
117   USE dynvor_tam      ! vorticity term                   (adjoint of dyn_vor routine)
118   USE dynhpg_tam      ! hydrostatic pressure grad.       (adjoint of dyn_hpg routine)
119   USE dynldf_tam      ! lateral momentum diffusion       (adjoint of dyn_ldf routine)
120   USE dynzdf_tam      ! vertical diffusion               (adjoint of dyn_zdf routine)
121   USE dynspg_tam      ! surface pressure gradient        (adjoint of dyn_spg routine)
122   USE dynnxt_tam      ! time-stepping                    (adjoint of dyn_nxt routine)
123
124!   USE bdy_par_tam
125!   USE bdydta_tam
126
127   USE divcur_tam      ! hor. divergence and curl      (adjoint of div & cur routines)
128   USE cla_div_tam     ! cross land: hor. divergence      (adjoint of div_cla routine)
129   USE wzvmod_tam      ! vertical velocity                (adjoint of wzv     routine)
130
131   USE zdfkpp_tam     ! KPP vertical mixing
132
133   USE zpshde_tam      ! partial step: hor. derivative     (adjoint of zps_hde routine)
134
135!!   USE diaobs_tam      ! obs-minus-model (adjoint of assimilation)   (adjoint of dia_obs routine)
136
137   USE trj_tam
138
139   USE stpctl_tam      ! time stepping control            (adjoint of stp_ctl routine)
140
141   USE gridrandom, ONLY: & 
142                       ! Random Gaussian noise on grids
143      & grid_random, &
144      & grid_rd_sd
145   USE dotprodfld,     ONLY: & ! Computes dot product for 3D and 2D fields
146      & dot_product
147   USE tstool_tam, ONLY: &
148      & prntst_adj,      & 
149      & prntst_tlm,      &
150      & stdemp,          & 
151      & stdu,stdv,       &
152      & stdt,            &   
153      & stds, stdssh,    &
154      & stdr
155
156   USE paresp, ONLY:     & 
157                       ! Weights for an energy-type scalar product
158      & wesp_t,          & 
159      & wesp_s,          & 
160      & wesp_u,          & 
161      & wesp_ssh           
162
163   USE istate_tam      !: Initial state setting          (istate_init routine)
164   USE sol_oce_tam
165   USE trc_oce_tam
166
167   USE step
168
169#if defined key_agrif
170#error 'agrif not yet implemented in nemotam'
171#endif
172
173   IMPLICIT NONE
174   PRIVATE
175
176   PUBLIC stp_tan_cpd
177   !! * Substitutions
178#  include "domzgr_substitute.h90"
179#  include "zdfddm_substitute.h90"
180
181CONTAINS
182#if defined key_agrif
183   #error 'agrif not yet implemented in nemotam'
184  ! SUBROUTINE stp_tan( )
185#else
186   SUBROUTINE stp_tan_cpd( kstp )
187#endif
188      !!----------------------------------------------------------------------
189      !!                     ***  ROUTINE stp_tan  ***
190      !!                     
191      !! ** Purpose of the direct routine:
192      !!              - Time stepping of OPA (momentum and active tracer eqs.)
193      !!              - Time stepping of LIM (dynamic and thermodynamic eqs.)
194      !!              - Tme stepping  of TRC (passive tracer eqs.)
195      !!
196      !! ** Method of the direct routine:
197      !!              -1- Update forcings and data 
198      !!              -2- Update ocean physics
199      !!              -3- Compute the t and s trends
200      !!              -4- Update t and s
201      !!              -5- Compute the momentum trends
202      !!              -6- Update the horizontal velocity
203      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
204      !!              -8- Outputs and diagnostics
205      !!----------------------------------------------------------------------
206      !! * Arguments
207#if defined key_agrif   
208      INTEGER               :: kstp   ! ocean time-step index
209#else
210      INTEGER, INTENT( in ) :: kstp   ! ocean time-step index
211#endif     
212
213      !! * local declarations
214      INTEGER ::   indic    ! error indicator if < 0
215      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zta_tmp, zsa_tmp
216      !! ---------------------------------------------------------------------
217
218
219#if defined key_agrif
220      kstp = nit000 + Agrif_Nb_Step()
221      !      IF ( Agrif_Root() .and. lwp) Write(*,*) '---'
222      !      IF (lwp) Write(*,*) 'Grid Number',Agrif_Fixed(),' time step ',kstp
223#endif   
224      indic = 1                    ! reset to no error condition
225
226      CALL day_tam( kstp, 0 )             ! Calendar
227
228      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
229      ! Update data, open boundaries, surface boundary condition (including sea-ice)
230      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
231
232      ! update 3D temperature data ... not needed in tangent
233
234      ! update 3D salinity data ... not needed in tangent (to be investigated, see sbc_ssr)
235
236      CALL sbc_tan ( kstp ) ! Sea boundary condition (including sea-ice)
237
238      ! update dynamic and tracer data at open boundaries ... not needed for the time being, to investigate whenever we do obc.
239
240      ! compute phase velocities at open boundaries ... not needed for the time being, to investigate whenever we do obc.
241!      IF( lk_bdy     )   CALL bdy_dta_tan( kstp )         ! update dynamic and tracer data at unstructured open boundary
242
243      ! Output the initial state and forcings ... not needed in tangent
244
245      ! saving direct variables ua,va, ta, sa before entering in tracer
246      zta_tmp (:,:,:) = ta (:,:,:)
247      zsa_tmp (:,:,:) = sa (:,:,:)
248      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
249      ! Ocean physics update
250      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
251      !-----------------------------------------------------------------------
252      !  VERTICAL PHYSICS
253      !-----------------------------------------------------------------------
254      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
255      !-----------------------------------------------------------------------
256
257      CALL bn2_tan( tb, sb, tb_tl, sb_tl, rn2_tl )              ! before Brunt-Vaisala frequency
258
259      !                                                     ! Vertical eddy viscosity and diffusivity coefficients
260      ! Richardson number dependent Kz  ... not available
261      ! TKE closure scheme for Kz ... not available
262      ! KPP closure scheme for Kz ... not available
263
264      ! Constant Kz (reset avt, avm[uv] to the background value)... not available
265
266      ! lk_zdfcst_tan:  Constant Kz read from the reference trajectory
267
268
269      ! ! increase diffusivity at rivers mouths... not needed in tangent
270
271      ! enhanced vertical eddy diffusivity ... not needed in tangent with lk_zdfcst_tan
272
273      ! double diffusive mixing ... not needed in tangent with lk_zdfcst_tan
274
275      ! bottom friction... not needed in tangent with lk_zdfcst_tan
276
277      ! mixed layer depth... not needed in tangent with lk_zdfcst_tan
278
279
280      !-----------------------------------------------------------------------
281      !  LATERAL PHYSICS
282      !-----------------------------------------------------------------------
283      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
284      !-----------------------------------------------------------------------
285
286      ! before slope of the lateral mixing... not needed in tangent with lk_zdfcst_tan
287
288#if defined key_traldf_c2d
289      ! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_tan
290#endif
291
292
293#if defined key_top
294      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
295      ! Passive Tracer Model
296      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
297      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
298      !-----------------------------------------------------------------------
299
300      ! time-stepping... not needed in tangent for the time being
301
302#endif
303
304      ta (:,:,:) = zta_tmp (:,:,:)
305      sa (:,:,:) = zsa_tmp (:,:,:)
306      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
307      ! Active tracers
308      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
309      ! N.B. ua, va arrays are used as workspace in this section
310      !-----------------------------------------------------------------------
311
312      ta_tl(:,:,:) = 0.e0            ! set tracer trends to zero
313      sa_tl(:,:,:) = 0.e0
314
315      ! apply tracer assimilation increment ... not needed in tangent
316      ! apply bias directly to tn, sn, tb, sb ... not needed in tangent
317      CALL tra_sbc_tan( kstp )       ! surface boundary condition
318
319      IF( ln_traqsr      )   CALL tra_qsr_tan( kstp )       ! penetrative solar radiation qsr
320
321      IF( lk_trabbc      )   CALL tra_bbc_tan( kstp )       ! bottom heat flux
322
323      ! diffusive bottom boundary layer scheme ... currently not available
324      ! advective (and/or diffusive) bottom boundary layer scheme ... currently not available
325
326      IF( lk_tradmp      )   CALL tra_dmp_tan( kstp )       ! internal damping trends
327
328
329!      CALL tra_adv_tan( kstp )       ! horizontal & vertical advection
330
331!      IF( n_cla == 1     )   CALL tra_cla_tan( kstp )       ! Cross Land Advection (Update Hor. advection)
332
333!!      IF( lk_zdfkpp )        CALL tra_kpp_tan( kstp )       ! KPP non-local tracer fluxes
334
335!      CALL tra_ldf_tan( kstp )       ! lateral mixing
336
337#if defined key_agrif
338      ! tracers sponge ... not available
339#endif
340      CALL tra_zdf_tan( kstp )       ! vertical mixing
341
342
343      ! update the new (t,s) fields by non
344      ! penetrative convective adjustment ... not available
345
346!      CALL tra_nxt_tan( kstp )           ! tracer fields at next time step
347
348
349      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
350         CALL eos_tan( ta, sa, ta_tl, sa_tl, rhd_tl, rhop_tl )      ! Time-filtered in situ density used in dynhpg module
351         IF( ln_zps    )       CALL zps_hde_tan( kstp, ta, sa, ta_tl, sa_tl, rhd_tl,&   ! Partial steps: time filtered hor. gradient
352            &                                     gtu_tl, gsu_tl, gru_tl,   &   ! of t, s, rd at the bottom ocean level
353            &                                     gtv_tl, gsv_tl, grv_tl )
354      ELSE                                                  ! centered hpg (default case)
355         CALL eos_tan( tb, sb, tb_tl, sb_tl, rhd_tl, rhop_tl )       ! now (swap=before) in situ density for dynhpg module
356         IF( ln_zps    )       CALL zps_hde_tan( kstp, tb, sb, tb_tl, sb_tl, rhd_tl,&   ! Partial steps: now horizontal gradient
357            &                                     gtu_tl, gsu_tl, gru_tl,   &   ! of t, s, rd at the bottom ocean level
358            &                                     gtv_tl, gsv_tl, grv_tl )
359      ENDIF
360
361      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
362      ! Dynamics
363      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
364      ! N.B. ta, sa arrays are used as workspace in this section
365      !-----------------------------------------------------------------------
366
367      ua_tl(:,:,:) = 0.e0               ! set dynamics trends to zero
368      va_tl(:,:,:) = 0.e0
369
370      ! apply dynamics assimilation increment ... not needed in tangent
371
372      CALL dyn_adv_tan( kstp )           ! advection (vector or flux form)
373      CALL dyn_vor_tan( kstp )           ! vorticity term including Coriolis
374      CALL dyn_ldf_tan( kstp )           ! lateral mixing
375#if defined key_agrif
376      ! momemtum sponge ... not available
377#endif
378
379      CALL dyn_hpg_tan( kstp )           ! horizontal gradient of Hydrostatic pressure
380
381      CALL dyn_zdf_tan( kstp )           ! vertical diffusion
382
383!PRINT*,'dynspg'
384      IF( lk_dynspg_rl ) THEN
385         ! surface pressure gradient at open boundaries ... not available
386      ENDIF
387      indic=0
388      !i
389      CALL dyn_spg_tan( kstp, indic )    ! surface pressure gradient
390
391      CALL dyn_nxt_tan( kstp )           ! lateral velocity at next time step
392
393      ! vertical mesh at next time step ... not available
394      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
395      ! Computation of diagnostic variables
396      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
397      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
398      !-----------------------------------------------------------------------
399      CALL div_cur_tan( kstp )                 ! Horizontal divergence & Relative vorticity
400      IF( n_cla == 1 ) CALL div_cla_tan( kstp )                 ! Cross Land Advection (Update Hor. divergence)
401      CALL wzv_tan( kstp )                     ! Vertical velocity
402
403      CALL trj_rea( kstp, 1) ! ... Read basic state trajectory at end of current step
404
405      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
406      ! Control, and restarts
407      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
408      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
409      !-----------------------------------------------------------------------
410      !                                                     ! Time loop: control and print
411      CALL stp_ctl_tan( kstp, indic, 0 )
412      IF( indic < 0          )   CALL ctl_stop( 'step_tan: indic < 0' )
413
414      ! close input  ocean restart file ... not needed in tangent
415      ! write output ocean restart file... not needed in tangent
416      ! write open boundary restart file... not needed in tangent
417
418      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
419      ! diagnostics and outputs
420      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
421      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
422      !-----------------------------------------------------------------------
423
424      IF ( nstop == 0 ) THEN                                ! Diagnostics
425         ! drifting Floats... not needed in tangent
426         ! trends: dynamics ... not needed in tangent
427         ! trends: active tracers... not needed in tangent
428         ! trends: Mixed-layer ... not needed in tangent
429         ! trends: vorticity budget... not needed in tangent
430         ! Surface pressure diagnostics... not needed in tangent
431         ! Thermocline depth (20 degres isotherm depth)... not needed in tangent
432         ! basin averaged diagnostics... not needed in tangent
433         ! dynamical heigh diagnostics... not needed in tangent
434         ! Fresh water budget diagnostics... not needed in tangent
435!!!         IF( lk_diaobs  )   CALL dia_obs_tan( kstp )                 ! obs-minus-model (assimilation) diagnostics NOTE: may be better off outside this module
436         ! Poleward TRansports diagnostics... not needed in tangent
437
438
439         !                                                 ! Outputs
440         ! ocean model: outputs... not needed in tangent
441
442      ENDIF
443
444      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
445      ! Coupled mode
446      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
447
448      ! coupled mode : field exchanges ... not available for the next 20 years
449      !
450   END SUBROUTINE stp_tan_cpd
451
452#endif
453END MODULE step_tam_cpd
Note: See TracBrowser for help on using the repository browser.