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.
tranxt_tam.F90 in branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA – NEMO

source: branches/2012/dev_r3604_LEGI8_TAM/NEMOGCM/NEMO/OPATAM_SRC/TRA/tranxt_tam.F90 @ 3611

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

Add TAM code and ORCA2_TAM configuration - see Ticket #1007
  • Property svn:executable set to *
File size: 34.3 KB
Line 
1MODULE tranxt_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  tranxt_tam  ***
5   !! Ocean active tracers:  time stepping on temperature and salinity
6   !!               Tangent and Adjoint module
7   !!======================================================================
8   !! History of the direct module:
9   !!            7.0  !  91-11  (G. Madec)  Original code
10   !!                 !  93-03  (M. Guyon)  symetrical conditions
11   !!                 !  96-02  (G. Madec & M. Imbard)  opa release 8.0
12   !!            8.0  !  96-04  (A. Weaver)  Euler forward step
13   !!            8.2  !  99-02  (G. Madec, N. Grima)  semi-implicit pressure grad.
14   !!  NEMO      1.0  !  2002-08  (G. Madec)  F90: Free form and module
15   !!             -   !  2002-11  (C. Talandier, A-M Treguier) Open boundaries
16   !!             -   !  2005-04  (C. Deltel) Add Asselin trend in the ML budget
17   !!            2.0  !  2006-02  (L. Debreu, C. Mazauric) Agrif implementation
18   !!            3.0  !  2008-06  (G. Madec)  time stepping always done in trazd
19   !!            3.1  !  2009-02  (G. Madec, R. Benshila)  re-introduce the vvl option
20   !! History of the TAM module:
21   !!            2.0  !  2008-09  (A. Vidard) tam of the 2006-02 version
22   !!            3.0  !  2008-11  (A. Vidard) tam of the 2008-06 version
23   !!             -   !  2009-01  (A. Weaver) corrections to test
24   !!            3.2  !  2010-04  (F. Vigilant)  version 3.2
25   !!----------------------------------------------------------------------
26
27   !!----------------------------------------------------------------------
28   !!   tra_nxt_tan    : time stepping on temperature and salinity (tangent)
29   !!   tra_nxt_adj    : time stepping on temperature and salinity (adjoint)
30   !!----------------------------------------------------------------------
31   USE par_kind
32   USE par_oce
33   USE dynhpg
34   USE oce_tam
35   USE zdf_oce
36   USE dom_oce
37   USE tranxt
38   USE in_out_manager
39   USE lbclnk
40   USE lbclnk_tam
41   USE gridrandom
42   USE dotprodfld
43   USE paresp
44   USE tstool_tam                   !           salinity
45   USE traqsr
46   USE wrk_nemo
47   USE timing
48
49   IMPLICIT NONE
50   PRIVATE
51
52   !! * Routine accessibility
53   PUBLIC   tra_nxt_tan     ! routine called by step_tam.F90
54   PUBLIC   tra_nxt_adj     ! routine called by step_tam.F90
55   PUBLIC   tra_nxt_adj_tst ! routine called by tst.F90
56
57   REAL(wp) :: rbcp
58   !! * Substitutions
59#  include "domzgr_substitute.h90"
60
61CONTAINS
62
63   SUBROUTINE tra_nxt_tan( kt )
64      !!----------------------------------------------------------------------
65      !!                   ***  ROUTINE tranxt_tan  ***
66      !!
67      !! ** Purpose of the direct routine:
68      !!             Apply the boundary condition on the after temperature
69      !!             and salinity fields, achieved the time stepping by adding
70      !!             the Asselin filter on now fields and swapping the fields.
71      !!
72      !! ** Method  :   At this stage of the computation, ta and sa are the
73      !!             after temperature and salinity as the time stepping has
74      !!             been performed in trazdf_imp or trazdf_exp module.
75      !!
76      !!              - Apply lateral boundary conditions on (ta,sa)
77      !!             at the local domain   boundaries through lbc_lnk call,
78      !!             at the radiative open boundaries (lk_obc=T),
79      !!             at the relaxed   open boundaries (lk_bdy=T), and
80      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
81      !!
82      !!              - Update lateral boundary conditions on AGRIF children
83      !!             domains (lk_agrif=T)
84      !!
85      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
86      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
87      !!----------------------------------------------------------------------
88      !!
89      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
90      INTEGER             :: jn, jk
91      !!----------------------------------------------------------------------
92      !
93      IF( nn_timing == 1 )  CALL timing_start('tra_nxt_tan')
94      !
95      IF( kt == nit000 ) THEN
96         IF(lwp) WRITE(numout,*)
97         IF(lwp) WRITE(numout,*) 'tra_nxt_tan : achieve the time stepping by Asselin filter and array swap'
98         IF(lwp) WRITE(numout,*) '~~~~~~~'
99      ENDIF
100      rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)
101
102      ! Update after tracer on domain lateral boundaries
103      !
104      CALL lbc_lnk( tsa_tl(:,:,:,jp_tem), 'T', 1.0_wp )      ! local domain boundaries  (T-point, unchanged sign)
105      CALL lbc_lnk( tsa_tl(:,:,:,jp_sal), 'T', 1.0_wp )
106      !
107      ! set time step size (Euler/Leapfrog)
108      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
109      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
110      ENDIF
111
112      IF ( neuler == 0 .AND. kt == nit000 ) THEN
113         DO jn = 1, jpts
114            DO jk = 1, jpkm1
115               tsn_tl(:,:,jk,jn) = tsa_tl(:,:,jk,jn)
116            END DO
117         END DO
118      ELSE
119         ! Leap-Frog + Asselin filter time stepping
120         IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl_tan( kt )      ! variable volume level (vvl)
121         ELSE                  ;   CALL tra_nxt_fix_tan( kt, nit000, 'TRA', tsb_tl, tsn_tl, tsa_tl, jpts )      ! fixed    volume level
122         ENDIF
123      END IF
124      !
125      IF( nn_timing == 1 )  CALL timing_stop('tra_nxt_tan')
126      !
127   END SUBROUTINE tra_nxt_tan
128
129   SUBROUTINE tra_nxt_fix_tan( kt, kit000, cdtype, ptb_tl, ptn_tl, pta_tl, kjpt )
130      !!----------------------------------------------------------------------
131      !!                   ***  ROUTINE tra_nxt_fix_tan  ***
132      !!
133      !! ** Purpose :   fixed volume: apply the Asselin time filter and
134      !!                swap the tracer fields.
135      !!
136      !! ** Method  : - Apply a Asselin time filter on now fields.
137      !!              - save in (ta,sa) an average over the three time levels
138      !!             which will be used to compute rdn and thus the semi-implicit
139      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
140      !!              - swap tracer fields to prepare the next time_step.
141      !!                This can be summurized for tempearture as:
142      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
143      !!             ztm = 0                   otherwise
144      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
145      !!             tn  = ta
146      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
147      !!
148      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
149      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
150      !!----------------------------------------------------------------------
151      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
152      INTEGER         , INTENT(in   )                               ::   kit000      ! first time step index
153      CHARACTER(len=3), INTENT(in   )                               ::   cdtype      ! =TRA or TRC (tracer indicator)
154      INTEGER         , INTENT(in   )                               ::   kjpt        ! number of tracers
155      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb_tl      ! before tracer fields
156      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn_tl      ! now tracer fields
157      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta_tl      ! tracer trend
158      !!
159      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
160      REAL(wp) ::   ztntl, ztdtl, ztn     ! temporary scalars
161      LOGICAL  ::   ll_tra_hpg
162      !!----------------------------------------------------------------------
163
164      IF( kt == kit000 ) THEN
165         IF(lwp) WRITE(numout,*)
166         IF(lwp) WRITE(numout,*) 'tra_nxt_fix_tan : time stepping ', cdtype
167         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
168      ENDIF
169      ztntl = 0._wp
170      ztdtl = 0._wp
171      !
172      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
173      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
174      ENDIF
175      !
176      DO jn = 1, kjpt
177         !
178         DO jk = 1, jpkm1
179            DO jj = 1, jpj
180               DO ji = 1, jpi
181                  ztntl = ptn_tl(ji,jj,jk,jn)
182                  ztdtl = pta_tl(ji,jj,jk,jn) - 2. * ztntl + ptb_tl(ji,jj,jk,jn)      !  time laplacian on tracers
183                  !
184                  ptb_tl(ji,jj,jk,jn) = ztntl + atfp * ztdtl                          ! ptb <-- filtered ptn
185                  ptn_tl(ji,jj,jk,jn) = pta_tl(ji,jj,jk,jn)                           ! ptn <-- pta
186                  !
187                  IF( ll_tra_hpg )   pta_tl(ji,jj,jk,jn) = ztntl + rbcp * ztdtl       ! pta <-- Brown & Campana average
188               END DO
189           END DO
190         END DO
191         !
192      END DO
193      !
194   END SUBROUTINE tra_nxt_fix_tan
195
196   SUBROUTINE tra_nxt_vvl_tan( kt )
197      !!----------------------------------------------------------------------
198      !!                   ***  ROUTINE tra_nxt_vvl_tan  ***
199      !!
200      !! ** Purpose :   Time varying volume: apply the Asselin time filter
201      !!                and swap the tracer fields.
202      !!
203      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
204      !!              - save in (ta,sa) a thickness weighted average over the three
205      !!             time levels which will be used to compute rdn and thus the semi-
206      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
207      !!              - swap tracer fields to prepare the next time_step.
208      !!                This can be summurized for tempearture as:
209      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T
210      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )
211      !!             ztm = 0                                otherwise
212      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
213      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
214      !!             tn  = ta
215      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
216      !!
217      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
218      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
219      !!----------------------------------------------------------------------
220      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
221      !!
222      INTEGER  ::   ji, jj, jk             ! dummy loop indices
223      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
224      REAL(wp) ::   zsm , zsc_f , zsf , zsca, zscn, zscb   !    -         -
225      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
226      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
227      !!----------------------------------------------------------------------
228
229      IF( kt == nit000 ) THEN
230         IF(lwp) WRITE(numout,*)
231         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl_tan : time stepping'
232         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
233      ENDIF
234
235      IF(lwp) WRITE(numout,*) "key_vvl net available in tangent yet"
236      CALL abort
237      !
238   END SUBROUTINE tra_nxt_vvl_tan
239
240   SUBROUTINE tra_nxt_adj( kt )
241      !!----------------------------------------------------------------------
242      !!                   ***  ROUTINE tranxt_adj  ***
243      !!
244      !! ** Purpose of the direct routine:
245      !!             Apply the boundary condition on the after temperature
246      !!             and salinity fields, achieved the time stepping by adding
247      !!             the Asselin filter on now fields and swapping the fields.
248      !!
249      !! ** Method  :   At this stage of the computation, ta and sa are the
250      !!             after temperature and salinity as the time stepping has
251      !!             been performed in trazdf_imp or trazdf_exp module.
252      !!
253      !!              - Apply lateral boundary conditions on (ta,sa)
254      !!             at the local domain   boundaries through lbc_lnk call,
255      !!             at the radiative open boundaries (lk_obc=T),
256      !!             at the relaxed   open boundaries (lk_bdy=T), and
257      !!             at the AGRIF zoom     boundaries (lk_agrif=T)
258      !!
259      !!              - Update lateral boundary conditions on AGRIF children
260      !!             domains (lk_agrif=T)
261      !!
262      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
263      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
264      !!----------------------------------------------------------------------
265      !!
266      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
267      INTEGER             ::   jn, jk
268      !!----------------------------------------------------------------------
269      !
270      IF( nn_timing == 1 )  CALL timing_start( 'tra_nxt_adj')
271      !
272      IF( kt == nitend ) THEN
273         IF(lwp) WRITE(numout,*)
274         IF(lwp) WRITE(numout,*) 'tra_nxt_adj : achieve the time stepping by Asselin filter and array swap'
275         IF(lwp) WRITE(numout,*) '~~~~~~~'
276         rbcp = 0.25 * (1. + atfp) * (1. + atfp) * ( 1. - atfp)
277      ENDIF
278
279      ! set time step size (Euler/Leapfrog)
280      r2dtra(:) =  2.* rdttra(:) ! initialization
281      IF( neuler == 0 .AND. kt == nit000 ) THEN   ;   r2dtra(:) =     rdttra(:)      ! at nit000             (Euler)
282      ELSEIF( kt <= nit000 + 1 )           THEN   ;   r2dtra(:) = 2.* rdttra(:)      ! at nit000 or nit000+1 (Leapfrog)
283      ENDIF
284
285
286      IF( neuler == 0 .AND. kt == nit000 ) THEN
287         DO jn = 1, jpts
288            DO jk = jpkm1, 1, -1
289               tsa_ad(:,:,jk,jn) = tsa_ad(:,:,jk,jn) + tsn_ad(:,:,jk,jn)
290               tsn_ad(:,:,jk,jn) = 0._wp
291            END DO
292         END DO
293      ELSE
294         !! Leap-Frog + Asselin filter time stepping
295         IF( lk_vvl )   THEN   ;   CALL tra_nxt_vvl_adj( kt )      ! variable volume level (vvl)
296         ELSE                  ;   CALL tra_nxt_fix_adj( kt, nit000, 'TRA', tsb_ad, tsn_ad, tsa_ad, jpts )      ! fixed    volume level
297         ENDIF
298      ENDIF
299
300      ! Update after tracer on domain lateral boundaries
301      !
302      CALL lbc_lnk_adj( tsa_ad(:,:,:,jp_sal), 'T', 1.0_wp )
303      CALL lbc_lnk_adj( tsa_ad(:,:,:,jp_tem), 'T', 1.0_wp )      ! local domain boundaries  (T-point, unchanged sign)
304      !
305      !
306      IF( nn_timing == 1 )  CALL timing_stop('tra_nxt_adj')
307      !
308   END SUBROUTINE tra_nxt_adj
309
310   SUBROUTINE tra_nxt_fix_adj( kt, kit000, cdtype, ptb_ad, ptn_ad, pta_ad, kjpt )
311      !!----------------------------------------------------------------------
312      !!                   ***  ROUTINE tra_nxt_fix_adj  ***
313      !!
314      !! ** Purpose :   fixed volume: apply the Asselin time filter and
315      !!                swap the tracer fields.
316      !!
317      !! ** Method  : - Apply a Asselin time filter on now fields.
318      !!              - save in (ta,sa) an average over the three time levels
319      !!             which will be used to compute rdn and thus the semi-implicit
320      !!             hydrostatic pressure gradient (ln_dynhpg_imp = T)
321      !!              - swap tracer fields to prepare the next time_step.
322      !!                This can be summurized for tempearture as:
323      !!             ztm = (ta+2tn+tb)/4       ln_dynhpg_imp = T
324      !!             ztm = 0                   otherwise
325      !!             tb  = tn + atfp*[ tb - 2 tn + ta ]
326      !!             tn  = ta
327      !!             ta  = ztm       (NB: reset to 0 after eos_bn2 call)
328      !!
329      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
330      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
331      !!----------------------------------------------------------------------
332      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
333      INTEGER         , INTENT(in   )                               ::   kit000      ! first time step index
334      CHARACTER(len=3), INTENT(in   )                               ::   cdtype      ! =TRA or TRC (tracer indicator)
335      INTEGER         , INTENT(in   )                               ::   kjpt        ! number of tracers
336      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptb_ad      ! before tracer fields
337      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   ptn_ad      ! now tracer fields
338      REAL(wp)        , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt)  ::   pta_ad      ! tracer trend
339      !!
340      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
341      REAL(wp) ::   ztnad, ztdad, ztn     ! temporary scalars
342      LOGICAL  ::   ll_tra_hpg
343      !!----------------------------------------------------------------------
344
345      IF( kt == kit000 ) THEN
346         IF(lwp) WRITE(numout,*)
347         IF(lwp) WRITE(numout,*) 'tra_nxt_fix_adj : time stepping ', cdtype
348         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
349      ENDIF
350      !
351      IF( cdtype == 'TRA' )  THEN   ;   ll_tra_hpg = ln_dynhpg_imp    ! active  tracers case  and  semi-implicit hpg
352      ELSE                          ;   ll_tra_hpg = .FALSE.          ! passive tracers case or NO semi-implicit hpg
353      ENDIF
354      ztnad = 0._wp
355      ztdad = 0._wp
356      !
357      DO jn = 1, kjpt
358         !
359         DO jk = jpkm1, 1, -1
360            DO jj = jpj, 1, -1
361               DO ji = jpi, 1, -1
362                  IF( ll_tra_hpg ) THEN
363                     ztnad = ztnad + pta_ad(ji,jj,jk,jn)
364                     ztdad = ztdad + rbcp * pta_ad(ji,jj,jk,jn)
365                     pta_ad(ji,jj,jk,jn) = 0._wp
366                  END IF
367                  pta_ad(ji,jj,jk,jn) = pta_ad(ji,jj,jk,jn) + ptn_ad(ji,jj,jk,jn)
368                  ptn_ad(ji,jj,jk,jn) = 0._wp
369                  ztdad = ztdad + atfp * ptb_ad(ji,jj,jk,jn)
370                  ztnad = ztnad + ptb_ad(ji,jj,jk,jn)
371                  ptb_ad(ji,jj,jk,jn) = 0._wp
372                  ptb_ad(ji,jj,jk,jn) = ztdad
373                  pta_ad(ji,jj,jk,jn) = pta_ad(ji,jj,jk,jn) + ztdad
374                  ztnad =  ztnad - 2._wp * ztdad
375                  ptn_ad(ji,jj,jk,jn) = ptn_ad(ji,jj,jk,jn) + ztnad
376                  ztdad = 0._wp
377                  ztnad = 0._wp
378               END DO
379           END DO
380         END DO
381         !
382      END DO
383      !
384   END SUBROUTINE tra_nxt_fix_adj
385
386   SUBROUTINE tra_nxt_vvl_adj( kt )
387      !!----------------------------------------------------------------------
388      !!                   ***  ROUTINE tra_nxt_vvl_adj  ***
389      !!
390      !! ** Purpose :   Time varying volume: apply the Asselin time filter
391      !!                and swap the tracer fields.
392      !!
393      !! ** Method  : - Apply a thickness weighted Asselin time filter on now fields.
394      !!              - save in (ta,sa) a thickness weighted average over the three
395      !!             time levels which will be used to compute rdn and thus the semi-
396      !!             implicit hydrostatic pressure gradient (ln_dynhpg_imp = T)
397      !!              - swap tracer fields to prepare the next time_step.
398      !!                This can be summurized for tempearture as:
399      !!             ztm = (e3t_a*ta+2*e3t_n*tn+e3t_b*tb)   ln_dynhpg_imp = T
400      !!                  /(e3t_a   +2*e3t_n   +e3t_b   )
401      !!             ztm = 0                                otherwise
402      !!             tb  = ( e3t_n*tn + atfp*[ e3t_b*tb - 2 e3t_n*tn + e3t_a*ta ] )
403      !!                  /( e3t_n    + atfp*[ e3t_b    - 2 e3t_n    + e3t_a    ] )
404      !!             tn  = ta
405      !!             ta  = zt        (NB: reset to 0 after eos_bn2 call)
406      !!
407      !! ** Action  : - (tb,sb) and (tn,sn) ready for the next time step
408      !!              - (ta,sa) time averaged (t,s)   (ln_dynhpg_imp = T)
409      !!----------------------------------------------------------------------
410      INTEGER, INTENT(in) ::   kt    ! ocean time-step index
411      !!
412      INTEGER  ::   ji, jj, jk             ! dummy loop indices
413      REAL(wp) ::   ztm , ztc_f , ztf , ztca, ztcn, ztcb   ! temporary scalar
414      REAL(wp) ::   zsm , zsc_f , zsf , zsca, zscn, zscb   !    -         -
415      REAL(wp) ::   ze3mr, ze3fr                           !    -         -
416      REAL(wp) ::   ze3t_b, ze3t_n, ze3t_a, ze3t_f         !    -         -
417      !!----------------------------------------------------------------------
418
419      IF( kt == nitend ) THEN
420         IF(lwp) WRITE(numout,*)
421         IF(lwp) WRITE(numout,*) 'tra_nxt_vvl_adj : time stepping'
422         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
423      ENDIF
424
425      IF(lwp) WRITE(numout,*) "key_vvl net available in tangent yet"
426      CALL abort
427      !
428   END SUBROUTINE tra_nxt_vvl_adj
429
430   SUBROUTINE tra_nxt_adj_tst( kumadt )
431      !!-----------------------------------------------------------------------
432      !!
433      !!          ***  ROUTINE tra_nxt_adj_tst : TEST OF tra_nxt_adj  ***
434      !!
435      !! ** Purpose : Test the adjoint routine.
436      !!
437      !! ** Method  : Verify the scalar product
438      !!
439      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
440      !!
441      !!              where  L   = tangent routine
442      !!                     L^T = adjoint routine
443      !!                     W   = diagonal matrix of scale factors
444      !!                     dx  = input perturbation (random field)
445      !!                     dy  = L dx
446      !!
447       !! History :
448      !!        ! 08-08 (A. Vidard)
449      !!-----------------------------------------------------------------------
450      !! * Modules used
451
452      !! * Arguments
453      INTEGER, INTENT(IN) :: &
454         & kumadt            ! Output unit
455
456      INTEGER ::  &
457         & ji,    &          ! dummy loop indices
458         & jj,    &
459         & jk,    &
460         & jn
461
462      LOGICAL ::  &          ! local variable for time scheme
463         & ll_dynhpg_imp
464
465      INTEGER, DIMENSION(jpi,jpj) :: &
466         & iseed_2d        ! 2D seed for the random number generator
467
468      !! * Local declarations
469      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
470         & zsb_tlin,     &! Tangent input : before salinity
471         & ztb_tlin,     &! Tangent input : before temperature
472         & zsa_tlin,     &! Tangent input : after salinity
473         & zta_tlin,     &! Tangent input : after temperature
474         & zsn_tlin,     &! Tangent input : now salinity
475         & ztn_tlin,     &! Tangent input : now temperature
476         & zsb_tlout,    &! Tangent output: before salinity
477         & ztb_tlout,    &! Tangent output: before temperature
478         & zsa_tlout,    &! Tangent output: after salinity
479         & zta_tlout,    &! Tangent output: after temperature
480         & zsn_tlout,    &! Tangent output: now salinity
481         & ztn_tlout,    &! Tangent output: now temperature
482         & zsb_adin,     &! Adjoint input : before salinity
483         & ztb_adin,     &! Adjoint input : before temperature
484         & zsa_adin,     &! Adjoint input : after salinity
485         & zta_adin,     &! Adjoint input : after temperature
486         & zsn_adin,     &! Adjoint input : now salinity
487         & ztn_adin,     &! Adjoint input : now temperature
488         & zsb_adout,    &! Adjoint output: before salinity
489         & ztb_adout,    &! Adjoint output: before temperature
490         & zsa_adout,    &! Adjoint output: after salinity
491         & zta_adout,    &! Adjoint output: after temperature
492         & zsn_adout,    &! Adjoint output: now salinity
493         & ztn_adout,    &! Adjoint output: now temperature
494         & zr             ! 3D field
495
496      REAL(KIND=wp) ::       &
497         & zsp1,             & ! scalar product involving the tangent routine
498         & zsp1_1,           & ! scalar product involving the tangent routine
499         & zsp1_2,           & ! scalar product involving the tangent routine
500         & zsp1_3,           & ! scalar product involving the tangent routine
501         & zsp1_4,           & ! scalar product involving the tangent routine
502         & zsp1_5,           & ! scalar product involving the tangent routine
503         & zsp1_6,           & ! scalar product involving the tangent routine
504         & zsp2,             & ! scalar product involving the adjoint routine
505         & zsp2_1,           & ! scalar product involving the adjoint routine
506         & zsp2_2,           & ! scalar product involving the adjoint routine
507         & zsp2_3,           & ! scalar product involving the adjoint routine
508         & zsp2_4,           & ! scalar product involving the adjoint routine
509         & zsp2_5,           & ! scalar product involving the adjoint routine
510         & zsp2_6              ! scalar product involving the adjoint routine
511      CHARACTER(LEN=14) ::   &
512         & cl_name
513
514      ALLOCATE( &
515         & zsb_tlin(jpi,jpj,jpk),    &! Tangent input : before salinity
516         & ztb_tlin(jpi,jpj,jpk),    &! Tangent input : before temperature
517         & zsa_tlin(jpi,jpj,jpk),    &! Tangent input : after salinity
518         & zta_tlin(jpi,jpj,jpk),    &! Tangent input : after temperature
519         & zsn_tlin(jpi,jpj,jpk),    &! Tangent input : now salinity
520         & ztn_tlin(jpi,jpj,jpk),    &! Tangent input : now temperature
521         & zsb_tlout(jpi,jpj,jpk),   &! Tangent output: before salinity
522         & ztb_tlout(jpi,jpj,jpk),   &! Tangent output: before temperature
523         & zsa_tlout(jpi,jpj,jpk),   &! Tangent output: after salinity
524         & zta_tlout(jpi,jpj,jpk),   &! Tangent output: after temperature
525         & zsn_tlout(jpi,jpj,jpk),   &! Tangent output: now salinity
526         & ztn_tlout(jpi,jpj,jpk),   &! Tangent output: now temperature
527         & zsb_adin(jpi,jpj,jpk),    &! Adjoint input : before salinity
528         & ztb_adin(jpi,jpj,jpk),    &! Adjoint input : before temperature
529         & zsa_adin(jpi,jpj,jpk),    &! Adjoint input : after salinity
530         & zta_adin(jpi,jpj,jpk),    &! Adjoint input : after temperature
531         & zsn_adin(jpi,jpj,jpk),    &! Adjoint input : now salinity
532         & ztn_adin(jpi,jpj,jpk),    &! Adjoint input : now temperature
533         & zsb_adout(jpi,jpj,jpk),   &! Adjoint output: before salinity
534         & ztb_adout(jpi,jpj,jpk),   &! Adjoint output: before temperature
535         & zsa_adout(jpi,jpj,jpk),   &! Adjoint output: after salinity
536         & zta_adout(jpi,jpj,jpk),   &! Adjoint output: after temperature
537         & zsn_adout(jpi,jpj,jpk),   &! Adjoint output: now salinity
538         & ztn_adout(jpi,jpj,jpk),   &! Adjoint output: now temperature
539         & zr       (jpi,jpj,jpk)    &! 3D field
540         & )
541
542      ll_dynhpg_imp = ln_dynhpg_imp     ! store namelist define time scheme
543
544      DO jn = 1, 2
545
546         IF ( jn .EQ. 1) ln_dynhpg_imp = .TRUE.
547         IF ( jn .EQ. 2) ln_dynhpg_imp = .FALSE.
548
549         !==================================================================
550         ! 1) dx = ( tb_tl, tn_tl, ta_tl,       dy = ( tb_tl, tn_tl, ta_tl,
551         !           sb_tl, sn_tl, sa_tl  ) and        sb_tl, sn_tl, sa_tl )
552         !==================================================================
553
554         !--------------------------------------------------------------------
555         ! Reset the tangent and adjoint variables
556         !--------------------------------------------------------------------
557         tsb_tl(:,:,:,:) = 0.0_wp
558         tsa_tl(:,:,:,:) = 0.0_wp
559         tsn_tl(:,:,:,:) = 0.0_wp
560         tsb_ad(:,:,:,:) = 0.0_wp
561         tsa_ad(:,:,:,:) = 0.0_wp
562         tsn_ad(:,:,:,:) = 0.0_wp
563         zsb_tlin(:,:,:) = 0.0_wp
564         ztb_tlin(:,:,:) = 0.0_wp
565         zsa_tlin(:,:,:) = 0.0_wp
566         zta_tlin(:,:,:) = 0.0_wp
567         zsn_tlin(:,:,:) = 0.0_wp
568         ztn_tlin(:,:,:) = 0.0_wp
569
570         r2dtra(:) =  2.* rdttra(:) ! initialization
571
572         CALL grid_random(  zr, 'T', 0.0_wp, stds )
573         DO jk = 1, jpk
574            DO jj = nldj, nlej
575               DO ji = nldi, nlei
576                  zsb_tlin(ji,jj,jk) = zr(ji,jj,jk)
577               END DO
578            END DO
579         END DO
580
581         CALL grid_random(  zr, 'T', 0.0_wp, stdt )
582         DO jk = 1, jpk
583            DO jj = nldj, nlej
584               DO ji = nldi, nlei
585                  ztb_tlin(ji,jj,jk) = zr(ji,jj,jk)
586               END DO
587            END DO
588         END DO
589
590         CALL grid_random(  zr, 'T', 0.0_wp, stds )
591         DO jk = 1, jpk
592            DO jj = nldj, nlej
593               DO ji = nldi, nlei
594                  zsa_tlin(ji,jj,jk) = zr(ji,jj,jk)
595               END DO
596            END DO
597         END DO
598
599         CALL grid_random(  zr, 'T', 0.0_wp, stdt )
600         DO jk = 1, jpk
601            DO jj = nldj, nlej
602               DO ji = nldi, nlei
603                  zta_tlin(ji,jj,jk) = zr(ji,jj,jk)
604               END DO
605            END DO
606         END DO
607
608         CALL grid_random(  zr, 'T', 0.0_wp, stds )
609         DO jk = 1, jpk
610            DO jj = nldj, nlej
611               DO ji = nldi, nlei
612                  zsn_tlin(ji,jj,jk) = zr(ji,jj,jk)
613               END DO
614            END DO
615         END DO
616
617         CALL grid_random(  zr, 'T', 0.0_wp, stdt )
618         DO jk = 1, jpk
619            DO jj = nldj, nlej
620               DO ji = nldi, nlei
621                  ztn_tlin(ji,jj,jk) = zr(ji,jj,jk)
622               END DO
623            END DO
624         END DO
625
626         tsb_tl(:,:,:,jp_sal) = zsb_tlin(:,:,:)
627         tsb_tl(:,:,:,jp_tem) = ztb_tlin(:,:,:)
628         tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:)
629         tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
630         tsn_tl(:,:,:,jp_sal) = zsn_tlin(:,:,:)
631         tsn_tl(:,:,:,jp_tem) = ztn_tlin(:,:,:)
632
633         CALL tra_nxt_tan( nit000 + 1 )
634
635         zsa_tlout(:,:,:) = tsa_tl(:,:,:,jp_sal)
636         zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
637         zsb_tlout(:,:,:) = tsb_tl(:,:,:,jp_sal)
638         ztb_tlout(:,:,:) = tsb_tl(:,:,:,jp_tem)
639         zsn_tlout(:,:,:) = tsn_tl(:,:,:,jp_sal)
640         ztn_tlout(:,:,:) = tsn_tl(:,:,:,jp_tem)
641
642         !--------------------------------------------------------------------
643         ! Initialize the adjoint variables: dy^* = W dy
644         !--------------------------------------------------------------------
645
646         DO jk = 1, jpk
647            DO jj = nldj, nlej
648               DO ji = nldi, nlei
649                  zsa_adin(ji,jj,jk)     = zsa_tlout(ji,jj,jk) &
650                       &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
651                       &                   * tmask(ji,jj,jk) * wesp_s(jk)
652                  zta_adin(ji,jj,jk)     = zta_tlout(ji,jj,jk) &
653                       &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
654                       &                   * tmask(ji,jj,jk) * wesp_t(jk)
655                  zsb_adin(ji,jj,jk)     = zsb_tlout(ji,jj,jk) &
656                       &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
657                       &                   * tmask(ji,jj,jk) * wesp_s(jk)
658                  ztb_adin(ji,jj,jk)     = ztb_tlout(ji,jj,jk) &
659                       &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
660                       &                   * tmask(ji,jj,jk) * wesp_t(jk)
661                  zsn_adin(ji,jj,jk)     = zsn_tlout(ji,jj,jk) &
662                       &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
663                       &                   * tmask(ji,jj,jk) * wesp_s(jk)
664                  ztn_adin(ji,jj,jk)     = ztn_tlout(ji,jj,jk) &
665                       &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
666                       &                   * tmask(ji,jj,jk) * wesp_t(jk)
667               END DO
668            END DO
669         END DO
670
671         !--------------------------------------------------------------------
672         ! Compute the scalar product: ( L dx )^T W dy
673         !--------------------------------------------------------------------
674
675         zsp1_1 = DOT_PRODUCT( zsa_tlout    , zsa_adin     )
676         zsp1_2 = DOT_PRODUCT( zta_tlout    , zta_adin     )
677         zsp1_3 = DOT_PRODUCT( zsb_tlout    , zsb_adin     )
678         zsp1_4 = DOT_PRODUCT( ztb_tlout    , ztb_adin     )
679         zsp1_5 = DOT_PRODUCT( zsn_tlout    , zsn_adin     )
680         zsp1_6 = DOT_PRODUCT( ztn_tlout    , ztn_adin     )
681         zsp1   = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4 + zsp1_5 + zsp1_6
682
683         !--------------------------------------------------------------------
684         ! Call the adjoint routine: dx^* = L^T dy^*
685         !--------------------------------------------------------------------
686
687         tsa_ad(:,:,:,jp_sal)     = zsa_adin(:,:,:)
688         tsa_ad(:,:,:,jp_tem)     = zta_adin(:,:,:)
689         tsb_ad(:,:,:,jp_sal)     = zsb_adin(:,:,:)
690         tsb_ad(:,:,:,jp_tem)     = ztb_adin(:,:,:)
691         tsn_ad(:,:,:,jp_sal)     = zsn_adin(:,:,:)
692         tsn_ad(:,:,:,jp_tem)     = ztn_adin(:,:,:)
693
694         CALL tra_nxt_adj ( nit000 + 1 )
695
696         zsb_adout(:,:,:) = tsb_ad(:,:,:,jp_sal)
697         ztb_adout(:,:,:) = tsb_ad(:,:,:,jp_tem)
698         zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal)
699         zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
700         zsn_adout(:,:,:) = tsn_ad(:,:,:,jp_sal)
701         ztn_adout(:,:,:) = tsn_ad(:,:,:,jp_tem)
702
703         !--------------------------------------------------------------------
704         ! Compute the scalar product: dx^T L^T W dy
705         !--------------------------------------------------------------------
706
707         zsp2_1 = DOT_PRODUCT( zsb_tlin  , zsb_adout   )
708         zsp2_2 = DOT_PRODUCT( ztb_tlin  , ztb_adout   )
709         zsp2_3 = DOT_PRODUCT( zsa_tlin  , zsa_adout   )
710         zsp2_4 = DOT_PRODUCT( zta_tlin  , zta_adout   )
711         zsp2_5 = DOT_PRODUCT( zsn_tlin  , zsn_adout   )
712         zsp2_6 = DOT_PRODUCT( ztn_tlin  , ztn_adout   )
713
714         zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4 + zsp2_5 + zsp2_6
715
716         ! Compare the scalar products
717
718         ! 14 char:'12345678901234'
719         IF ( jn .EQ. 1) cl_name = 'tra_nxt_adj T1'
720         IF ( jn .EQ. 2) cl_name = 'tra_nxt_adj T2'
721         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
722
723      ENDDO
724
725      ln_dynhpg_imp = ll_dynhpg_imp   ! restore initial value of ln_dynhpg_imp
726
727      DEALLOCATE( &
728         & zsb_tlin,     &
729         & ztb_tlin,     &
730         & zsa_tlin,     &
731         & zta_tlin,     &
732         & zsn_tlin,     &
733         & ztn_tlin,     &
734         & zsb_tlout,    &
735         & ztb_tlout,    &
736         & zsa_tlout,    &
737         & zta_tlout,    &
738         & zsn_tlout,    &
739         & ztn_tlout,    &
740         & zsb_adin,     &
741         & ztb_adin,     &
742         & zsa_adin,     &
743         & zta_adin,     &
744         & zsn_adin,     &
745         & ztn_adin,     &
746         & zsb_adout,    &
747         & ztb_adout,    &
748         & zsa_adout,    &
749         & zta_adout,    &
750         & zsn_adout,    &
751         & ztn_adout,    &
752         & zr            &
753         & )
754
755     END SUBROUTINE tra_nxt_adj_tst
756   !!======================================================================
757#endif
758END MODULE tranxt_tam
Note: See TracBrowser for help on using the repository browser.