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/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/step_tam.F90 @ 3658

Last change on this file since 3658 was 3658, checked in by pabouttier, 11 years ago

Missing allocation for sbcssr_tam variables - see Ticket #1022

  • Property svn:executable set to *
File size: 38.1 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   !!                  !  07-12  (P.-A. Bouttier) Phasing with 3.4 version
31   !!----------------------------------------------------------------------
32   !!   stp_tam        : OPA system time-stepping (tangent linear)
33   !!   stp_adj        : OPA system time-stepping (adjoint)
34   !!----------------------------------------------------------------------
35
36   USE step_oce_tam
37#if defined key_agrif
38#error 'agrif not yet implemented in nemotam'
39#endif
40
41   IMPLICIT NONE
42   PRIVATE
43
44   PUBLIC stp_tan,      &
45      &   stp_adj,      & ! called by simvar.F90
46      &   stp_adj_tst
47
48   !! * Substitutions
49#  include "domzgr_substitute.h90"
50#  include "zdfddm_substitute.h90"
51
52CONTAINS
53   SUBROUTINE stp_tan( kstp )
54      !!----------------------------------------------------------------------
55      !!                     ***  ROUTINE stp_tan  ***
56      !!
57      !! ** Purpose of the direct routine:
58      !!              - Time stepping of OPA (momentum and active tracer eqs.)
59      !!              - Time stepping of LIM (dynamic and thermodynamic eqs.)
60      !!              - Tme stepping  of TRC (passive tracer eqs.)
61      !!
62      !! ** Method of the direct routine:
63      !!              -1- Update forcings and data
64      !!              -2- Update ocean physics
65      !!              -3- Compute the t and s trends
66      !!              -4- Update t and s
67      !!              -5- Compute the momentum trends
68      !!              -6- Update the horizontal velocity
69      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
70      !!              -8- Outputs and diagnostics
71      !!----------------------------------------------------------------------
72      !! * Arguments
73      INTEGER, INTENT( in ) :: kstp   ! ocean time-step index
74
75      !! * local declarations
76      INTEGER ::   indic    ! error indicator if < 0
77      !! ---------------------------------------------------------------------
78
79      indic = 0                    ! reset to no error condition
80      IF ( kstp == nit000 )    CALL tl_trj_wri(nit000-1)
81      IF ( kstp /= nit000 )    CALL day_tam( kstp, 0 )             ! Calendar (day was already called at nit000 in day_init)
82
83                               CALL iom_setkt( kstp )                          ! say to iom that we are at time step kstp
84      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
85      ! Update data, open boundaries, surface boundary condition (including sea-ice)
86      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
87
88      ! update 3D temperature data ... not needed in tangent
89
90      ! update 3D salinity data ... not needed in tangent (to be investigated, see sbc_ssr)
91                             CALL sbc_tan ( kstp ) ! Sea boundary condition (including sea-ice)
92
93      ! Output the initial state and forcings ... not needed in tangent
94
95      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
96      !  Ocean dynamics : ssh, wn, hdiv, rot                                 !
97      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
98
99                             CALL ssh_wzv_tan( kstp )         ! after ssh & vertical velocity
100
101      ! saving direct variables ua,va, ta, sa before entering in tracer
102
103      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
104      ! Ocean physics update                (ua, va, ta, sa used as workspace)
105      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
106                             CALL bn2_tan( tsb, tsb_tl, rn2b_tl )        ! now Brunt-Vaisala frequency
107                             CALL bn2_tan( tsn, tsn_tl, rn2_tl )         ! now Brunt-Vaisala frequency
108
109      !-----------------------------------------------------------------------
110      !  VERTICAL PHYSICS
111      !-----------------------------------------------------------------------
112                             CALL zdf_bfr_tan ( kstp ) ! bottom friction...
113
114      !                                                ! Vertical eddy viscosity and diffusivity coefficients
115      ! Richardson number dependent Kz  ... not available
116      ! TKE closure scheme for Kz ... not available
117      ! KPP closure scheme for Kz ... not available
118      ! Constant Kz (reset avt, avm[uv] to the background value)...
119      ! increase diffusivity at rivers mouths... not needed in tangent
120      ! enhanced vertical eddy diffusivity ... not needed in tangent with lk_zdfcst_tan
121      ! double diffusive mixing ... not needed in tangent with lk_zdfcst_tan
122      ! mixed layer depth... not needed in tangent with lk_zdfcst_tan
123      !-----------------------------------------------------------------------
124      !  LATERAL PHYSICS
125      !-----------------------------------------------------------------------
126      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
127      !-----------------------------------------------------------------------
128      ! before slope of the lateral mixing... not needed in tangent with lk_zdfcst_tan
129#if defined key_traldf_c2d
130      ! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_tan
131#endif
132      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
133      ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
134      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
135      !      not needed in tangent
136#if defined key_top
137      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
138      ! Passive Tracer Model
139      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
140      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
141      !-----------------------------------------------------------------------
142
143      ! time-stepping... not needed in tangent for the time being
144
145#endif
146
147      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
148      ! Active tracers                              (ua, va used as workspace)
149      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
150                             tsa_tl(:,:,:,:) = 0.0_wp            ! set tracer trends to zero
151
152      ! apply tracer assimilation increment   ... not needed in tangent
153                             CALL tra_sbc_tan( kstp )       ! surface boundary condition
154      IF( ln_traqsr      )   CALL tra_qsr_tan( kstp )       ! penetrative solar radiation qsr
155      IF( ln_trabbc      )   CALL tra_bbc_tan( kstp )       ! bottom heat flux
156      IF( lk_trabbl  )       CALL tra_bbl_tan( kstp )   ! diffusive bottom boundary layer scheme
157      !! advective (and/or diffusive) bottom boundary layer scheme ... currently not available
158      IF( ln_tradmp      )   CALL tra_dmp_tan( kstp )       ! internal damping trends
159                             CALL tra_adv_tan( kstp )       ! horizontal & vertical advection
160                             CALL tra_ldf_tan( kstp )       ! lateral mixing
161                             CALL tra_zdf_tan( kstp )       ! vertical mixing
162
163      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg (time stepping then eos)
164      ! update the new (t,s) fields by non
165      ! penetrative convective adjustment ... not available
166                             CALL tra_nxt_tan( kstp )                              ! tracer fields at next time step
167                             CALL eos_tan( tsa, tsa_tl,          &                 ! Time-filtered in situ density used in dynhpg module
168                                         & rhd_tl, rhop_tl         )
169         IF( ln_zps    )     CALL zps_hde_tan( kstp, jpts, tsa, tsa_tl, rhd_tl,     &
170            &                                  gtsu_tl, gru_tl, gtsv_tl,    &   ! Partial steps: time filtered hor. gradient
171            &                                  grv_tl )                 ! of t, s, rd at the bottom ocean level
172      ELSE                                              ! centered hpg (default case)
173                             CALL eos_tan( tsn, tsn_tl, &                 ! now (swap=before) in situ density for dynhpg module
174            &                              rhd_tl, rhop_tl )
175
176         IF( ln_zps    )     CALL zps_hde_tan( kstp, jpts, tsn, tsn_tl, rhd_tl,     &
177            &                                  gtsu_tl, gru_tl, gtsv_tl,    &   ! Partial steps: time filtered hor. gradient
178            &                                  grv_tl )
179                             CALL tra_nxt_tan( kstp )                              ! tracer fields at next time step
180      ENDIF
181
182      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
183      ! Dynamics                                    (ta, sa used as workspace)
184      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
185                             ua_tl(:,:,:) = 0.0_wp               ! set dynamics trends to zero
186                             va_tl(:,:,:) = 0.0_wp
187
188      ! apply dynamics assimilation increment ... not needed in tangent
189                             CALL dyn_adv_tan( kstp )       ! advection (vector or flux form)
190                             CALL dyn_vor_tan( kstp )       ! vorticity term including Coriolis
191                             CALL dyn_ldf_tan( kstp )       ! lateral mixing
192                             CALL dyn_hpg_tan( kstp )           ! horizontal gradient of Hydrostatic pressure
193                             CALL dyn_bfr_tan( kstp )           ! bottom friction
194                             CALL dyn_zdf_tan( kstp )           ! vertical diffusion
195                             CALL dyn_spg_tan( kstp, indic )    ! surface pressure gradient
196                             CALL dyn_nxt_tan( kstp )           ! lateral velocity at next time step
197                             CALL ssh_nxt_tan( kstp )           ! sea surface height at next time step
198
199      ! vertical mesh at next time step ... not available
200      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
201      ! Computation of diagnostic variables
202      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
203      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
204      !-----------------------------------------------------------------------
205      CALL tl_trj_wri( kstp )
206      CALL trj_rea( kstp, 1) ! ... Read basic state trajectory at end of current step
207
208      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
209      ! Control, and restarts
210      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
211      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
212      !-----------------------------------------------------------------------
213      !                                                     ! Time loop: control and print
214!*B This fails with cgmod. To be revised      CALL stp_ctl_tan( kstp, indic, 0 )
215      IF( indic < 0          )   CALL ctl_stop( 'step_tan: indic < 0' )
216      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
217      ! diagnostics and outputs
218      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
219      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
220      !-----------------------------------------------------------------------
221      IF ( nstop == 0 ) THEN                                ! Diagnostics
222         ! not needed in tangent
223      ENDIF
224   END SUBROUTINE stp_tan
225
226   SUBROUTINE stp_adj( kstp )
227      !!----------------------------------------------------------------------
228      !!                     ***  ROUTINE stp_adj  ***
229      !!
230      !! ** Purpose of the direct routine:
231      !!              - Time stepping of OPA (momentum and active tracer eqs.)
232      !!              - Time stepping of LIM (dynamic and thermodynamic eqs.)
233      !!              - Tme stepping  of TRC (passive tracer eqs.)
234      !!
235      !! ** Method of the direct routine:
236      !!              -1- Update forcings and data
237      !!              -2- Update ocean physics
238      !!              -3- Compute the t and s trends
239      !!              -4- Update t and s
240      !!              -5- Compute the momentum trends
241      !!              -6- Update the horizontal velocity
242      !!              -7- Compute the diagnostics variables (rd,N2, div,cur,w)
243      !!              -8- Outputs and diagnostics
244      !!----------------------------------------------------------------------
245      !! * Arguments
246      INTEGER, INTENT( in ) :: kstp   ! ocean time-step index
247      !! * local declarations
248      INTEGER ::   indic    ! error indicator if < 0
249      !! ---------------------------------------------------------------------
250
251      indic = 1                    ! reset to no error condition
252                             CALL day_tam( kstp, 1 )            ! Calendar
253
254      ! ... Read basic state trajectory at end of previous step
255                             CALL trj_rea( ( kstp - 1 ), -1 )
256
257      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
258      ! Dynamics                                    (ta, sa used as workspace)
259      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
260      ! apply dynamics assimilation increment ... not needed in adjoint
261                             indic=0
262      ! vertical mesh at next time step ... not available
263                             CALL ssh_nxt_adj( kstp )         ! sea surface height at next time step
264                             CALL dyn_nxt_adj( kstp )           ! lateral velocity at next time step
265                             CALL dyn_spg_adj( kstp, indic )    ! surface pressure gradient
266                             CALL dyn_zdf_adj( kstp )           ! vertical diffusion
267                             CALL dyn_bfr_adj( kstp )           ! bottom friction
268                             CALL dyn_hpg_adj( kstp )           ! horizontal gradient of Hydrostatic pressure
269                             CALL dyn_ldf_adj( kstp )           ! lateral mixing
270                             CALL dyn_vor_adj( kstp )           ! vorticity term including Coriolis
271                             CALL dyn_adv_adj( kstp )           ! advection (vector or flux form)
272
273                             ua_ad(:,:,:)  = 0.0_wp                ! set tracer trends to zero
274                             va_ad(:,:,:)  = 0.0_wp
275
276      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
277      ! Active tracers                              (ua, va used as workspace)
278      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
279
280      IF( ln_dynhpg_imp  ) THEN                             ! semi-implicit hpg
281         IF( ln_zps    )     CALL zps_hde_adj( kstp, jpts, tsa, gtsu, gtsv, tsa_ad,      &   ! Partial steps: time filtered hor. gradient
282            &                                  rhd_ad, gtsu_ad, gru_ad, gtsv_ad,         &   ! of t, s, rd at the bottom ocean level
283            &                                  grv_ad )
284                             CALL eos_adj( tsa, tsa_ad, rhd_ad, rhop_ad )                    ! Time-filtered in situ density used in dynhpg module
285                             CALL tra_nxt_adj( kstp )                                        ! tracer fields at next time step
286      ELSE                                                                                   ! centered hpg (default case)
287                             CALL tra_nxt_adj( kstp )                                        ! tracer fields at next time step
288         IF( ln_zps    )     CALL zps_hde_adj( kstp, jpts, tsn, gtsu, gtsv, tsn_ad,      &   ! Partial steps: time filtered hor. gradient
289            &                                  rhd_ad, gtsu_ad, gru_ad, gtsv_ad,         &   ! of t, s, rd at the bottom ocean level
290            &                                  grv_ad )
291                             CALL eos_adj( tsn, tsn_ad, rhd_ad, rhop_ad )  ! now (swap=before) in situ density for dynhpg module
292      ENDIF
293
294      ! update the new (t,s) fields by non
295      ! penetrative convective adjustment ... not available
296
297                             CALL tra_zdf_adj( kstp )       ! vertical mixing
298                             CALL tra_ldf_adj( kstp )       ! lateral mixing
299                             CALL tra_adv_adj( kstp )      ! horizontal & vertical advection
300      IF( ln_tradmp      )   CALL tra_dmp_adj( kstp )       ! internal damping trends
301      !! advective (and/or diffusive) bottom boundary layer scheme ... currently not available
302      IF( lk_trabbl  )       CALL tra_bbl_adj( kstp )   ! diffusive bottom boundary layer scheme
303      IF( ln_trabbc      )   CALL tra_bbc_adj( kstp )       ! bottom heat flux
304
305      IF( ln_traqsr      )   CALL tra_qsr_adj( kstp )       ! penetrative solar radiation qsr
306
307                             CALL tra_sbc_adj( kstp )       ! surface boundary condition
308
309                             tsa_ad(:,:,:,:) = 0.0_wp            ! set tracer trends to zero
310
311      ! apply tracer assimilation increment ... not needed in adjoint
312#if defined key_top
313      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
314      ! Passive Tracer Model
315      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
316      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
317      !-----------------------------------------------------------------------
318
319      ! time-stepping... not needed in adjoint for the time being
320
321#endif
322      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
323      ! Ocean physics update                (ua, va, ta, sa used as workspace)
324      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
325      !-----------------------------------------------------------------------
326      !  LATERAL PHYSICS
327      !-----------------------------------------------------------------------
328#if defined key_traldf_c2d
329      ! eddy induced velocity coefficient... not needed in tangent with lk_zdfcst_adj
330#endif
331      ! before slope of the lateral mixing... not needed in adjoint with lk_zdfcst_adj
332      !-----------------------------------------------------------------------
333      !  VERTICAL PHYSICS
334      !-----------------------------------------------------------------------
335                             CALL zdf_bfr_adj ( kstp )
336      !! N.B. ua, va, ta, sa arrays are used as workspace in this section
337      !!-----------------------------------------------------------------------
338      !! mixed layer depth... not needed in adjoint with lk_zdfcst_adj
339      !! bottom friction... not needed in adjoint with lk_zdfcst_adj
340      !! double diffusive mixing ... not needed in adjoint with lk_zdfcst_adj
341      !! enhanced vertical eddy diffusivity ... not needed in adjoint with lk_zdfcst_adj
342      !! ! increase diffusivity at rivers mouths... not needed in tangent
343      !! lk_zdfcst_adj:  Constant Kz read from the reference trajectory
344      !! Constant Kz (reset avt, avm[uv] to the background value)... not available
345      !! KPP closure scheme for Kz ... not available
346      !! TKE closure scheme for Kz ... not available
347      !! Richardson number dependent Kz  ... not available
348      !!                                                                 ! Vertical eddy viscosity and diffusivity coefficients
349                             CALL bn2_adj( tsn, tsn_ad, rn2_ad )         ! now Brunt-Vaisala frequency
350                             CALL bn2_adj( tsb, tsb_ad, rn2b_ad )        ! now Brunt-Vaisala frequency
351      !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
352      !!  Ocean dynamics : ssh, wn, hdiv, rot                                 !
353      !!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
354                             CALL ssh_wzv_adj( kstp )         ! after ssh & vertical velocity
355      !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
356      !! Update data, open boundaries and Forcings
357      !!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
358      !! Output the initial state and forcings ... not needed in adjoint
359                             CALL sbc_adj ( kstp ) ! Sea boundary condition (including sea-ice)
360      ! update 3D salinity data ... not needed in tangent
361      ! update 3D temperature data ... not needed in adjoint
362
363      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
364      ! Control, and restarts
365      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
366      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
367      !-----------------------------------------------------------------------
368      !                                                     ! Time loop: control and print
369                             CALL stp_ctl_adj( kstp, indic, 0 )
370      IF( indic < 0      )   CALL ctl_stop( 'step_adj: indic < 0' )
371
372      ! close input  ocean restart file ... not needed in adjoint
373      ! write output ocean restart file... not needed in adjoint
374      ! write open boundary restart file... not needed in adjoint
375      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
376      ! diagnostics and outputs
377      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
378      ! N.B. ua, va, ta, sa arrays are used as workspace in this section
379      !-----------------------------------------------------------------------
380
381      IF ( nstop == 0 ) THEN                                ! Diagnostics
382         ! not needed in adjoint
383      ENDIF
384
385      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
386      ! Coupled mode
387      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
388
389      ! coupled mode : field exchanges ... not available for the next 20 years
390      !
391      !
392      IF( nn_timing == 1 .AND.  kstp == nit000  )   CALL timing_reset
393      !
394   END SUBROUTINE stp_adj
395
396   SUBROUTINE stp_adj_tst( kumadt )
397      !!-----------------------------------------------------------------------
398      !!
399      !!                  ***  ROUTINE example_adj_tst ***
400      !!
401      !! ** Purpose : Test the adjoint routine.
402      !!
403      !! ** Method  : Verify the scalar product
404      !!
405      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
406      !!
407      !!              where  L   = tangent routine
408      !!                     L^T = adjoint routine
409      !!                     W   = diagonal matrix of scale factors
410      !!                     dx  = input perturbation (random field)
411      !!                     dy  = L dx
412      !!
413      !!
414      !! History :
415      !!        ! 09-03 (F. Vigilant)
416      !!-----------------------------------------------------------------------
417      !! * Modules used
418      USE sbcssr_tam, ONLY: &
419         & qrp_tl, qrp_ad,  &
420         & erp_tl, erp_ad
421      USE sbcfwb_tam, ONLY: &
422         & a_fwb_tl,        & ! for before year.
423         & a_fwb_ad  ! for before year.
424
425      !! * Arguments
426      INTEGER, INTENT(IN) :: &
427         & kumadt               ! Output unit
428
429      !! * Local declarations
430      INTEGER ::             &
431         & ji,               &  ! dummy loop indices
432         & jj,               &
433         & jk,               &
434         & istp
435      INTEGER, DIMENSION(jpi,jpj) :: &
436         & iseed_2d             ! 2D seed for the random number generator
437      REAL(KIND=wp) ::       &
438         & zsp1    ,         &  ! scalar product involving the tangent routine
439         & zsp1_U  ,         &
440         & zsp1_V  ,         &
441         & zsp1_T  ,         &
442         & zsp1_S  ,         &
443         & zsp1_SSH,         &
444         & zsp2    ,         &  ! scalar product involving the adjoint routine
445         & zsp2_U  ,         &
446         & zsp2_V  ,         &
447         & zsp2_T  ,         &
448         & zsp2_S  ,         &
449         & zsp2_SSH
450      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
451         & zun_tlin ,        & ! Tangent input
452         & zvn_tlin ,        & ! Tangent input
453         & ztn_tlin ,        & ! Tangent input
454         & zsn_tlin ,        & ! Tangent input
455         & zun_tlout,        & ! Tangent output
456         & zvn_tlout,        & ! Tangent output
457         & ztn_tlout,        & ! Tangent output
458         & zsn_tlout,        & ! Tangent outpu
459         & zun_adin ,        & ! Adjoint input
460         & zvn_adin ,        & ! Adjoint input
461         & ztn_adin ,        & ! Adjoint input
462         & zsn_adin ,        & ! Adjoint input
463         & zun_adout,        & ! Adjoint output
464         & zvn_adout,        & ! Adjoint output
465         & ztn_adout,        & ! Adjoint output
466         & zsn_adout,        & ! Adjoint output
467         & z3r                 ! 3D random field
468      REAL(KIND=wp), DIMENSION(:,:), ALLOCATABLE :: &
469         & zsshn_tlin  ,     & ! Tangent input
470         & zsshn_tlout ,     & ! Tangent input
471         & zsshn_adin  ,     & ! Adjoint output
472         & zsshn_adout ,     & ! Adjoint output
473         & z2r                 ! 2D random field
474
475      CHARACTER(LEN=14) ::   &
476         & cl_name
477
478      INTEGER ::             &
479         & jpert
480      INTEGER, PARAMETER :: jpertmax = 6
481
482      ! Allocate memory
483
484      ALLOCATE( &
485         & zun_tlin(   jpi,jpj,jpk), zvn_tlin(   jpi,jpj,jpk), ztn_tlin(  jpi,jpj,jpk), &
486         & zsn_tlin(   jpi,jpj,jpk), zsshn_tlin( jpi,jpj    ), zun_tlout( jpi,jpj,jpk), &
487         & zvn_tlout(  jpi,jpj,jpk), ztn_tlout(  jpi,jpj,jpk), zsn_tlout( jpi,jpj,jpk), &
488         & zsshn_tlout(jpi,jpj    ), zun_adin(   jpi,jpj,jpk), zvn_adin(  jpi,jpj,jpk), &
489         & ztn_adin(   jpi,jpj,jpk), zsn_adin(   jpi,jpj,jpk), zsshn_adin(jpi,jpj    ), &
490         & zun_adout(  jpi,jpj,jpk), zvn_adout(  jpi,jpj,jpk), ztn_adout( jpi,jpj,jpk), &
491         & zsn_adout(  jpi,jpj,jpk), zsshn_adout(jpi,jpj    ), z2r(       jpi,jpj    ), &
492         & z3r(        jpi,jpj,jpk)                                                     &
493         & )
494
495      !==================================================================
496      ! 1) dx = ( un_tl, vn_tl, tn_tl, sn_tl, sshn_tl ) and
497      !    dy = ( ..... )
498      !==================================================================
499
500
501      !DO jpert = 1, jpertmax
502      jpert = 6
503         !--------------------------------------------------------------------
504         ! Reset the tangent and adjoint variables
505         !--------------------------------------------------------------------
506         zun_tlin   (:,:,:) = 0.0_wp
507         zvn_tlin   (:,:,:) = 0.0_wp
508         ztn_tlin   (:,:,:) = 0.0_wp
509         zsn_tlin   (:,:,:) = 0.0_wp
510         zsshn_tlin (  :,:) = 0.0_wp
511         zun_tlout  (:,:,:) = 0.0_wp
512         zvn_tlout  (:,:,:) = 0.0_wp
513         ztn_tlout  (:,:,:) = 0.0_wp
514         zsn_tlout  (:,:,:) = 0.0_wp
515         zsshn_tlout(  :,:) = 0.0_wp
516         zun_adin   (:,:,:) = 0.0_wp
517         zvn_adin   (:,:,:) = 0.0_wp
518         ztn_adin   (:,:,:) = 0.0_wp
519         zsn_adin   (:,:,:) = 0.0_wp
520         zsshn_adin (  :,:) = 0.0_wp
521         zun_adout  (:,:,:) = 0.0_wp
522         zvn_adout  (:,:,:) = 0.0_wp
523         ztn_adout  (:,:,:) = 0.0_wp
524         zsn_adout  (:,:,:) = 0.0_wp
525         zsshn_adout(  :,:) = 0.0_wp
526         z3r        (:,:,:) = 0.0_wp
527         z2r        (  :,:) = 0.0_wp
528
529         !--------------------------------------------------------------------
530         ! Initialize the tangent input with random noise: dx
531         !--------------------------------------------------------------------
532
533         IF ( (jpert == 1) .OR. (jpert == jpertmax) ) THEN
534            CALL grid_random(  z3r, 'U', 0.0_wp, stdu )
535            DO jk = 1, jpk
536               DO jj = nldj, nlej
537                  DO ji = nldi, nlei
538                     zun_tlin(ji,jj,jk) = z3r(ji,jj,jk)
539                  END DO
540               END DO
541            END DO
542         ENDIF
543         IF ( (jpert == 2) .OR. (jpert == jpertmax) ) THEN
544            CALL grid_random(  z3r, 'V', 0.0_wp, stdv )
545            DO jk = 1, jpk
546               DO jj = nldj, nlej
547                  DO ji = nldi, nlei
548                     zvn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
549                  END DO
550               END DO
551            END DO
552         ENDIF
553         IF ( (jpert == 3) .OR. (jpert == jpertmax) ) THEN
554            CALL grid_random(  z3r, 'T', 0.0_wp, stdt )
555            DO jk = 1, jpk
556               DO jj = nldj, nlej
557                  DO ji = nldi, nlei
558                     ztn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
559                  END DO
560               END DO
561            END DO
562         ENDIF
563         IF ( (jpert == 4) .OR. (jpert == jpertmax) ) THEN
564            CALL grid_random(  z3r, 'T', 0.0_wp, stds )
565            DO jk = 1, jpk
566               DO jj = nldj, nlej
567                  DO ji = nldi, nlei
568                     zsn_tlin(ji,jj,jk) = z3r(ji,jj,jk)
569                  END DO
570               END DO
571            END DO
572         ENDIF
573         IF ( (jpert == 5) .OR. (jpert == jpertmax) ) THEN
574            CALL grid_random(  z2r, 'T', 0.0_wp, stdssh )
575            DO jj = nldj, nlej
576               DO ji = nldi, nlei
577                  zsshn_tlin(ji,jj) = z2r(ji,jj)
578               END DO
579            END DO
580         ENDIF
581
582         IF (jpert>0) THEN
583
584            !CALL     oce_tam_init( 1 )    ! allocate/initialize tl variables
585            !CALL sbc_oce_tam_init( 1 )
586            !CALL sol_oce_tam_init( 1 )
587            !CALL trc_oce_tam_init( 1 )
588
589            qrp_tl = 0.0_wp
590            erp_tl = 0.0_wp
591
592            emp_tl(:,:) = 0.0_wp
593#if defined key_tradmp
594            strdmp_tl = 0.0_wp
595            ttrdmp_tl = 0.0_wp
596#endif
597            a_fwb_tl = 0.0_wp
598
599            istp = nit000 - 1
600            CALL  trj_rea( istp, 1 )
601            !--------------------------------------------------------------------
602            ! Initialize the tangent variables: dy^* = W dy
603            !--------------------------------------------------------------------
604
605            un_tl  (:,:,:) = zun_tlin   (:,:,:)
606            vn_tl  (:,:,:) = zvn_tlin   (:,:,:)
607            tsn_tl  (:,:,:,jp_tem) = ztn_tlin   (:,:,:)
608            tsn_tl  (:,:,:,jp_sal) = zsn_tlin   (:,:,:)
609            sshn_tl(  :,:) = zsshn_tlin (  :,:)
610
611            !CALL     oce_tam_deallocate( 2 )    ! deallocate adj variables
612            !CALL sbc_oce_tam_deallocate( 2 )
613            !CALL sol_oce_tam_deallocate( 2 )
614
615            !-----------------------------------------------------------------------
616            !  Initialization of the dynamics and tracer fields for the tangent
617            !-----------------------------------------------------------------------
618
619            CALL istate_init_tan
620
621            DO istp = nit000, nitend, 1
622               CALL stp_tan( istp )
623            END DO
624
625            zun_tlout  ( :,:,:) = un_tl   (:,:,:)
626            zvn_tlout  ( :,:,:) = vn_tl   (:,:,:)
627            ztn_tlout  ( :,:,:) = tsn_tl   (:,:,:,jp_tem)
628            zsn_tlout  ( :,:,:) = tsn_tl   (:,:,:,jp_sal)
629            zsshn_tlout(   :,:) = sshn_tl (  :,:)
630
631            !--------------------------------------------------------------------
632            ! Initialize the adjoint variables: dy^* = W dy
633            !--------------------------------------------------------------------
634
635            DO jk = 1, jpk
636               DO jj = nldj, nlej
637                  DO ji = nldi, nlei
638                     zun_adin(ji,jj,jk) = zun_tlout(ji,jj,jk) &
639                          &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) &
640                          &               * umask(ji,jj,jk) * wesp_u
641                  END DO
642               END DO
643            END DO
644
645            DO jk = 1, jpk
646               DO jj = nldj, nlej
647                  DO ji = nldi, nlei
648                     zvn_adin(ji,jj,jk) = zvn_tlout(ji,jj,jk) &
649                          &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) &
650                          &               * vmask(ji,jj,jk) * wesp_u
651                  END DO
652               END DO
653            END DO
654
655            DO jk = 1, jpk
656               DO jj = nldj, nlej
657                  DO ji = nldi, nlei
658                     ztn_adin(ji,jj,jk) = ztn_tlout(ji,jj,jk) &
659                          &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
660                          &               * tmask(ji,jj,jk) * wesp_t(jk)
661                  END DO
662               END DO
663            END DO
664
665            DO jk = 1, jpk
666               DO jj = nldj, nlej
667                  DO ji = nldi, nlei
668                     zsn_adin(ji,jj,jk) = zsn_tlout(ji,jj,jk) &
669                          &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
670                          &               * tmask(ji,jj,jk) * wesp_s(jk)
671                  END DO
672               END DO
673            END DO
674
675            DO jj = nldj, nlej
676               DO ji = nldi, nlei
677                  zsshn_adin(ji,jj) = zsshn_tlout(ji,jj) &
678                       &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,1) &
679                       &               * tmask(ji,jj,1) * wesp_ssh
680               END DO
681            END DO
682
683            !--------------------------------------------------------------------
684            ! Compute the scalar product: ( L dx )^T W dy
685            !--------------------------------------------------------------------
686
687            zsp1_U    = DOT_PRODUCT( zun_tlout  , zun_adin    )
688            zsp1_V    = DOT_PRODUCT( zvn_tlout  , zvn_adin    )
689            zsp1_T    = DOT_PRODUCT( ztn_tlout  , ztn_adin    )
690            zsp1_S    = DOT_PRODUCT( zsn_tlout  , zsn_adin    )
691            zsp1_SSH  = DOT_PRODUCT( zsshn_tlout, zsshn_adin  )
692
693            zsp1      = zsp1_U + zsp1_V + zsp1_T + zsp1_S + zsp1_SSH
694
695            !--------------------------------------------------------------------
696            ! Call the adjoint routine: dx^* = L^T dy^*
697            !--------------------------------------------------------------------
698            !CALL     oce_tam_init( 2 )    ! allocate/initialize adj variables
699            !CALL sbc_oce_tam_init( 2 )
700            !CALL sol_oce_tam_init( 2 )
701            !CALL trc_oce_tam_init( 2 )
702
703            qrp_ad = 0.0_wp
704            erp_ad = 0.0_wp
705            emp_ad(:,:) = 0.0_wp
706#if defined key_tradmp
707            strdmp_ad = 0.0_wp
708            ttrdmp_ad = 0.0_wp
709#endif
710            a_fwb_ad = 0.0_wp
711
712            un_ad  (:,:,:) = zun_adin   (:,:,:)
713            vn_ad  (:,:,:) = zvn_adin   (:,:,:)
714            tsn_ad  (:,:,:,jp_tem) = ztn_adin   (:,:,:)
715            tsn_ad  (:,:,:,jp_sal) = zsn_adin   (:,:,:)
716            sshn_ad(  :,:) = zsshn_adin (  :,:)
717
718            !CALL     oce_tam_deallocate( 1 )    !deallocate tl variables
719            !CALL sbc_oce_tam_deallocate( 1 )
720            !CALL sol_oce_tam_deallocate( 1 )
721            !CALL trc_oce_tam_deallocate( 1 )
722
723            DO istp = nitend, nit000, -1
724               CALL stp_adj(istp)
725            END DO
726
727            CALL istate_init_adj
728
729            zun_adout  ( :,:,:) = un_ad   (:,:,:)
730            zvn_adout  ( :,:,:) = vn_ad   (:,:,:)
731            ztn_adout  ( :,:,:) = tsn_ad   (:,:,:,jp_tem)
732            zsn_adout  ( :,:,:) = tsn_ad   (:,:,:,jp_sal)
733            zsshn_adout(   :,:) = sshn_ad (  :,:)
734
735            zsp2_U    = DOT_PRODUCT( zun_tlin  , zun_adout )
736            zsp2_V    = DOT_PRODUCT( zvn_tlin  , zvn_adout )
737            zsp2_T    = DOT_PRODUCT( ztn_tlin  , ztn_adout )
738            zsp2_S    = DOT_PRODUCT( zsn_tlin  , zsn_adout )
739            zsp2_SSH  = DOT_PRODUCT( zsshn_tlin, zsshn_adout )
740
741            zsp2      = zsp2_U + zsp2_V + zsp2_T + zsp2_S + zsp2_SSH
742
743            !CALL     oce_tam_deallocate( 2 )    !deallocate adj variables
744            !CALL sbc_oce_tam_deallocate( 2 )
745            !CALL sol_oce_tam_deallocate( 2 )
746            !CALL trc_oce_tam_deallocate( 2 )
747
748            ! 14 char:'12345678901234'
749            SELECT CASE (jpert)
750            CASE(1)
751               cl_name = 'step_adj     U'
752            CASE(2)
753               cl_name = 'step_adj     V'
754            CASE(3)
755               cl_name = 'step_adj     T'
756            CASE(4)
757               cl_name = 'step_adj     S'
758            CASE(5)
759               cl_name = 'step_adj   SSH'
760            CASE(jpertmax)
761               cl_name = 'step_adj      '
762            END SELECT
763            CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
764         endif
765      !END DO
766
767      DEALLOCATE(                                 &
768         & zun_tlin   , zvn_tlin   , ztn_tlin   , &
769         & zsn_tlin   , zsshn_tlin , zun_tlout  , &
770         & zvn_tlout  , ztn_tlout  , zsn_tlout  , &
771         & zsshn_tlout, zun_adin   , zvn_adin   , &
772         & ztn_adin   , zsn_adin   , zsshn_adin , &
773         & zun_adout  , zvn_adout  , ztn_adout  , &
774         & zsn_adout  , zsshn_adout, z2r        , &
775         & z3r                                    &
776         & )
777
778   END SUBROUTINE stp_adj_tst
779
780      !!----------------------------------------------------------------------
781
782   !!======================================================================
783#endif
784END MODULE step_tam
Note: See TracBrowser for help on using the repository browser.