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.
tradmp_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/tradmp_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: 30.0 KB
Line 
1MODULE tradmp_tam
2#ifdef key_tam
3   !!======================================================================
4   !!                       ***  MODULE  tradmp_tam  ***
5   !! Ocean physics: internal restoring trend on active tracers (T and S)
6   !!                Tangent and Adjoint Module
7   !!======================================================================
8   !! History of the direct module:
9   !!            OPA  ! 1991-03  (O. Marti, G. Madec)  Original code
10   !!                 ! 1992-06  (M. Imbard)  doctor norme
11   !!                 ! 1996-01  (G. Madec)  statement function for e3
12   !!                 ! 1997-05  (G. Madec)  macro-tasked on jk-slab
13   !!                 ! 1998-07  (M. Imbard, G. Madec) ORCA version
14   !!            7.0  ! 2001-02  (M. Imbard)  cofdis, Original code
15   !!            8.1  ! 2001-02  (G. Madec, E. Durand)  cleaning
16   !!  NEMO      1.0  ! 2002-08  (G. Madec, E. Durand)  free form + modules
17   !!                 ! 2003-08  (M. Balmaseda)  mods to allow eq. damping
18   !!            3.2  ! 2009-08  (G. Madec, C. Talandier)  DOCTOR norm for namelist parameter
19   !! History of the TAM:
20   !!                 ! 2008-09  (A. Vidard) tangent and adjoint module
21   !!                                       of the 03-08 version
22   !!        NEMO 3.2 ! 2010-04 (F. Vigilant) 3.2 version
23   !!        NEMO 3.4 ! 2012-07 (P.-A. Bouttier) 3.4 Version@
24   !!----------------------------------------------------------------------
25   !!----------------------------------------------------------------------
26   !!   tra_dmp      : update the tracer trend with the internal damping
27   !!   tra_dmp_init : initialization, namlist read, parameters control
28   !!   dtacof_zoom  : restoring coefficient for zoom domain
29   !!   dtacof       : restoring coefficient for global domain
30   !!   cofdis       : compute the distance to the coastline
31   !!----------------------------------------------------------------------
32   USE par_oce
33   USE oce_tam
34   USE dom_oce
35   USE zdf_oce
36   USE in_out_manager
37   USE phycst
38   USE lib_mpp
39   USE tradmp
40   USE prtctl
41   USE gridrandom
42   USE dotprodfld
43   USE dtatsd
44   USE zdfmxl
45   USE tstool_tam
46   USE paresp
47   USE lib_mpp
48   USE wrk_nemo
49   USE timing
50
51   IMPLICIT NONE
52   PRIVATE
53
54   PUBLIC tra_dmp_tan      ! routine called by step_tam.F90
55   PUBLIC tra_dmp_adj      ! routine called by step_tam.F90
56   PUBLIC tra_dmp_init_tam
57   PUBLIC tra_dmp_adj_tst  ! routine called by tst.F90
58   !LOGICAL, PUBLIC            ::   lk_tradmp = .TRUE.     !: internal damping flag
59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: &
60     & strdmp_tl,   & !: damping salinity trend (psu/s)
61     & ttrdmp_tl,   & !: damping temperature trend (psu/s)
62     & ttrdmp_ad,   & !: damping temperature trend (psu/s)
63     & strdmp_ad      !: damping salinity trend (psu/s)
64
65
66   LOGICAL :: lfirst = .TRUE.        !: flag for initialisation
67   !                                 !!* Namelist namtra_dmp : T & S newtonian damping *
68   INTEGER  ::   nn_hdmp =   -1      ! = 0/-1/'latitude' for damping over T and S
69   INTEGER  ::   nn_zdmp =    0      ! = 0/1/2 flag for damping in the mixed layer
70   REAL(wp) ::   rn_surf =   50._wp  ! surface time scale for internal damping        [days]
71   REAL(wp) ::   rn_bot  =  360._wp  ! bottom time scale for internal damping         [days]
72   REAL(wp) ::   rn_dep  =  800._wp  ! depth of transition between rn_surf and rn_bot [meters]
73   INTEGER  ::   nn_file =    2      ! = 1 create a damping.coeff NetCDF file
74   LOGICAL, PRIVATE, SAVE :: ll_alloctl = .FALSE.
75   LOGICAL, PRIVATE, SAVE :: ll_allocad = .FALSE.
76   !! * Substitutions
77#  include "domzgr_substitute.h90"
78#  include "vectopt_loop_substitute.h90"
79
80
81CONTAINS
82
83   INTEGER FUNCTION tra_dmp_alloc_tam( kmode )
84      !!----------------------------------------------------------------------
85      !!                ***  FUNCTION tra_dmp_alloc_tam  ***
86      !!----------------------------------------------------------------------
87      INTEGER, OPTIONAL :: kmode
88      INTEGER :: ierr(2)
89      integer :: jmode
90
91      IF ( PRESENT(kmode) ) THEN
92         jmode = kmode 
93      ELSE
94         jmode = 0
95      END IF
96
97      IF ( ( jmode == 0 ) .OR. ( jmode == 1 ) .AND. ( .NOT. ll_alloctl) ) THEN
98         ALLOCATE( strdmp_tl(jpi,jpj,jpk) , ttrdmp_tl(jpi,jpj,jpk), STAT= ierr(1) )
99      END IF
100      IF ( ( jmode == 0 ) .OR. ( jmode == 2 ) .AND. ( .NOT. ll_allocad) ) THEN
101         ALLOCATE( strdmp_ad(jpi,jpj,jpk) , ttrdmp_ad(jpi,jpj,jpk), STAT= ierr(2) )
102      END IF
103      tra_dmp_alloc_tam = SUM( ierr )
104      !
105      IF( lk_mpp            )   CALL mpp_sum ( tra_dmp_alloc_tam )
106      IF( tra_dmp_alloc_tam > 0 )   CALL ctl_warn('tra_dmp_alloc_tam: allocation of arrays failed')
107      !
108   END FUNCTION tra_dmp_alloc_tam
109
110   INTEGER FUNCTION tra_dmp_dealloc_tam( kmode )
111      !!----------------------------------------------------------------------
112      !!                ***  FUNCTION tra_dmp_alloc_tam  ***
113      !!----------------------------------------------------------------------
114      INTEGER, OPTIONAL :: kmode
115      INTEGER :: ierr(2)
116      integer :: jmode
117
118      IF ( PRESENT( kmode) ) THEN
119         jmode = kmode 
120      ELSE
121         jmode = 0
122      END IF
123
124      IF ( ( jmode == 0 ) .OR. ( jmode == 1 ) .AND. ( ll_alloctl) ) THEN
125         DEALLOCATE( strdmp_tl, ttrdmp_tl, STAT= ierr(1) )
126      END IF
127      IF ( ( jmode == 0 ) .OR. ( jmode == 2 ) .AND. ( ll_allocad) ) THEN
128         DEALLOCATE( strdmp_ad, ttrdmp_ad, STAT= ierr(2) )
129      END IF
130      tra_dmp_dealloc_tam = SUM( ierr )
131      !
132      IF( lk_mpp            )   CALL mpp_sum ( tra_dmp_dealloc_tam )
133      IF( tra_dmp_dealloc_tam > 0 )   CALL ctl_warn('tra_dmp_dealloc_tam: deallocation of arrays failed')
134      !
135   END FUNCTION tra_dmp_dealloc_tam
136
137   SUBROUTINE tra_dmp_tan( kt )
138      !!----------------------------------------------------------------------
139      !!                   ***  ROUTINE tra_dmp_tan  ***
140      !!
141      !! ** Purpose of direct routine:
142      !!      Compute the tracer trend due to a newtonian damping
143      !!      of the tracer field towards given data field and add it to the
144      !!      general tracer trends.
145      !!
146      !! ** Method  :   Newtonian damping towards t_dta and s_dta computed
147      !!      and add to the general tracer trends:
148      !!                     ta = ta + resto * (t_dta - tb)
149      !!                     sa = sa + resto * (s_dta - sb)
150      !!         The trend is computed either throughout the water column
151      !!      (nlmdmp=0) or in area of weak vertical mixing (nlmdmp=1) or
152      !!      below the well mixed layer (nlmdmp=2)
153      !!
154      !! ** Action  : - (ta,sa)   tracer trends updated with the damping trend
155      !!
156      !! ASSUME key_zdfcst_tam
157      !!----------------------------------------------------------------------
158      !!
159      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
160      !!
161      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
162      REAL(wp) ::   ztest, ztatl, zsatl       ! temporary scalars
163      !!----------------------------------------------------------------------
164      !
165      IF( nn_timing == 1 )  CALL timing_start( 'tra_dmp_tan')
166      !
167      ! 1. Newtonian damping trends on tracer fields
168      ! --------------------------------------------
169      !    compute the newtonian damping trends depending on nmldmp
170
171      SELECT CASE ( nn_zdmp )
172      !
173      CASE( 0 )                ! newtonian damping throughout the water column
174         DO jk = 1, jpkm1
175            DO jj = 2, jpjm1
176               DO ji = fs_2, fs_jpim1   ! vector opt.
177                  ztatl = - resto(ji,jj,jk) *  tsb_tl(ji,jj,jk,jp_tem)
178                  zsatl = - resto(ji,jj,jk) *  tsb_tl(ji,jj,jk,jp_sal)
179                  ! add the trends to the general tracer trends
180                  tsa_tl(ji,jj,jk,jp_tem) = tsa_tl(ji,jj,jk,jp_tem) + ztatl
181                  tsa_tl(ji,jj,jk,jp_sal) = tsa_tl(ji,jj,jk,jp_sal) + zsatl
182                  ! save the salinity trend (used in flx to close the salt budget)
183                  strdmp_tl(ji,jj,jk) = zsatl
184                  ttrdmp_tl(ji,jj,jk) = ztatl
185               END DO
186            END DO
187         END DO
188         !
189      CASE ( 1 )                ! no damping in the turbocline (avt > 5 cm2/s)
190         DO jk = 1, jpkm1
191            DO jj = 2, jpjm1
192               DO ji = fs_2, fs_jpim1   ! vector opt.
193                  ztest = avt(ji,jj,jk) - 5.e-4_wp
194      !! ASSUME key_zdfcst_tam
195                  IF( ztest < 0. ) THEN
196                     ztatl = - resto(ji,jj,jk) * tsb_tl(ji,jj,jk,jp_tem)
197                     zsatl = - resto(ji,jj,jk) * tsb_tl(ji,jj,jk,jp_sal)
198                  ELSE
199                     ztatl = 0.0_wp
200                     zsatl = 0.0_wp
201                  ENDIF
202                  ! add the trends to the general tracer trends
203                  tsa_tl(ji,jj,jk,jp_tem) = tsa_tl(ji,jj,jk,jp_tem) + ztatl
204                  tsa_tl(ji,jj,jk,jp_sal) = tsa_tl(ji,jj,jk,jp_sal) + zsatl
205                  ! save the salinity trend (used in flx to close the salt budget)
206                  strdmp_tl(ji,jj,jk) = zsatl
207                  ttrdmp_tl(ji,jj,jk) = ztatl
208               END DO
209            END DO
210         END DO
211         !
212      CASE ( 2 )                ! no damping in the mixed layer
213         DO jk = 1, jpkm1
214            DO jj = 2, jpjm1
215               DO ji = fs_2, fs_jpim1   ! vector opt.
216                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
217                     ztatl = - resto(ji,jj,jk) * tsb_tl(ji,jj,jk,jp_tem)
218                     zsatl = - resto(ji,jj,jk) * tsb_tl(ji,jj,jk,jp_sal)
219                  ELSE
220                     ztatl = 0.e0
221                     zsatl = 0.e0
222                  ENDIF
223                  ! add the trends to the general tracer trends
224                  tsa_tl(ji,jj,jk,jp_tem) = tsa_tl(ji,jj,jk,jp_tem) + ztatl
225                  tsa_tl(ji,jj,jk,jp_sal) = tsa_tl(ji,jj,jk,jp_sal) + zsatl
226                  ! save the salinity trend (used in flx to close the salt budget)
227                  strdmp_tl(ji,jj,jk) = zsatl
228                  ttrdmp_tl(ji,jj,jk) = ztatl
229               END DO
230            END DO
231         END DO
232         !
233      END SELECT
234      !
235      IF( nn_timing == 1 )  CALL timing_stop( 'tra_dmp_tan')
236      !
237   END SUBROUTINE tra_dmp_tan
238   SUBROUTINE tra_dmp_adj( kt )
239      !!----------------------------------------------------------------------
240      !!                   ***  ROUTINE tra_dmp_adj  ***
241      !!
242      !! ** Purpose of direct routine:
243      !!      Compute the tracer trend due to a newtonian damping
244      !!      of the tracer field towards given data field and add it to the
245      !!      general tracer trends.
246      !!
247      !! ** Method  :   Newtonian damping towards t_dta and s_dta computed
248      !!      and add to the general tracer trends:
249      !!                     ta = ta + resto * (t_dta - tb)
250      !!                     sa = sa + resto * (s_dta - sb)
251      !!         The trend is computed either throughout the water column
252      !!      (nlmdmp=0) or in area of weak vertical mixing (nlmdmp=1) or
253      !!      below the well mixed layer (nlmdmp=2)
254      !!
255      !! ** Action  : - (ta,sa)   tracer trends updated with the damping trend
256      !!              - save the trends in (ttrd,strd) ('key_trdtra')
257      !! ASSUME key_zdfcst_tam
258      !!----------------------------------------------------------------------
259      !!
260      INTEGER, INTENT( in ) ::   kt   ! ocean time-step index
261      !!
262      INTEGER  ::   ji, jj, jk, jn            ! dummy loop indices
263      REAL(wp) ::   ztest, ztaad, zsaad       ! temporary scalars
264      !!----------------------------------------------------------------------
265      !
266      IF( nn_timing == 1 )  CALL timing_start( 'tra_dmp_adj')
267      !
268      ! 1. Newtonian damping trends on tracer fields
269      ! --------------------------------------------
270      !    compute the newtonian damping trends depending on nmldmp
271      ztaad = 0.0_wp
272      zsaad = 0.0_wp
273
274      SELECT CASE ( nn_zdmp )
275      !
276      CASE( 0 )                ! newtonian damping throughout the water column
277         DO jk = 1, jpkm1
278            DO jj = 2, jpjm1
279               DO ji = fs_2, fs_jpim1   ! vector opt.
280                  ! save the salinity trend (used in flx to close the salt budget)
281                  ztaad = ztaad + ttrdmp_ad(ji,jj,jk)
282                  ttrdmp_ad(ji,jj,jk) = 0.0_wp
283                  zsaad = zsaad + strdmp_ad(ji,jj,jk)
284                  strdmp_ad(ji,jj,jk) = 0.0_wp
285                  ! add the trends to the general tracer trends
286                  ztaad = tsa_ad(ji,jj,jk,jp_tem) + ztaad
287                  zsaad = tsa_ad(ji,jj,jk,jp_sal) + zsaad
288
289                  tsb_ad(ji,jj,jk,jp_tem)  = tsb_ad(ji,jj,jk,jp_tem) - ztaad * resto(ji,jj,jk)
290                  tsb_ad(ji,jj,jk,jp_sal)  = tsb_ad(ji,jj,jk,jp_sal) - zsaad * resto(ji,jj,jk)
291                  ztaad = 0.0_wp
292                  zsaad = 0.0_wp
293                END DO
294            END DO
295         END DO
296         !
297      CASE ( 1 )                ! no damping in the turbocline (avt > 5 cm2/s)
298         DO jk = 1, jpkm1
299            DO jj = 2, jpjm1
300               DO ji = fs_2, fs_jpim1   ! vector opt.
301                  ! save the salinity trend (used in flx to close the salt budget)
302                  ztaad = ztaad + ttrdmp_ad(ji,jj,jk)
303                  ttrdmp_ad(ji,jj,jk) = 0.0_wp
304                  zsaad = zsaad + strdmp_ad(ji,jj,jk)
305                  strdmp_ad(ji,jj,jk) = 0.0_wp
306                  ! add the trends to the general tracer trends
307                  ztaad = tsa_ad(ji,jj,jk,jp_tem) + ztaad
308                  zsaad = tsa_ad(ji,jj,jk,jp_sal) + zsaad
309                  ztest = avt(ji,jj,jk) - 5.e-4
310      !! ASSUME key_zdfcst_tam
311                  IF( ztest < 0. ) THEN
312                     tsb_ad(ji,jj,jk,jp_tem)  = tsb_ad(ji,jj,jk,jp_tem) - ztaad * resto(ji,jj,jk)
313                     tsb_ad(ji,jj,jk,jp_sal)  = tsb_ad(ji,jj,jk,jp_sal) - zsaad * resto(ji,jj,jk)
314                     ztaad = 0.0_wp
315                     zsaad = 0.0_wp
316                  ELSE
317                     ztaad = 0.0_wp
318                     zsaad = 0.0_wp
319                  ENDIF
320               END DO
321            END DO
322         END DO
323         !
324      CASE ( 2 )                ! no damping in the mixed layer
325         DO jk = 1, jpkm1
326            DO jj = 2, jpjm1
327               DO ji = fs_2, fs_jpim1   ! vector opt.
328                  ! save the salinity trend (used in flx to close the salt budget)
329                  ztaad = ztaad + ttrdmp_ad(ji,jj,jk)
330                  ttrdmp_ad(ji,jj,jk) = 0.0_wp
331                  zsaad = zsaad + strdmp_ad(ji,jj,jk)
332                  strdmp_ad(ji,jj,jk) = 0.0_wp
333                  ! add the trends to the general tracer trends
334                  ztaad = tsa_ad(ji,jj,jk,jp_tem) + ztaad
335                  zsaad = tsa_ad(ji,jj,jk,jp_sal) + zsaad
336                  IF( fsdept(ji,jj,jk) >= hmlp (ji,jj) ) THEN
337                     tsb_ad(ji,jj,jk,jp_tem)  = tsb_ad(ji,jj,jk,jp_tem) - ztaad * resto(ji,jj,jk)
338                     tsb_ad(ji,jj,jk,jp_sal)  = tsb_ad(ji,jj,jk,jp_sal) - zsaad * resto(ji,jj,jk)
339                     ztaad = 0.0_wp
340                     zsaad = 0.0_wp
341                  ELSE
342                     ztaad = 0.0_wp
343                     zsaad = 0.0_wp
344                  ENDIF
345               END DO
346            END DO
347         END DO
348         !
349      END SELECT
350      !
351      IF( nn_timing == 1 )  CALL timing_stop( 'tra_dmp_adj')
352      !
353   END SUBROUTINE tra_dmp_adj
354   SUBROUTINE tra_dmp_init_tam
355      !!----------------------------------------------------------------------
356      !!                  ***  ROUTINE tra_dmp_init_tam  ***
357      !!
358      !! ** Purpose :   Initialization for the newtonian damping
359      !!
360      !! ** Method  :   read the nammbf namelist and check the parameters
361      !!      called by tra_dmp at the first timestep (nit000)
362      !!----------------------------------------------------------------------
363
364      !NAMELIST/namtra_dmp/ nn_hdmp, nn_file, nn_zdmp, rn_surf, rn_bot, rn_dep
365      !!----------------------------------------------------------------------
366
367      IF (lfirst) THEN
368
369         !REWIND ( numnam )                  ! Read Namelist namtdp : temperature and salinity damping term
370         !READ   ( numnam, namtra_dmp )
371         !IF( lzoom )   nn_zdmp = 0           ! restoring to climatology at closed north or south boundaries
372
373         !IF(lwp) THEN                       ! Namelist print
374            !WRITE(numout,*)
375            !WRITE(numout,*) 'tra_dmp_init_tam : T and S newtonian damping'
376            !WRITE(numout,*) '~~~~~~~~~~~'
377            !WRITE(numout,*) '      Namelist namtra_dmp : set damping parameter'
378            !WRITE(numout,*) '      T and S damping option         nn_hdmp = ', nn_hdmp
379            !WRITE(numout,*) '      mixed layer damping option     nn_zdmp = ', nn_zdmp, '(zoom: forced to 0)'
380            !WRITE(numout,*) '      surface time scale (days)      rn_surf = ', rn_surf
381            !WRITE(numout,*) '      bottom time scale (days)       rn_bot  = ', rn_bot
382            !WRITE(numout,*) '      depth of transition (meters)   rn_dep  = ', rn_dep
383            !WRITE(numout,*) '      create a damping.coeff file    nn_file = ', nn_file
384         !ENDIF
385
386         !! ** modified to allow damping at the equator
387         IF( ln_tradmp ) THEN               ! initialization for T-S damping
388         !
389            IF( tra_dmp_alloc_tam() /= 0 )   CALL ctl_stop( 'STOP', 'tra_dmp_init_tam: unable to allocate arrays' )
390            SELECT CASE ( nn_hdmp )
391            CASE (  -1  )   ;   IF(lwp) WRITE(numout,*) '   tracer damping in the Med & Red seas only'
392            CASE ( 0:90 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping poleward of', nn_hdmp, ' degrees'
393            CASE DEFAULT
394               IF(lwp) WRITE(numout,*) '          tracer damping disabled', nn_hdmp
395               !         WRITE(ctmp1,*) '          bad flag value for ndmp = ', ndmp
396               !         CALL ctl_stop(ctmp1)
397            END SELECT
398
399            SELECT CASE ( nn_zdmp )
400            CASE ( 0 )   ;   IF(lwp) WRITE(numout,*) '   tracer damping throughout the water column'
401            CASE ( 1 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the turbocline (avt > 5 cm2/s)'
402            CASE ( 2 )   ;   IF(lwp) WRITE(numout,*) '   no tracer damping in the mixed layer'
403            CASE DEFAULT
404               WRITE(ctmp1,*) '   bad flag value for nn_zdmp = ', nn_zdmp
405               CALL ctl_stop(ctmp1)
406            END SELECT
407            !
408            !* We don't need data for TAM *!
409            !IF( .NOT. ln_tsd_tradmp )
410               !CALL ctl_warn( 'tra_dmp_init_tam: read T-S data not initialized, we force ln_tsd_tradmp=T' )
411            !ENDIF
412            !                          ! Damping coefficients initialization
413            !IF( lzoom ) THEN   ;   CALL dtacof_zoom( resto )
414            !ELSE               ;   CALL dtacof( nn_hdmp, rn_surf, rn_bot, rn_dep, nn_file, 'TRA', resto )
415            !ENDIF
416            !
417            lfirst = .FALSE.
418         END IF
419         ttrdmp_tl(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
420         ttrdmp_ad(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
421         strdmp_tl(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
422         strdmp_ad(:,:,:) = 0.0_wp       ! internal damping salinity trend (used in ocesbc)
423      END IF
424   END SUBROUTINE tra_dmp_init_tam
425
426   SUBROUTINE tra_dmp_adj_tst ( kumadt )
427      !!-----------------------------------------------------------------------
428      !!
429      !!          ***  ROUTINE tra_sbc_adj_tst : TEST OF tra_sbc_adj  ***
430      !!
431      !! ** Purpose : Test the adjoint routine.
432      !!
433      !! ** Method  : Verify the scalar product
434      !!
435      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
436      !!
437      !!              where  L   = tangent routine
438      !!                     L^T = adjoint routine
439      !!                     W   = diagonal matrix of scale factors
440      !!                     dx  = input perturbation (random field)
441      !!                     dy  = L dx
442      !!
443       !! History :
444      !!        ! 08-08 (A. Vidard)
445      !!-----------------------------------------------------------------------
446      !! * Modules used
447
448      !! * Arguments
449      INTEGER, INTENT(IN) :: &
450         & kumadt             ! Output unit
451
452      INTEGER ::  &
453         & istp,  &
454         & jstp,  &
455         & ji,    &        ! dummy loop indices
456         & jj,    &
457         & jk
458
459      !! * Local declarations
460      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
461         & zsb_tlin,     &! Tangent input : before salinity
462         & ztb_tlin,     &! Tangent input : before temperature
463         & zsa_tlin,     &! Tangent input : after salinity
464         & zta_tlin,     &! Tangent input : after temperature
465         & zsa_tlout,    &! Tangent output: after salinity
466         & zta_tlout,    &! Tangent output: after temperature
467         & zstrdmp_tlout,&! Tangent output: salinity trend
468         & zttrdmp_tlout,&! Tangent output: temperature trend
469         & zsb_adout,    &! Adjoint output : before salinity
470         & ztb_adout,    &! Adjoint output : before temperature
471         & zsa_adout,    &! Adjoint output : after salinity
472         & zta_adout,    &! Adjoint output : after temperature
473         & zsa_adin,     &! Adjoint input : after salinity
474         & zta_adin,     &! Adjoint input : after temperature
475         & zstrdmp_adin, &! Adjoint input : salinity trend
476         & zttrdmp_adin, &! Tangent output: temperature trend
477         & z3r            ! 3D field
478
479      REAL(KIND=wp) ::       &
480         & zsp1,             & ! scalar product involving the tangent routine
481         & zsp1_1,           & ! scalar product involving the tangent routine
482         & zsp1_2,           & ! scalar product involving the tangent routine
483         & zsp1_3,           & ! scalar product involving the tangent routine
484         & zsp1_4,           & ! scalar product involving the tangent routine
485         & zsp2,             & ! scalar product involving the adjoint routine
486         & zsp2_1,           & ! scalar product involving the adjoint routine
487         & zsp2_2,           & ! scalar product involving the adjoint routine
488         & zsp2_3,           & ! scalar product involving the adjoint routine
489         & zsp2_4              ! scalar product involving the adjoint routine
490      CHARACTER(LEN=14) ::   &
491         & cl_name
492
493      ALLOCATE( &
494         & ztb_tlin(jpi,jpj,jpk),     &
495         & zsb_tlin(jpi,jpj,jpk),     &
496         & zsa_tlin(jpi,jpj,jpk),     &
497         & zta_tlin(jpi,jpj,jpk),     &
498         & zsa_tlout(jpi,jpj,jpk),    &
499         & zta_tlout(jpi,jpj,jpk),    &
500         & zstrdmp_tlout(jpi,jpj,jpk),&
501         & zttrdmp_tlout(jpi,jpj,jpk),&
502         & ztb_adout(jpi,jpj,jpk),    &
503         & zsb_adout(jpi,jpj,jpk),    &
504         & zsa_adout(jpi,jpj,jpk),    &
505         & zta_adout(jpi,jpj,jpk),    &
506         & zsa_adin(jpi,jpj,jpk),     &
507         & zta_adin(jpi,jpj,jpk),     &
508         & zttrdmp_adin(jpi,jpj,jpk), &
509         & zstrdmp_adin(jpi,jpj,jpk), &
510         & z3r         (jpi,jpj,jpk)  &
511         & )
512      ! Test for time steps nit000 and nit000 + 1 (the matrix changes)
513
514      DO jstp = nit000, nit000 + 2
515         istp = jstp
516         IF ( jstp == nit000 +2 ) istp = nitend
517
518      ! Initialize the reference state
519      avt(:,:,:) = 1.e-1
520      !=============================================================
521      ! 1) dx = ( T ) and dy = ( T )
522      !=============================================================
523
524      !--------------------------------------------------------------------
525      ! Reset the tangent and adjoint variables
526      !--------------------------------------------------------------------
527      zsb_tlin(:,:,:)      = 0.0_wp
528      ztb_tlin(:,:,:)      = 0.0_wp
529      zsa_tlin(:,:,:)      = 0.0_wp
530      zta_tlin(:,:,:)      = 0.0_wp
531      zsa_tlout(:,:,:)     = 0.0_wp
532      zta_tlout(:,:,:)     = 0.0_wp
533      zstrdmp_tlout(:,:,:) = 0.0_wp
534      zttrdmp_tlout(:,:,:) = 0.0_wp
535      zsb_adout(:,:,:)     = 0.0_wp
536      ztb_adout(:,:,:)     = 0.0_wp
537      zsa_adout(:,:,:)     = 0.0_wp
538      zta_adout(:,:,:)     = 0.0_wp
539      zsa_adin(:,:,:)      = 0.0_wp
540      zta_adin(:,:,:)      = 0.0_wp
541      zstrdmp_adin(:,:,:)  = 0.0_wp
542      zttrdmp_adin(:,:,:)  = 0.0_wp
543
544      tsa_ad(:,:,:,:)      = 0.0_wp
545      tsb_ad(:,:,:,:)      = 0.0_wp
546      strdmp_ad(:,:,:)  = 0.0_wp
547      ttrdmp_ad(:,:,:)  = 0.0_wp
548      tsa_tl(:,:,:,:)      = 0.0_wp
549      tsb_tl(:,:,:,:)      = 0.0_wp
550      strdmp_tl(:,:,:)  = 0.0_wp
551      ttrdmp_tl(:,:,:)  = 0.0_wp
552
553      CALL grid_random(  z3r, 'T', 0.0_wp, stds )
554      DO jk = 1, jpk
555         DO jj = nldj, nlej
556            DO ji = nldi, nlei
557               zsb_tlin(ji,jj,jk) = z3r(ji,jj,jk)
558            END DO
559         END DO
560      END DO
561      CALL grid_random(  z3r, 'T', 0.0_wp, stdt )
562      DO jk = 1, jpk
563         DO jj = nldj, nlej
564            DO ji = nldi, nlei
565               ztb_tlin(ji,jj,jk) = z3r(ji,jj,jk)
566            END DO
567         END DO
568      END DO
569
570      CALL grid_random(  z3r, 'T', 0.0_wp, stds )
571      DO jk = 1, jpk
572         DO jj = nldj, nlej
573            DO ji = nldi, nlei
574               zsa_tlin(ji,jj,jk) = z3r(ji,jj,jk)
575            END DO
576         END DO
577      END DO
578
579      CALL grid_random(  z3r, 'T', 0.0_wp, stdt )
580      DO jk = 1, jpk
581         DO jj = nldj, nlej
582            DO ji = nldi, nlei
583               zta_tlin(ji,jj,jk) = z3r(ji,jj,jk)
584            END DO
585         END DO
586      END DO
587
588      tsb_tl(:,:,:,jp_sal) = zsb_tlin(:,:,:)
589      tsb_tl(:,:,:,jp_tem) = ztb_tlin(:,:,:)
590      tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:)
591      tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
592
593      CALL tra_dmp_tan( istp )
594
595      zsa_tlout(:,:,:)     = tsa_tl(:,:,:,jp_sal)
596      zta_tlout(:,:,:)     = tsa_tl(:,:,:,jp_tem)
597      zstrdmp_tlout(:,:,:) = strdmp_tl(:,:,:)
598      zttrdmp_tlout(:,:,:) = ttrdmp_tl(:,:,:)
599
600      !--------------------------------------------------------------------
601      ! Initialize the adjoint variables: dy^* = W dy
602      !--------------------------------------------------------------------
603
604      DO jk = 1, jpk
605        DO jj = nldj, nlej
606           DO ji = nldi, nlei
607              zsa_adin(ji,jj,jk)     = zsa_tlout(ji,jj,jk) &
608                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
609                 &                   * tmask(ji,jj,jk) * wesp_s(jk)
610              zta_adin(ji,jj,jk)     = zta_tlout(ji,jj,jk) &
611                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
612                 &                   * tmask(ji,jj,jk) * wesp_t(jk)
613              zstrdmp_adin(ji,jj,jk) = zstrdmp_tlout(ji,jj,jk) &
614                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
615                 &                   * tmask(ji,jj,jk) * wesp_s(jk)
616              zttrdmp_adin(ji,jj,jk) = zttrdmp_tlout(ji,jj,jk) &
617                 &                   * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
618                 &                   * tmask(ji,jj,jk) * wesp_t(jk)
619            END DO
620         END DO
621      END DO
622
623      !--------------------------------------------------------------------
624      ! Compute the scalar product: ( L dx )^T W dy
625      !--------------------------------------------------------------------
626
627      zsp1_1 = DOT_PRODUCT( zsa_tlout    , zsa_adin     )
628      zsp1_2 = DOT_PRODUCT( zta_tlout    , zta_adin     )
629      zsp1_3 = DOT_PRODUCT( zstrdmp_tlout, zstrdmp_adin )
630      zsp1_4 = DOT_PRODUCT( zttrdmp_tlout, zttrdmp_adin )
631      zsp1   = zsp1_1 + zsp1_2 + zsp1_3 + zsp1_4
632
633      !--------------------------------------------------------------------
634      ! Call the adjoint routine: dx^* = L^T dy^*
635      !--------------------------------------------------------------------
636
637      tsa_ad(:,:,:,jp_sal)     = zsa_adin(:,:,:)
638      tsa_ad(:,:,:,jp_tem)     = zta_adin(:,:,:)
639      strdmp_ad(:,:,:) = zstrdmp_adin(:,:,:)
640      ttrdmp_ad(:,:,:) = zttrdmp_adin(:,:,:)
641
642      CALL tra_dmp_adj( istp )
643
644      zsb_adout(:,:,:) = tsb_ad(:,:,:,jp_sal)
645      ztb_adout(:,:,:) = tsb_ad(:,:,:,jp_tem)
646      zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal)
647      zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
648
649      !--------------------------------------------------------------------
650      ! Compute the scalar product: dx^T L^T W dy
651      !--------------------------------------------------------------------
652
653      zsp2_1 = DOT_PRODUCT( zsb_tlin  , zsb_adout   )
654      zsp2_2 = DOT_PRODUCT( ztb_tlin  , ztb_adout   )
655      zsp2_3 = DOT_PRODUCT( zsa_tlin  , zsa_adout   )
656      zsp2_4 = DOT_PRODUCT( zta_tlin  , zta_adout   )
657
658      zsp2 = zsp2_1 + zsp2_2 + zsp2_3 + zsp2_4
659
660      ! Compare the scalar products
661
662      ! 14 char:'12345678901234'
663      IF ( istp == nit000 ) THEN
664         cl_name = 'tra_dmp_adj T1'
665      ELSEIF ( istp == nit000 +1 ) THEN
666         cl_name = 'tra_dmp_adj T2'
667      ELSEIF ( istp == nitend ) THEN
668         cl_name = 'tra_dmp_adj T3'
669      END IF
670      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
671
672      END DO
673
674      DEALLOCATE( &
675         & zsb_tlin,     &
676         & ztb_tlin,     &
677         & zsa_tlin,     &
678         & zta_tlin,     &
679         & zsa_tlout,    &
680         & zta_tlout,    &
681         & zstrdmp_tlout,&
682         & zttrdmp_tlout,&
683         & zsb_adout,    &
684         & ztb_adout,    &
685         & zsa_adout,    &
686         & zta_adout,    &
687         & zsa_adin,     &
688         & zta_adin,     &
689         & zstrdmp_adin, &
690         & zttrdmp_adin, &
691         & z3r           &
692         & )
693
694   END SUBROUTINE tra_dmp_adj_tst
695
696#else
697   !!----------------------------------------------------------------------
698   !!   Default key                                     NO internal damping
699   !!----------------------------------------------------------------------
700   LOGICAL , PUBLIC, PARAMETER ::   lk_tradmp = .FALSE.    !: internal damping flag
701CONTAINS
702   SUBROUTINE tra_dmp_tan( kt )        ! Empty routine
703      WRITE(*,*) 'tra_dmp_tan: You should not have seen this print! error?', kt
704   END SUBROUTINE tra_dmp_tan
705   SUBROUTINE tra_dmp_adj( kt )        ! Empty routine
706      WRITE(*,*) 'tra_dmp_adj: You should not have seen this print! error?', kt
707   END SUBROUTINE tra_dmp_adj
708   SUBROUTINE tra_dmp_adj_tst( kt )        ! Empty routine
709      WRITE(*,*) 'tra_dmp_adj_tst: You should not have seen this print! error?', kt
710   END SUBROUTINE tra_dmp_adj_tst
711
712   !!======================================================================
713#endif
714END MODULE tradmp_tam
Note: See TracBrowser for help on using the repository browser.