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 @ 3627

Last change on this file since 3627 was 3627, checked in by rblod, 10 years ago

Correct long lines in TAM 4.3 see ticket #1010

  • Property svn:executable set to *
File size: 31.1 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) )   &
266                                  &   / zwt(ji,jj,jk) * tmask(ji,jj,jk)
267               END DO
268            END DO
269         END DO
270      !                                            ! ================= !
271      END DO                                          !  end tracer loop  !
272      !                                               ! ================= !
273      !
274      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zwd, zws )
275      !
276      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp_tan')
277      !
278   END SUBROUTINE tra_zdf_imp_tan
279   SUBROUTINE tra_zdf_imp_adj( kt, kit000, cdtype, p2dt, ptb_ad, pta_ad, kjpt )
280      !!----------------------------------------------------------------------
281      !!                  ***  ROUTINE tra_zdf_imp_adj  ***
282      !!
283      !! ** Purpose of the direct routine:
284      !!     Compute the trend due to the vertical tracer diffusion
285      !!     including the vertical component of lateral mixing (only for 2nd
286      !!     order operator, for fourth order it is already computed and add
287      !!     to the general trend in traldf.F) and add it to the general trend
288      !!     of the tracer equations.
289      !!
290      !! ** Method of the direct routine :
291      !!      The vertical component of the lateral diffusive trends
292      !!      is provided by a 2nd order operator rotated along neutral or geo-
293      !!      potential surfaces to which an eddy induced advection can be
294      !!      added. It is computed using before fields (forward in time) and
295      !!      isopycnal or geopotential slopes computed in routine ldfslp.
296      !!
297      !!      Second part: vertical trend associated with the vertical physics
298      !!      ===========  (including the vertical flux proportional to dk[t]
299      !!                  associated with the lateral mixing, through the
300      !!                  update of avt)
301      !!      The vertical diffusion of tracers (t & s) is given by:
302      !!             difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) )
303      !!      It is computed using a backward time scheme (t=ta).
304      !!      Surface and bottom boundary conditions: no diffusive flux on
305      !!      both tracers (bottom, applied through the masked field avt).
306      !!      Add this trend to the general trend ta,sa :
307      !!         ta = ta + dz( avt dz(t) )
308      !!         (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T )
309      !!
310      !!      Third part: recover avt resulting from the vertical physics
311      !!      ==========  alone, for further diagnostics (for example to
312      !!                  compute the turbocline depth in zdfmxl.F90).
313      !!         avt = zavt
314      !!         (avs = zavs if lk_zdfddm=T )
315      !!
316      !! ** Remarks on the adjoint routine : - key_vvl is not available in adjoint yet.
317      !!    Once it will be this routine wil need to be rewritten
318      !!                                     - simplified version, slopes (wslp[ij])
319      !!     assumed to be constant (read from the trajectory). same for av[ts]
320      !!
321      !!---------------------------------------------------------------------
322      INTEGER                              , INTENT(in   ) ::   kt       ! ocean time-step index
323      INTEGER                              , INTENT(in   ) ::   kit000          ! first time step index
324      CHARACTER(len=3)                     , INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
325      INTEGER                              , INTENT(in   ) ::   kjpt     ! number of tracers
326      REAL(wp), DIMENSION(        jpk     ), INTENT(in   ) ::   p2dt     ! vertical profile of tracer time-step
327      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptb_ad      ! before and now tracer fields
328      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_ad      ! tracer trend
329      !! * Local declarations
330      INTEGER  ::   ji, jj, jk, jn               ! dummy loop indices
331      REAL(wp) ::   zavi, zrhsad, znvvl,   & ! temporary scalars
332         ze3tb, ze3tn, ze3ta, zvsfvvl        ! variable vertical scale factors
333      REAL(wp), POINTER, DIMENSION(:,:,:) ::   &
334         zwi, zwt, zws, zwd                     ! workspace arrays
335      !!---------------------------------------------------------------------
336      !
337      IF( nn_timing == 1 )  CALL timing_start('tra_zdf_imp_adj')
338      !
339      CALL wrk_alloc( jpi, jpj, jpk, zwi, zwt, zws, zwd )
340      !
341      IF( kt == nitend ) THEN
342         IF(lwp)WRITE(numout,*)
343         IF(lwp)WRITE(numout,*) 'tra_zdf_imp_adj : implicit vertical mixing on', cdtype
344         IF(lwp)WRITE(numout,*) '~~~~~~~~~~~~~~~ '
345      ENDIF
346
347      ! I.1 Variable volume : to take into account vertical variable vertical scale factors
348      ! -------------------
349      ! ... not available  in tangent yet
350      ! II. Vertical trend associated with the vertical physics
351      ! =======================================================
352      !     (including the vertical flux proportional to dk[t] associated
353      !      with the lateral mixing, through the avt update)
354      !     dk[ avt dk[ (t,s) ] ] diffusive trends
355      DO jn = 1, kjpt                                 !  tracer loop  !
356         !                                            ! ============= !
357         !
358         !  Matrix construction
359         ! --------------------
360         ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
361         !
362         IF(  ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR.   &
363            & ( cdtype == 'TRC' .AND. jn == 1 )  )  THEN
364            !
365            ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
366            IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN   ;   zwt(:,:,2:jpk) = avt  (:,:,2:jpk)
367            ELSE                                            ;   zwt(:,:,2:jpk) = fsavs(:,:,2:jpk)
368            ENDIF
369            zwt(:,:,1) = 0._wp
370
371#if defined key_ldfslp
372      ! update and save of avt (and avs if double diffusive mixing)
373            IF ( ln_traldf_grif ) THEN
374               IF ( lwp ) WRITE(numout, *) 'Griffies operator for lateral tracer diffusion not avaible in TAM yet'
375               CALL abort
376            ELSE IF( l_traldf_rot ) THEN
377               DO jk = 2, jpkm1
378                  DO jj = 2, jpjm1
379                     DO ji = fs_2, fs_jpim1   ! vector opt.
380                        zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk)                       &   ! vertical mixing coef. due to lateral mixing
381                           & * (  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   &
382                           &    + wslpj(ji,jj,jk) * wslpj(ji,jj,jk)  )
383                     END DO
384                  END DO
385               END DO
386            ENDIF
387#endif
388            ! Diagonal, inferior, superior  (including the bottom boundary condition via avt masked)
389            DO jk = 1, jpkm1
390               DO jj = 2, jpjm1
391                  DO ji = fs_2, fs_jpim1   ! vector opt.
392                     ze3ta = 1._wp                                ! after scale factor at T-point
393                     ze3tn = fse3t(ji,jj,jk)                      ! now   scale factor at T-point
394                     zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk  ) / ( ze3tn * fse3w(ji,jj,jk  ) )
395                     zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )
396                     zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)
397                  END DO
398               END DO
399            END DO
400
401            !! Matrix inversion from the first level
402            !!----------------------------------------------------------------------
403            !   solve m.x = y  where m is a tri diagonal matrix ( jpk*jpk )
404            !
405            !        ( zwd1 zws1   0    0    0  )( zwx1 ) ( zwy1 )
406            !        ( zwi2 zwd2 zws2   0    0  )( zwx2 ) ( zwy2 )
407            !        (  0   zwi3 zwd3 zws3   0  )( zwx3 )=( zwy3 )
408            !        (        ...               )( ...  ) ( ...  )
409            !        (  0    0    0   zwik zwdk )( zwxk ) ( zwyk )
410            !
411            !   m is decomposed in the product of an upper and lower triangular matrix
412            !   The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
413            !   The second member is in 2d array zwy
414            !   The solution is in 2d array zwx
415            !   The 3d arry zwt is a work space array
416            !   zwy is used and then used as a work space array : its value is modified!
417            ! first recurrence:   Tk = Dk - Ik Sk-1 / Tk-1   (increasing k)
418            DO jj = 2, jpjm1
419               DO ji = fs_2, fs_jpim1
420                  zwt(ji,jj,1) = zwd(ji,jj,1)
421               END DO
422            END DO
423            DO jk = 2, jpkm1
424               DO jj = 2, jpjm1
425                  DO ji = fs_2, fs_jpim1
426                     zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1)  /zwt(ji,jj,jk-1)
427                  END DO
428               END DO
429            END DO
430         END IF
431         ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
432         ! Save the masked temperature after in ta
433         ! (c a u t i o n:  temperature not its trend, Leap-frog scheme done it will not be done in tranxt)
434         DO jk = 1, jpk-2
435            DO jj = 2, jpjm1
436               DO ji = fs_2, fs_jpim1
437                  pta_ad(ji,jj,jk+1,jn) = pta_ad(ji,jj,jk+1,jn) - zws(ji,jj,jk) * pta_ad(ji,jj,jk,jn)   &
438                                    &   / zwt(ji,jj,jk) * tmask(ji,jj,jk)
439                  pta_ad(ji,jj,jk,jn)   = pta_ad(ji,jj,jk,jn) / zwt(ji,jj,jk) * tmask(ji,jj,jk)
440               END DO
441            END DO
442         END DO
443         DO jj = 2, jpjm1
444            DO ji = fs_2, fs_jpim1
445               pta_ad(ji,jj,jpkm1,jn) = pta_ad(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
446            END DO
447         END DO
448         ! second recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
449         DO jk = jpkm1, 2, -1
450            DO jj = 2, jpjm1
451               DO ji = fs_2, fs_jpim1
452                  ze3tb = 1._wp
453                  ze3tn = 1._wp
454                  zrhsad = zrhsad + pta_ad(ji,jj,jk,jn)
455                  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)
456                  pta_ad(ji,jj,jk,jn) = 0.0_wp
457                  ptb_ad(ji,jj,jk,jn) = ptb_ad(ji,jj,jk,jn) + ze3tb * zrhsad
458                  pta_ad(ji,jj,jk,jn) = pta_ad(ji,jj,jk,jn) + p2dt(jk) * ze3tn * zrhsad
459                  zrhsad = 0.0_wp
460               END DO
461            END DO
462         END DO
463         DO jj = 2, jpjm1
464            DO ji = fs_2, fs_jpim1
465               ze3tb = 1._wp
466               ze3tn = 1._wp
467               ptb_ad(ji,jj,1,jn) = ptb_ad(ji,jj,1,jn) + ze3tb * pta_ad(ji,jj,1,jn)
468               pta_ad(ji,jj,1,jn) = pta_ad(ji,jj,1,jn) * p2dt(1) * ze3tn
469            END DO
470         END DO
471      END DO
472      !
473      CALL wrk_dealloc( jpi, jpj, jpk, zwi, zwt, zws, zwd )
474      !
475      IF( nn_timing == 1 )  CALL timing_stop('tra_zdf_imp_adj')
476      !
477   END SUBROUTINE tra_zdf_imp_adj
478   SUBROUTINE tra_zdf_imp_adj_tst( kumadt )
479      !!-----------------------------------------------------------------------
480      !!
481      !!                  ***  ROUTINE tra_zdf_imp_adj_tst ***
482      !!
483      !! ** Purpose : Test the adjoint routine.
484      !!
485      !! ** Method  : Verify the scalar product
486      !!
487      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
488      !!
489      !!              where  L   = tangent routine
490      !!                     L^T = adjoint routine
491      !!                     W   = diagonal matrix of scale factors
492      !!                     dx  = input perturbation (random field)
493      !!                     dy  = L dx
494      !!
495      !!
496      !! History :
497      !!        ! 08-08 (A. Vidard)
498      !!-----------------------------------------------------------------------
499      !! * Modules used
500
501      !! * Arguments
502      INTEGER, INTENT(IN) :: &
503         & kumadt             ! Output unit
504
505      !! * Local declarations
506      INTEGER ::  &
507         & istp,  &
508         & jstp,  &
509         & ji,    &        ! dummy loop indices
510         & jj,    &
511         & jk
512      REAL(KIND=wp) ::   &
513         & zsp1,         & ! scalar product involving the tangent routine
514         & zsp2            ! scalar product involving the adjoint routine
515      REAL(KIND=wp), DIMENSION(:,:,:), ALLOCATABLE :: &
516         & zta_tlin ,     & ! Tangent input
517         & ztb_tlin ,     & ! Tangent input
518         & zsa_tlin ,     & ! Tangent input
519         & zsb_tlin ,     & ! Tangent input
520         & zta_tlout,     & ! Tangent output
521         & zsa_tlout,     & ! Tangent output
522         & zta_adin ,     & ! Adjoint input
523         & zsa_adin ,     & ! Adjoint input
524         & zta_adout,     & ! Adjoint output
525         & ztb_adout,     & ! Adjoint output
526         & zsa_adout,     & ! Adjoint output
527         & zsb_adout,     & ! Adjoint output
528         & zr             ! 3D random field
529      CHARACTER(LEN=14) :: cl_name
530      ! Allocate memory
531
532      ALLOCATE( &
533         & zta_tlin( jpi,jpj,jpk),     &
534         & zsa_tlin( jpi,jpj,jpk),     &
535         & ztb_tlin( jpi,jpj,jpk),     &
536         & zsb_tlin( jpi,jpj,jpk),     &
537         & zta_tlout(jpi,jpj,jpk),     &
538         & zsa_tlout(jpi,jpj,jpk),     &
539         & zta_adin( jpi,jpj,jpk),     &
540         & zsa_adin( jpi,jpj,jpk),     &
541         & zta_adout(jpi,jpj,jpk),     &
542         & zsa_adout(jpi,jpj,jpk),     &
543         & ztb_adout(jpi,jpj,jpk),     &
544         & zsb_adout(jpi,jpj,jpk),     &
545         & zr(       jpi,jpj,jpk)      &
546         & )
547      !==================================================================
548      ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
549      !    dy = ( hdivb_tl, hdivn_tl )
550      !==================================================================
551
552      ! initialization (normally done in traldf)
553      l_traldf_rot = .TRUE.
554
555      ! Test for time steps nit000 and nit000 + 1 (the matrix changes)
556
557      DO jstp = nit000, nit000 + 2
558         istp = jstp
559         IF ( jstp == nit000+2 ) istp = nitend
560
561      !--------------------------------------------------------------------
562      ! Reset the tangent and adjoint variables
563      !--------------------------------------------------------------------
564          zta_tlin( :,:,:) = 0.0_wp
565          ztb_tlin( :,:,:) = 0.0_wp
566          zsa_tlin( :,:,:) = 0.0_wp
567          zsb_tlin( :,:,:) = 0.0_wp
568          zta_tlout(:,:,:) = 0.0_wp
569          zsa_tlout(:,:,:) = 0.0_wp
570          zta_adin( :,:,:) = 0.0_wp
571          zsa_adin( :,:,:) = 0.0_wp
572          zta_adout(:,:,:) = 0.0_wp
573          zsa_adout(:,:,:) = 0.0_wp
574          ztb_adout(:,:,:) = 0.0_wp
575          zsb_adout(:,:,:) = 0.0_wp
576          zr(       :,:,:) = 0.0_wp
577          tsb_ad(:,:,:,:)     = 0.0_wp
578          tsb_ad(:,:,:,:)     = 0.0_wp
579      !--------------------------------------------------------------------
580      ! Initialize the tangent input with random noise: dx
581      !--------------------------------------------------------------------
582
583      CALL grid_random(  zr, 'T', 0.0_wp, stdt )
584      DO jk = 1, jpk
585        DO jj = nldj, nlej
586           DO ji = nldi, nlei
587              zta_tlin(ji,jj,jk) = zr(ji,jj,jk)
588            END DO
589         END DO
590      END DO
591      CALL grid_random(  zr, 'T', 0.0_wp, stdt )
592      DO jk = 1, jpk
593        DO jj = nldj, nlej
594           DO ji = nldi, nlei
595              ztb_tlin(ji,jj,jk) = zr(ji,jj,jk)
596            END DO
597         END DO
598      END DO
599      CALL grid_random(  zr, 'T', 0.0_wp, stds )
600      DO jk = 1, jpk
601        DO jj = nldj, nlej
602           DO ji = nldi, nlei
603              zsa_tlin(ji,jj,jk) = zr(ji,jj,jk)
604            END DO
605         END DO
606      END DO
607      CALL grid_random(  zr, 'T', 0.0_wp, stds )
608      DO jk = 1, jpk
609        DO jj = nldj, nlej
610           DO ji = nldi, nlei
611              zsb_tlin(ji,jj,jk) = zr(ji,jj,jk)
612            END DO
613         END DO
614      END DO
615
616
617      tsa_tl(:,:,:,jp_tem) = zta_tlin(:,:,:)
618      tsa_tl(:,:,:,jp_sal) = zsa_tlin(:,:,:)
619      tsb_tl(:,:,:,jp_tem) = ztb_tlin(:,:,:)
620      tsb_tl(:,:,:,jp_sal) = zsb_tlin(:,:,:)
621      CALL tra_zdf_imp_tan ( istp, nit000, 'TRA', r2dtra, tsb_tl, tsa_tl, jpts )
622      zta_tlout(:,:,:) = tsa_tl(:,:,:,jp_tem)
623      zsa_tlout(:,:,:) = tsa_tl(:,:,:,jp_sal)
624
625      !--------------------------------------------------------------------
626      ! Initialize the adjoint variables: dy^* = W dy
627      !--------------------------------------------------------------------
628
629      DO jk = 1, jpk
630        DO jj = nldj, nlej
631           DO ji = nldi, nlei
632              zta_adin(ji,jj,jk) = zta_tlout(ji,jj,jk) &
633                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
634                 &               * tmask(ji,jj,jk) * wesp_t(jk)
635              zsa_adin(ji,jj,jk) = zsa_tlout(ji,jj,jk) &
636                 &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
637                 &               * tmask(ji,jj,jk) * wesp_s(jk)
638            END DO
639         END DO
640      END DO
641      !--------------------------------------------------------------------
642      ! Compute the scalar product: ( L dx )^T W dy
643      !--------------------------------------------------------------------
644
645      zsp1 = DOT_PRODUCT( zta_tlout, zta_adin ) &
646         & + DOT_PRODUCT( zsa_tlout, zsa_adin )
647
648      !--------------------------------------------------------------------
649      ! Call the adjoint routine: dx^* = L^T dy^*
650      !--------------------------------------------------------------------
651
652      tsa_ad(:,:,:,jp_tem) = zta_adin(:,:,:)
653      tsa_ad(:,:,:,jp_sal) = zsa_adin(:,:,:)
654
655      CALL tra_zdf_imp_adj ( istp, nit000, 'TRA', r2dtra, tsb_ad, tsa_ad, jpts )
656
657      zta_adout(:,:,:) = tsa_ad(:,:,:,jp_tem)
658      zsa_adout(:,:,:) = tsa_ad(:,:,:,jp_sal)
659      ztb_adout(:,:,:) = tsb_ad(:,:,:,jp_tem)
660      zsb_adout(:,:,:) = tsb_ad(:,:,:,jp_sal)
661      zsp2 = DOT_PRODUCT( zta_tlin, zta_adout ) &
662         & + DOT_PRODUCT( zsa_tlin, zsa_adout ) &
663         & + DOT_PRODUCT( ztb_tlin, ztb_adout ) &
664         & + DOT_PRODUCT( zsb_tlin, zsb_adout )
665
666      ! 14 char:'12345678901234'
667      IF ( istp == nit000 ) THEN
668         cl_name = 'trazdfimpadjT1'
669      ELSEIF ( istp == nit000 +1 ) THEN
670         cl_name = 'trazdfimpadjT2'
671      ELSEIF ( istp == nitend ) THEN
672         cl_name = 'trazdfimpadjT3'
673      END IF
674      CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
675
676      END DO
677
678      DEALLOCATE(   &
679         & zta_tlin,  &
680         & ztb_tlin,  &
681         & zsa_tlin,  &
682         & zsb_tlin,  &
683         & zta_tlout, &
684         & zsa_tlout, &
685         & zta_adin,  &
686         & zsa_adin,  &
687         & zta_adout, &
688         & ztb_adout, &
689         & zsa_adout, &
690         & zsb_adout, &
691         & zr       &
692         & )
693
694
695
696   END SUBROUTINE tra_zdf_imp_adj_tst
697#endif
698END MODULE trazdf_imp_tam
Note: See TracBrowser for help on using the repository browser.