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.
trazdf_imp_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/trazdf_imp_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: 31.0 KB
Line 
1MODULE trazdf_imp_tam
2#ifdef key_tam
3   !!==============================================================================
4   !!                 ***  MODULE  trazdf_imp_tam  ***
5   !! Ocean active tracers:  vertical component of the tracer mixing trend
6   !!                        Tangent and Adjoint Module
7   !!==============================================================================
8   !! History of the direct module:
9   !!            OPA  !  1990-10  (B. Blanke)  Original code
10   !!            7.0  !  1991-11  (G. Madec)
11   !!                 !  1992-06  (M. Imbard) correction on tracer trend loops
12   !!                 !  1996-01  (G. Madec) statement function for e3
13   !!                 !  1997-05  (G. Madec) vertical component of isopycnal
14   !!                 !  1997-07  (G. Madec) geopotential diffusion in s-coord
15   !!                 !  2000-08  (G. Madec) double diffusive mixing
16   !!   NEMO     1.0  !  2002-08  (G. Madec) F90: Free form and module
17   !!            2.0  !  2006-11  (G. Madec) New step reorganisation
18   !!            3.2  !  2009-03  (G. Madec)  heat and salt content trends
19
20   !! History of the T&A module:
21   !!        !  09-01 (A. Vidard) tam of the 06-11 version
22   !!----------------------------------------------------------------------
23   !!   tra_zdf_imp_tan : Update the tracer trend with the diagonal vertical
24   !!                 part of the mixing tensor (tangent).
25   !!   tra_zdf_imp_adj : Update the tracer trend with the diagonal vertical
26   !!                 part of the mixing tensor (adjoint).
27   !!----------------------------------------------------------------------
28   !! * Modules used
29   USE par_kind
30   USE par_oce
31   USE oce_tam
32   USE dom_oce
33   USE oce
34   USE zdf_oce
35   USE ldftra_oce
36   USE zdfddm
37   USE traldf_tam
38   USE in_out_manager
39   USE gridrandom
40   USE dotprodfld
41   USE tstool_tam
42   USE trc_oce
43   USE trc_oce_tam
44   USE ldftra
45   USE lib_mpp
46   USE wrk_nemo
47   USE timing
48   USE ldfslp
49   USE paresp
50
51   IMPLICIT NONE
52   PRIVATE
53
54   !! * Routine accessibility
55   PUBLIC tra_zdf_imp_tan       !  routine called by tra_zdf_tan.F90
56   PUBLIC tra_zdf_imp_adj       !  routine called by tra_zdf_adj.F90
57   PUBLIC tra_zdf_imp_adj_tst   !  routine called by tst.F90
58#if defined key_tst_tlm
59   PUBLIC tra_zdf_imp_tlm_tst   !  routine called by tamtst.F90
60#endif
61
62   !! * Substitutions
63#  include "domzgr_substitute.h90"
64#  include "ldftra_substitute.h90"
65#  include "zdfddm_substitute.h90"
66#  include "vectopt_loop_substitute.h90"
67   !!----------------------------------------------------------------------
68   !!----------------------------------------------------------------------
69   !!  OPA 9.0 , LOCEAN-IPSL (2005)
70   !! $Id: trazdf_imp.F90 1156 2008-06-26 16:06:45Z rblod $
71   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
72   !!----------------------------------------------------------------------
73CONTAINS
74
75   SUBROUTINE tra_zdf_imp_tan( kt, kit000, cdtype, p2dt, ptb_tl, pta_tl, kjpt )
76      !!----------------------------------------------------------------------
77      !!                  ***  ROUTINE tra_zdf_imp_tan  ***
78      !!
79      !! ** Purpose of the direct routine:
80      !!     Compute the trend due to the vertical tracer diffusion
81      !!     including the vertical component of lateral mixing (only for 2nd
82      !!     order operator, for fourth order it is already computed and add
83      !!     to the general trend in traldf.F) and add it to the general trend
84      !!     of the tracer equations.
85      !!
86      !! ** Method of the direct routine :
87      !!      The vertical component of the lateral diffusive trends
88      !!      is provided by a 2nd order operator rotated along neutral or geo-
89      !!      potential surfaces to which an eddy induced advection can be
90      !!      added. It is computed using before fields (forward in time) and
91      !!      isopycnal or geopotential slopes computed in routine ldfslp.
92      !!
93      !!      Second part: vertical trend associated with the vertical physics
94      !!      ===========  (including the vertical flux proportional to dk[t]
95      !!                  associated with the lateral mixing, through the
96      !!                  update of avt)
97      !!      The vertical diffusion of tracers (t & s) is given by:
98      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
99      !!      It is computed using a backward time scheme (t=ta).
100      !!      Surface and bottom boundary conditions: no diffusive flux on
101      !!      both tracers (bottom, applied through the masked field avt).
102      !!      Add this trend to the general trend ta,sa :
103      !!         ta = ta + dz( avt dz(t) )
104      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
105      !!
106      !!      Third part: recover avt resulting from the vertical physics
107      !!      ==========  alone, for further diagnostics (for example to
108      !!                  compute the turbocline depth in zdfmxl.F90).
109      !!         avt = zavt
110      !!         (avs = zavs if lk_zdfddm=T )
111      !!
112      !! ** Remarks on the tangent routine : - key_vvl is not available in tangent yet.
113      !!    Once it will be this routine wil need to be rewritten
114      !!                                     - simplified version, slopes (wslp[ij])
115      !!     assumed to be constant (read from the trajectory). same for av[ts]
116      !!
117      !!---------------------------------------------------------------------
118      !!
119      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
120      INTEGER                              , INTENT(in   ) ::   kit000   ! first time step index
121      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
122      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
123      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
124      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb_tl   ! before and now tracer fields
125      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_tl   ! tracer trend
126      !! * Local declarations
127      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
128      REAL(wp) ::   zavi, zrhstl, znvvl,   & ! temporary scalars
129         ze3tb, ze3tn, ze3ta, zvsfvvl        ! variable vertical scale factors
130      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
131         zwi, zwt, zwd, zws                     ! workspace arrays
132      !!---------------------------------------------------------------------
133      !
134      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_imp_tan')
135      !
136      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zwd, zws )
137      !
138      IF( kt == kit000 ) THEN
139         IF(lwp)WRITE(numout,*)
140         IF(lwp)WRITE(numout,*) 'tra_zdf_imp_tan : implicit vertical mixing on ', cdtype
141         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ '
142      ENDIF
143
144      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
145      ! -------------------
146      ! ... not available  in tangent yet
147      ! II. Vertical trend associated with the vertical physics
148      ! =======================================================
149      !     (including the vertical flux proportional to dk[t] associated
150      !      with the lateral mixing, through the avt update)
151      !     dk[ avt dk[ (t,s) ] ] diffusive trends
152      DO jn = 1, kjpt                                 !  tracer loop  !
153         !                                            ! ============= !
154         !
155         !  Matrix construction
156         ! --------------------
157         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
158         !
159         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR.   &
160            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
161            !
162            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
163            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt  (:,:,2:jpk)
164            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
165            ENDIF
166            zwt(:,:,1) = 0._wp
167            !
168            ! II.0 Matrix construction
169            ! ------------------------
170#if defined key_ldfslp
171            ! update and save of avt (and avs if double diffusive mixing)
172            IF ( ln_traldf_grif ) THEN
173               IF ( lwp ) WRITE(numout, *) 'Griffies operator for lateral tracer diffusion not avaible in TAM yet'
174               CALL abort
175            ELSE IF( l_traldf_rot ) THEN
176               DO jk = 2, jpkm1
177                  DO jj = 2, jpjm1
178                     DO ji = fs_2, fs_jpim1   ! vector opt.
179                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
180                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
181                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
182                     END DO
183                  END DO
184               END DO
185            ENDIF
186#endif
187            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
188            DO jk = 1, jpkm1
189               DO jj = 2, jpjm1
190                  DO ji = fs_2, fs_jpim1   ! vector opt.
191                     ze3ta = 1._wp                                ! after scale factor at T-point
192                     ze3tn = fse3t(ji,jj,jk)                      ! now   scale factor at T-point
193                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
194                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
195                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
196                  END DO
197               END DO
198            END DO
199            !
200            ! II.1. Vertical diffusion on t
201            ! ---------------------------
202            !
203            !! Matrix inversion from the first level
204            !!----------------------------------------------------------------------
205            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
206            !
207            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
208            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
209            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
210            !        (        ...               )( ...  ) ( ...  )
211            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
212            !
213            !   m is decomposed in the product of an upper and lower triangular matrix
214            !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
215            !   The second member is in 2d array zwy
216            !   The solution is in 2d array zwx
217            !   The 3d arry zwt is a work space array
218            !   zwy is used and then used as a work space array : its value is modified!
219
220            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
221            DO jj = 2, jpjm1
222               DO ji = fs_2, fs_jpim1
223                  zwt(ji,jj,1) = zwd(ji,jj,1)
224               END DO
225            END DO
226            DO jk = 2, jpkm1
227               DO jj = 2, jpjm1
228                  DO ji = fs_2, fs_jpim1
229                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
230                  END DO
231               END DO
232            END DO
233         END IF
234         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
235         DO jj = 2, jpjm1
236            DO ji = fs_2, fs_jpim1
237               ze3tb = 1._wp
238               ze3tn = 1._wp
239               pta_tl(ji,jj,1,jn) = ze3tb * ptb_tl(ji,jj,1,jn) + p2dt(1) * ze3tn * pta_tl(ji,jj,1,jn)
240            END DO
241         END DO
242
243         DO jk = 2, jpkm1
244            DO jj = 2, jpjm1
245               DO ji = fs_2, fs_jpim1
246                  ze3tb = 1._wp
247                  ze3tn = 1._wp
248                  zrhstl = ze3tb * ptb_tl(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta_tl(ji,jj,jk,jn)   ! zrhs=right hand side
249                  pta_tl(ji,jj,jk,jn) = zrhstl - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta_tl(ji,jj,jk-1,jn)
250               END DO
251            END DO
252         END DO
253
254         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
255         ! Save the masked temperature after in ta
256         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
257         DO jj = 2, jpjm1
258            DO ji = fs_2, fs_jpim1
259               pta_tl(ji,jj,jpkm1,jn) = pta_tl(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
260            END DO
261         END DO
262         DO jk = jpk-2, 1, -1
263            DO jj = 2, jpjm1
264               DO ji = fs_2, fs_jpim1
265                  pta_tl(ji,jj,jk,jn) = ( pta_tl(ji,jj,jk,jn) - zws(ji,jj,jk) * pta_tl(ji,jj,jk+1,jn) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
266               END DO
267            END DO
268         END DO
269      !                                            ! ================= !
270      END DO                                          !  end tracer loop  !
271      !                                               ! ================= !
272      !
273      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zwd, zws )
274      !
275      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp_tan')
276      !
277   END SUBROUTINE tra_zdf_imp_tan
278   SUBROUTINE tra_zdf_imp_adj( kt, kit000, cdtype, p2dt, ptb_ad, pta_ad, kjpt )
279      !!----------------------------------------------------------------------
280      !!                  ***  ROUTINE tra_zdf_imp_adj  ***
281      !!
282      !! ** Purpose of the direct routine:
283      !!     Compute the trend due to the vertical tracer diffusion
284      !!     including the vertical component of lateral mixing (only for 2nd
285      !!     order operator, for fourth order it is already computed and add
286      !!     to the general trend in traldf.F) and add it to the general trend
287      !!     of the tracer equations.
288      !!
289      !! ** Method of the direct routine :
290      !!      The vertical component of the lateral diffusive trends
291      !!      is provided by a 2nd order operator rotated along neutral or geo-
292      !!      potential surfaces to which an eddy induced advection can be
293      !!      added. It is computed using before fields (forward in time) and
294      !!      isopycnal or geopotential slopes computed in routine ldfslp.
295      !!
296      !!      Second part: vertical trend associated with the vertical physics
297      !!      ===========  (including the vertical flux proportional to dk[t]
298      !!                  associated with the lateral mixing, through the
299      !!                  update of avt)
300      !!      The vertical diffusion of tracers (t & s) is given by:
301      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
302      !!      It is computed using a backward time scheme (t=ta).
303      !!      Surface and bottom boundary conditions: no diffusive flux on
304      !!      both tracers (bottom, applied through the masked field avt).
305      !!      Add this trend to the general trend ta,sa :
306      !!         ta = ta + dz( avt dz(t) )
307      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
308      !!
309      !!      Third part: recover avt resulting from the vertical physics
310      !!      ==========  alone, for further diagnostics (for example to
311      !!                  compute the turbocline depth in zdfmxl.F90).
312      !!         avt = zavt
313      !!         (avs = zavs if lk_zdfddm=T )
314      !!
315      !! ** Remarks on the adjoint routine : - key_vvl is not available in adjoint yet.
316      !!    Once it will be this routine wil need to be rewritten
317      !!                                     - simplified version, slopes (wslp[ij])
318      !!     assumed to be constant (read from the trajectory). same for av[ts]
319      !!
320      !!---------------------------------------------------------------------
321      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
322      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
323      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
324      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
325      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
326      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptb_ad      ! before and now tracer fields
327      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_ad      ! tracer trend
328      !! * Local declarations
329      INTEGER  ::   ji, jj, jk, jn               ! dummy loop indices
330      REAL(wp) ::   zavi, zrhsad, znvvl,   & ! temporary scalars
331         ze3tb, ze3tn, ze3ta, zvsfvvl        ! variable vertical scale factors
332      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
333         zwi, zwt, zws, zwd                     ! workspace arrays
334      !!---------------------------------------------------------------------
335      !
336      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_imp_adj')
337      !
338      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zws, zwd )
339      !
340      IF( kt == nitend ) THEN
341         IF(lwp)WRITE(numout,*)
342         IF(lwp)WRITE(numout,*) 'tra_zdf_imp_adj : implicit vertical mixing on', cdtype
343         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~~~~ '
344      ENDIF
345
346      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
347      ! -------------------
348      ! ... not available  in tangent yet
349      ! II. Vertical trend associated with the vertical physics
350      ! =======================================================
351      !     (including the vertical flux proportional to dk[t] associated
352      !      with the lateral mixing, through the avt update)
353      !     dk[ avt dk[ (t,s) ] ] diffusive trends
354      DO jn = 1, kjpt                                 !  tracer loop  !
355         !                                            ! ============= !
356         !
357         !  Matrix construction
358         ! --------------------
359         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
360         !
361         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR.   &
362            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
363            !
364            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
365            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt  (:,:,2:jpk)
366            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
367            ENDIF
368            zwt(:,:,1) = 0._wp
369
370#if defined key_ldfslp
371      ! update and save of avt (and avs if double diffusive mixing)
372            IF ( ln_traldf_grif ) THEN
373               IF ( lwp ) WRITE(numout, *) 'Griffies operator for lateral tracer diffusion not avaible in TAM yet'
374               CALL abort
375            ELSE IF( l_traldf_rot ) THEN
376               DO jk = 2, jpkm1
377                  DO jj = 2, jpjm1
378                     DO ji = fs_2, fs_jpim1   ! vector opt.
379                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
380                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
381                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
382                     END DO
383                  END DO
384               END DO
385            ENDIF
386#endif
387            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
388            DO jk = 1, jpkm1
389               DO jj = 2, jpjm1
390                  DO ji = fs_2, fs_jpim1   ! vector opt.
391                     ze3ta = 1._wp                                ! after scale factor at T-point
392                     ze3tn = fse3t(ji,jj,jk)                      ! now   scale factor at T-point
393                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
394                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
395                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
396                  END DO
397               END DO
398            END DO
399
400            !! Matrix inversion from the first level
401            !!----------------------------------------------------------------------
402            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
403            !
404            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
405            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
406            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
407            !        (        ...               )( ...  ) ( ...  )
408            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
409            !
410            !   m is decomposed in the product of an upper and lower triangular matrix
411            !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
412            !   The second member is in 2d array zwy
413            !   The solution is in 2d array zwx
414            !   The 3d arry zwt is a work space array
415            !   zwy is used and then used as a work space array : its value is modified!
416            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
417            DO jj = 2, jpjm1
418               DO ji = fs_2, fs_jpim1
419                  zwt(ji,jj,1) = zwd(ji,jj,1)
420               END DO
421            END DO
422            DO jk = 2, jpkm1
423               DO jj = 2, jpjm1
424                  DO ji = fs_2, fs_jpim1
425                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
426                  END DO
427               END DO
428            END DO
429         END IF
430         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
431         ! Save the masked temperature after in ta
432         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
433         DO jk = 1, jpk-2
434            DO jj = 2, jpjm1
435               DO ji = fs_2, fs_jpim1
436                  pta_ad(ji,jj,jk+1,jn) = pta_ad(ji,jj,jk+1,jn) - zws(ji,jj,jk) * pta_ad(ji,jj,jk,jn) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
437                  pta_ad(ji,jj,jk,jn)   = pta_ad(ji,jj,jk,jn) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
438               END DO
439            END DO
440         END DO
441         DO jj = 2, jpjm1
442            DO ji = fs_2, fs_jpim1
443               pta_ad(ji,jj,jpkm1,jn) = pta_ad(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
444            END DO
445         END DO
446         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
447         DO jk = jpkm1, 2, -1
448            DO jj = 2, jpjm1
449               DO ji = fs_2, fs_jpim1
450                  ze3tb = 1._wp
451                  ze3tn = 1._wp
452                  zrhsad = zrhsad + pta_ad(ji,jj,jk,jn)
453                  pta_ad(ji,jj,jk-1,jn) = pta_ad(ji,jj,jk-1,jn) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta_ad(ji,jj,jk,jn)
454                  pta_ad(ji,jj,jk,jn) = 0.0_wp
455                  ptb_ad(ji,jj,jk,jn) = ptb_ad(ji,jj,jk,jn) + ze3tb * zrhsad
456                  pta_ad(ji,jj,jk,jn) = pta_ad(ji,jj,jk,jn) + p2dt(jk) * ze3tn * zrhsad
457                  zrhsad = 0.0_wp
458               END DO
459            END DO
460         END DO
461         DO jj = 2, jpjm1
462            DO ji = fs_2, fs_jpim1
463               ze3tb = 1._wp
464               ze3tn = 1._wp
465               ptb_ad(ji,jj,1,jn) = ptb_ad(ji,jj,1,jn) + ze3tb * pta_ad(ji,jj,1,jn)
466               pta_ad(ji,jj,1,jn) = pta_ad(ji,jj,1,jn) * p2dt(1) * ze3tn
467            END DO
468         END DO
469      END DO
470      !
471      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zws, zwd )
472      !
473      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp_adj')
474      !
475   END SUBROUTINE tra_zdf_imp_adj
476   SUBROUTINE tra_zdf_imp_adj_tst( kumadt )
477      !!-----------------------------------------------------------------------
478      !!
479      !!                  ***  ROUTINE tra_zdf_imp_adj_tst ***
480      !!
481      !! ** Purpose : Test the adjoint routine.
482      !!
483      !! ** Method  : Verify the scalar product
484      !!
485      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
486      !!
487      !!              where  L   = tangent routine
488      !!                     L^T = adjoint routine
489      !!                     W   = diagonal matrix of scale factors
490      !!                     dx  = input perturbation (random field)
491      !!                     dy  = L dx
492      !!
493      !!
494      !! History :
495      !!        ! 08-08 (A. Vidard)
496      !!-----------------------------------------------------------------------
497      !! * Modules used
498
499      !! * Arguments
500      INTEGER, INTENT(IN) :: &
501         & kumadt             ! Output unit
502
503      !! * Local declarations
504      INTEGER ::  &
505         & istp,  &
506         & jstp,  &
507         & ji,    &        ! dummy loop indices
508         & jj,    &
509         & jk
510      REAL(KIND=wp) ::   &
511         & zsp1,         & ! scalar product involving the tangent routine
512         & zsp2            ! scalar product involving the adjoint routine
513      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
514         & zta_tlin ,     & ! Tangent input
515         & ztb_tlin ,     & ! Tangent input
516         & zsa_tlin ,     & ! Tangent input
517         & zsb_tlin ,     & ! Tangent input
518         & zta_tlout,     & ! Tangent output
519         & zsa_tlout,     & ! Tangent output
520         & zta_adin ,     & ! Adjoint input
521         & zsa_adin ,     & ! Adjoint input
522         & zta_adout,     & ! Adjoint output
523         & ztb_adout,     & ! Adjoint output
524         & zsa_adout,     & ! Adjoint output
525         & zsb_adout,     & ! Adjoint output
526         & zr             ! 3D random field
527      CHARACTER(LEN=14) :: cl_name
528      ! Allocate memory
529
530      ALLOCATE( &
531         & zta_tlin( jpi,jpj,jpk),     &
532         & zsa_tlin( jpi,jpj,jpk),     &
533         & ztb_tlin( jpi,jpj,jpk),     &
534         & zsb_tlin( jpi,jpj,jpk),     &
535         & zta_tlout(jpi,jpj,jpk),     &
536         & zsa_tlout(jpi,jpj,jpk),     &
537         & zta_adin( jpi,jpj,jpk),     &
538         & zsa_adin( jpi,jpj,jpk),     &
539         & zta_adout(jpi,jpj,jpk),     &
540         & zsa_adout(jpi,jpj,jpk),     &
541         & ztb_adout(jpi,jpj,jpk),     &
542         & zsb_adout(jpi,jpj,jpk),     &
543         & zr(       jpi,jpj,jpk)      &
544         & )
545      !==================================================================
546      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
547      !    dy = ( hdivb_tl, hdivn_tl )
548      !==================================================================
549
550      ! initialization (normally done in traldf)
551      l_traldf_rot = .TRUE.
552
553      ! Test for time steps nit000 and nit000 + 1 (the matrix changes)
554
555      DO jstp = nit000, nit000 + 2
556         istp = jstp
557         IF ( jstp == nit000+2 ) istp = nitend
558
559      !--------------------------------------------------------------------
560      ! Reset the tangent and adjoint variables
561      !--------------------------------------------------------------------
562          zta_tlin( :,:,:) = 0.0_wp
563          ztb_tlin( :,:,:) = 0.0_wp
564          zsa_tlin( :,:,:) = 0.0_wp
565          zsb_tlin( :,:,:) = 0.0_wp
566          zta_tlout(:,:,:) = 0.0_wp
567          zsa_tlout(:,:,:) = 0.0_wp
568          zta_adin( :,:,:) = 0.0_wp
569          zsa_adin( :,:,:) = 0.0_wp
570          zta_adout(:,:,:) = 0.0_wp
571          zsa_adout(:,:,:) = 0.0_wp
572          ztb_adout(:,:,:) = 0.0_wp
573          zsb_adout(:,:,:) = 0.0_wp
574          zr(       :,:,:) = 0.0_wp
575          tsb_ad(:,:,:,:)     = 0.0_wp
576          tsb_ad(:,:,:,:)     = 0.0_wp
577      !--------------------------------------------------------------------
578      ! Initialize the tangent input with random noise: dx
579      !--------------------------------------------------------------------
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              zta_tlin(ji,jj,jk) = zr(ji,jj,jk)
586            END DO
587         END DO
588      END DO
589      CALL grid_random(  zr, 'T', 0.0_wp, stdt )
590      DO jk = 1, jpk
591        DO jj = nldj, nlej
592           DO ji = nldi, nlei
593              ztb_tlin(ji,jj,jk) = zr(ji,jj,jk)
594            END DO
595         END DO
596      END DO
597      CALL grid_random(  zr, 'T', 0.0_wp, stds )
598      DO jk = 1, jpk
599        DO jj = nldj, nlej
600           DO ji = nldi, nlei
601              zsa_tlin(ji,jj,jk) = zr(ji,jj,jk)
602            END DO
603         END DO
604      END DO
605      CALL grid_random(  zr, 'T', 0.0_wp, stds )
606      DO jk = 1, jpk
607        DO jj = nldj, nlej
608           DO ji = nldi, nlei
609              zsb_tlin(ji,jj,jk) = zr(ji,jj,jk)
610            END DO
611         END DO
612      END DO
613
614
615      tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
616      tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:)
617      tsb_tl(:,:,:,jp_tem) = ztb_tlin(:,:,:)
618      tsb_tl(:,:,:,jp_sal) = zsb_tlin(:,:,:)
619      CALL tra_zdf_imp_tan ( istp, nit000, 'TRA', r2dtra, tsb_tl, tsa_tl, jpts )
620      zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
621      zsa_tlout(:,:,:) = tsa_tl(:,:,:,jp_sal)
622
623      !--------------------------------------------------------------------
624      ! Initialize the adjoint variables: dy^* = W dy
625      !--------------------------------------------------------------------
626
627      DO jk = 1, jpk
628        DO jj = nldj, nlej
629           DO ji = nldi, nlei
630              zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
631                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
632                 &               * tmask(ji,jj,jk) * wesp_t(jk)
633              zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) &
634                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
635                 &               * tmask(ji,jj,jk) * wesp_s(jk)
636            END DO
637         END DO
638      END DO
639      !--------------------------------------------------------------------
640      ! Compute the scalar product: ( L dx )^T W dy
641      !--------------------------------------------------------------------
642
643      zsp1 = DOT_PRODUCT( zta_tlout, zta_adin ) &
644         & + DOT_PRODUCT( zsa_tlout, zsa_adin )
645
646      !--------------------------------------------------------------------
647      ! Call the adjoint routine: dx^* = L^T dy^*
648      !--------------------------------------------------------------------
649
650      tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
651      tsa_ad(:,:,:,jp_sal) = zsa_adin(:,:,:)
652
653      CALL tra_zdf_imp_adj ( istp, nit000, 'TRA', r2dtra, tsb_ad, tsa_ad, jpts )
654
655      zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
656      zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal)
657      ztb_adout(:,:,:) = tsb_ad(:,:,:,jp_tem)
658      zsb_adout(:,:,:) = tsb_ad(:,:,:,jp_sal)
659      zsp2 = DOT_PRODUCT( zta_tlin, zta_adout ) &
660         & + DOT_PRODUCT( zsa_tlin, zsa_adout ) &
661         & + DOT_PRODUCT( ztb_tlin, ztb_adout ) &
662         & + DOT_PRODUCT( zsb_tlin, zsb_adout )
663
664      ! 14 char:'12345678901234'
665      IF ( istp == nit000 ) THEN
666         cl_name = 'trazdfimpadjT1'
667      ELSEIF ( istp == nit000 +1 ) THEN
668         cl_name = 'trazdfimpadjT2'
669      ELSEIF ( istp == nitend ) THEN
670         cl_name = 'trazdfimpadjT3'
671      END IF
672      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
673
674      END DO
675
676      DEALLOCATE(   &
677         & zta_tlin,  &
678         & ztb_tlin,  &
679         & zsa_tlin,  &
680         & zsb_tlin,  &
681         & zta_tlout, &
682         & zsa_tlout, &
683         & zta_adin,  &
684         & zsa_adin,  &
685         & zta_adout, &
686         & ztb_adout, &
687         & zsa_adout, &
688         & zsb_adout, &
689         & zr       &
690         & )
691
692
693
694   END SUBROUTINE tra_zdf_imp_adj_tst
695#endif
696END MODULE trazdf_imp_tam
Note: See TracBrowser for help on using the repository browser.