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.
trabbl_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/trabbl_tam.F90 @ 3640

Last change on this file since 3640 was 3640, checked in by pabouttier, 10 years ago

Missing allocation/deallocation in TAM routines - See ticket #1013

  • Property svn:executable set to *
File size: 86.1 KB
Line 
1MODULE trabbl_tam
2   !!==============================================================================
3   !!                       ***  MODULE  trabbl  ***
4   !! Ocean physics :  advective and/or diffusive bottom boundary layer scheme
5   !!==============================================================================
6   !! History :  OPA  ! 1996-06  (L. Mortier)  Original code
7   !!            8.0  ! 1997-11  (G. Madec)    Optimization
8   !!   NEMO     1.0  ! 2002-08  (G. Madec)  free form + modules
9   !!             -   ! 2004-01  (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl
10   !!            3.3  ! 2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization
11   !!             -   ! 2010-04  (G. Madec)  Campin & Goosse advective bbl
12   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
13   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
14   !! History of the T&A module
15   !!   NEMO     3.2  ! 2011-02  (A. Vidard) Original version
16   !!            3.4  ! 2012-09  (A. Vidard) Update to 3.4
17   !!
18   !!----------------------------------------------------------------------
19#if   defined key_trabbl   ||   defined key_esopa
20   !!----------------------------------------------------------------------
21   !!   'key_trabbl'   or                             bottom boundary layer
22   !!----------------------------------------------------------------------
23   !!   tra_bbl_alloc : allocate trabbl arrays
24   !!   tra_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
25   !!   tra_bbl_dif   : generic routine to compute bbl diffusive trend
26   !!   tra_bbl_adv   : generic routine to compute bbl advective trend
27   !!   bbl           : computation of bbl diffu. flux coef. & transport in bottom boundary layer
28   !!   tra_bbl_init  : initialization, namelist read, parameters control
29   !!----------------------------------------------------------------------
30   USE oce            ! ocean dynamics and active tracers
31   USE oce_tam
32   USE dom_oce        ! ocean space and time domain
33   USE phycst         ! physical constant
34   USE eosbn2         ! equation of state
35   USE iom            ! IOM server
36   USE in_out_manager ! I/O manager
37   USE lbclnk         ! ocean lateral boundary conditions
38   USE prtctl         ! Print control
39   USE wrk_nemo       ! Memory Allocation
40   USE timing         ! Timing
41   USE trabbl
42   USE gridrandom
43   USE dotprodfld
44   USE tstool_tam
45
46   IMPLICIT NONE
47   PRIVATE
48
49   PUBLIC   tra_bbl_tan       !  routine called by step.F90
50   PUBLIC   tra_bbl_init_tam  !  routine called by opa.F90
51   PUBLIC   tra_bbl_dif_tan   !  routine called by trcbbl.F90
52   PUBLIC   tra_bbl_adv_tan   !  -          -          -              -
53   PUBLIC   bbl_tan           !  routine called by trcbbl.F90 and dtadyn.F90
54   PUBLIC   tra_bbl_adj       !  routine called by step.F90
55   PUBLIC   tra_bbl_dif_adj   !  routine called by trcbbl.F90
56   PUBLIC   tra_bbl_adv_adj   !  -          -          -              -
57   PUBLIC   bbl_adj           !  routine called by trcbbl.F90 and dtadyn.F90
58   PUBLIC   tra_bbl_adj_tst   !  routine called by tamtst
59
60   REAL(WP), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: utr_bbl_tl, vtr_bbl_tl
61   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: ahu_bbl_tl, ahv_bbl_tl
62   REAL(WP), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: utr_bbl_ad, vtr_bbl_ad
63   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: ahu_bbl_ad, ahv_bbl_ad
64
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ahu_bbl_0_tl, ahv_bbl_0_tl
66   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ahu_bbl_0_ad, ahv_bbl_0_ad
67
68   LOGICAL, PRIVATE :: ll_alloctl = .FALSE., ll_allocad = .FALSE.
69
70   !! * Substitutions
71#  include "domzgr_substitute.h90"
72#  include "vectopt_loop_substitute.h90"
73   !!----------------------------------------------------------------------
74   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
75   !! $Id$
76   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
77   !!----------------------------------------------------------------------
78CONTAINS
79
80   INTEGER FUNCTION tra_bbl_alloc_tam( kmode )
81      !!----------------------------------------------------------------------
82      !!                ***  FUNCTION tra_bbl_alloc_tam  ***
83      !!----------------------------------------------------------------------
84      INTEGER, OPTIONAL :: kmode
85      INTEGER, DIMENSION(2) :: ierr
86      INTEGER :: jmode
87
88      IF ( PRESENT( kmode ) ) THEN
89         jmode = kmode
90      ELSE
91         jmode = 0
92      END IF
93      ierr(:) = 0.0_wp
94
95      IF ( ( jmode == 0 ) .OR. ( jmode == 1 ) .AND. ( .NOT. ll_alloctl ) ) THEN
96         ALLOCATE( utr_bbl_tl  (jpi,jpj), vtr_bbl_tl  (jpi,jpj), &
97            &      ahu_bbl_tl  (jpi,jpj), ahv_bbl_tl  (jpi,jpj), &
98            &      ahu_bbl_0_tl  (jpi,jpj), ahv_bbl_0_tl  (jpi,jpj), &
99            &      STAT= ierr(1) )
100         ll_alloctl = .TRUE.
101      END IF
102         !
103      IF ( ( jmode == 0 ) .OR. ( jmode == 2 ) .AND. ( .NOT. ll_allocad ) ) THEN
104         ALLOCATE( utr_bbl_ad  (jpi,jpj), vtr_bbl_ad  (jpi,jpj), &
105            &      ahu_bbl_ad  (jpi,jpj), ahv_bbl_ad  (jpi,jpj), &
106            &      ahu_bbl_0_ad  (jpi,jpj), ahv_bbl_0_ad  (jpi,jpj), &
107            &      STAT= ierr(2) )
108         ll_allocad = .TRUE.
109      END IF
110      tra_bbl_alloc_tam = SUM(ierr)
111         !
112      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc_tam )
113      IF( tra_bbl_alloc_tam > 0 )   CALL ctl_warn('tra_bbl_alloc_tam: allocation of arrays failed.')
114   END FUNCTION tra_bbl_alloc_tam
115
116
117   INTEGER FUNCTION tra_bbl_dealloc_tam( kmode )
118      !!----------------------------------------------------------------------
119      !!                ***  FUNCTION tra_bbl_dealloc  ***
120      !!----------------------------------------------------------------------
121      INTEGER, OPTIONAL :: kmode
122      INTEGER, DIMENSION(2) :: ierr
123
124      IF ( .NOT. PRESENT( kmode ) ) kmode=0
125      ierr(:) = 0.0_wp
126
127      IF ( ( kmode == 0 ) .OR. ( kmode == 1 ) .AND. ( ll_alloctl ) ) THEN
128         DEALLOCATE( utr_bbl_tl, vtr_bbl_tl, &
129            &      ahu_bbl_0_tl, ahv_bbl_0_tl, &
130            &      STAT= ierr(1) )
131         ll_alloctl = .FALSE.
132      END IF
133         !
134      IF ( ( kmode == 0 ) .OR. ( kmode == 1 ) .AND. ( ll_allocad ) ) THEN
135         DEALLOCATE( utr_bbl_ad, vtr_bbl_ad, &
136            &      ahu_bbl_ad, ahv_bbl_ad, &
137            &      ahu_bbl_0_ad, ahv_bbl_0_ad, &
138            &      STAT= ierr(2) )
139         ll_allocad = .FALSE.
140      END IF
141      tra_bbl_dealloc_tam = SUM(ierr)
142         !
143      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_dealloc_tam )
144      IF( tra_bbl_dealloc_tam > 0 )   CALL ctl_warn('tra_bbl_dealloc_tam: allocation of arrays failed.')
145   END FUNCTION tra_bbl_dealloc_tam
146
147
148   SUBROUTINE tra_bbl_tan( kt )
149      !!----------------------------------------------------------------------
150      !!                  ***  ROUTINE bbl_tan  ***
151      !!
152      !! ** Purpose :   Compute the before tracer (t & s) trend associated
153      !!              with the bottom boundary layer and add it to the general
154      !!              trend of tracer equations.
155      !!
156      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
157      !!              diffusive and/or advective contribution to the tracer trend
158      !!              is added to the general tracer trend
159      !!----------------------------------------------------------------------
160      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
161      !!
162      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdttl, ztrdstl
163      !!----------------------------------------------------------------------
164      !
165      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_tan')
166      !
167      IF( l_bbl )  CALL bbl_tan( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
168
169      IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl
170         !
171         CALL tra_bbl_dif_tan( tsb, tsb_tl, tsa_tl, jpts )
172         !
173      END IF
174
175      IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl
176         !
177         CALL tra_bbl_adv_tan( tsb, tsb_tl, tsa_tl, jpts )
178         !
179      END IF
180      !
181      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_tan')
182      !
183   END SUBROUTINE tra_bbl_tan
184
185
186   SUBROUTINE tra_bbl_dif_tan( ptb, ptb_tl, pta_tl, kjpt )
187      !!----------------------------------------------------------------------
188      !!                  ***  ROUTINE tra_bbl_dif_tan  ***
189      !!
190      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
191      !!                advection terms.
192      !!
193      !! ** Method  :
194      !!        * diffusive bbl (nn_bbl_ldf=1) :
195      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
196      !!      along bottom slope gradient) an additional lateral 2nd order
197      !!      diffusion along the bottom slope is added to the general
198      !!      tracer trend, otherwise the additional trend is set to 0.
199      !!      A typical value of ahbt is 2000 m2/s (equivalent to
200      !!      a downslope velocity of 20 cm/s if the condition for slope
201      !!      convection is satified)
202      !!
203      !! ** Action  :   pta   increased by the bbl diffusive trend
204      !!
205      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
206      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
207      !!----------------------------------------------------------------------
208      !
209      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
210      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
211      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb_tl    ! before and now tracer fields
212      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_tl    ! tracer trend
213      !
214      INTEGER  ::   ji, jj, jn   ! dummy loop indices
215      INTEGER  ::   ik           ! local integers
216      REAL(wp) ::   zbtr         ! local scalars
217      REAL(wp), POINTER, DIMENSION(:,:) :: zptb, zptbtl
218      !!----------------------------------------------------------------------
219      !
220      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif_tan')
221      !
222      CALL wrk_alloc( jpi, jpj, zptb, zptbtl )
223      !
224      DO jn = 1, kjpt                                     ! tracer loop
225         !                                                ! ===========
226#  if defined key_vectopt_loop
227         DO jj = 1, 1   ! vector opt. (forced unrolling)
228            DO ji = 1, jpij
229#else
230         DO jj = 1, jpj
231            DO ji = 1, jpi
232#endif
233               ik = mbkt(ji,jj)                        ! bottom T-level index
234               zptbtl(ji,jj) = ptb_tl(ji,jj,ik,jn)              ! bottom before T and S
235               zptb  (ji,jj) = ptb   (ji,jj,ik,jn)     ! bottom before T and S
236            END DO
237         END DO
238         !                                                ! Compute the trend
239#  if defined key_vectopt_loop
240         DO jj = 1, 1   ! vector opt. (forced unrolling)
241            DO ji = jpi+1, jpij-jpi-1
242#  else
243         DO jj = 2, jpjm1
244            DO ji = 2, jpim1
245#  endif
246               ik = mbkt(ji,jj)                            ! bottom T-level index
247               zbtr = r1_e1e2t(ji,jj)  / fse3t(ji,jj,ik)
248               pta_tl(ji,jj,ik,jn) = pta_tl(ji,jj,ik,jn)                                                         &
249# if defined control_param
250                  &               + (   ahu_bbl_tl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )    &
251                  &                   - ahu_bbl_tl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )    &
252                  &                   + ahv_bbl_tl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )    &
253                  &                   - ahv_bbl_tl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr &
254#endif
255                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptbtl(ji+1,jj  ) - zptbtl(ji  ,jj  ) )   &
256                  &                   - ahu_bbl(ji-1,jj  ) * ( zptbtl(ji  ,jj  ) - zptbtl(ji-1,jj  ) )   &
257                  &                   + ahv_bbl(ji  ,jj  ) * ( zptbtl(ji  ,jj+1) - zptbtl(ji  ,jj  ) )   &
258                  &                   - ahv_bbl(ji  ,jj-1) * ( zptbtl(ji  ,jj  ) - zptbtl(ji  ,jj-1) )   ) * zbtr
259            END DO
260         END DO
261         !                                                  ! ===========
262      END DO                                                ! end tracer
263      !                                                     ! ===========
264      CALL wrk_dealloc( jpi, jpj, zptbtl, zptb )
265      !
266      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif_tan')
267      !
268   END SUBROUTINE tra_bbl_dif_tan
269
270
271   SUBROUTINE tra_bbl_adv_tan( ptb, ptb_tl, pta_tl, kjpt )
272      !!----------------------------------------------------------------------
273      !!                  ***  ROUTINE trc_bbl  ***
274      !!
275      !! ** Purpose :   Compute the before passive tracer trend associated
276      !!     with the bottom boundary layer and add it to the general trend
277      !!     of tracer equations.
278      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
279      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
280      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
281      !!                       transport proportional to the along-slope density gradient
282      !!
283      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
284      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
285      !!----------------------------------------------------------------------
286      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
287      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
288      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb_tl ! before and now tangent tracer fields
289      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_tl    ! tracer trend
290      !
291      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
292      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
293      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
294      REAL(wp) ::   zbtr, ztra               ! local scalars
295      REAL(wp) ::   ztratl                   !   -      -
296      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
297      REAL(wp) ::   zu_bbltl, zv_bbltl       !   -      -
298      !!----------------------------------------------------------------------
299      !
300      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv_tan')
301      !                                                          ! ===========
302      DO jn = 1, kjpt                                            ! tracer loop
303         !                                                       ! ===========
304# if defined key_vectopt_loop
305         DO jj = 1, 1
306            DO ji = 1, jpij-jpi-1   ! vector opt. (forced unrolling)
307# else
308         DO jj = 1, jpjm1
309            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
310# endif
311               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
312                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
313                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
314                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
315                  zu_bbl = ABS( utr_bbl(ji,jj) )
316                  zu_bbltl = SIGN( utr_bbl_tl(ji,jj), utr_bbl(ji,jj) )
317                  !
318                  !                                               ! up  -slope T-point (shelf bottom point)
319                  zbtr = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus)
320                  ztratl = ( zu_bbltl * ( ptb   (iid,jj,ikus,jn) - ptb   (iis,jj,ikus,jn) )   &
321                     &     + zu_bbl   * ( ptb_tl(iid,jj,ikus,jn) - ptb_tl(iis,jj,ikus,jn) ) ) * zbtr
322                  pta_tl(iis,jj,ikus,jn) = pta_tl(iis,jj,ikus,jn) + ztratl
323                  !
324                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
325                     zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk)
326                     ztratl = ( zu_bbltl * ( ptb   (iid,jj,jk+1,jn) - ptb   (iid,jj,jk,jn) ) &
327                        &     + zu_bbl   * ( ptb_tl(iid,jj,jk+1,jn) - ptb_tl(iid,jj,jk,jn) ) ) * zbtr
328                     pta_tl(iid,jj,jk,jn) = pta_tl(iid,jj,jk,jn) + ztratl
329                  END DO
330                  !
331                  zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud)
332                  ztratl = ( zu_bbltl * ( ptb   (iis,jj,ikus,jn) - ptb   (iid,jj,ikud,jn) ) &
333                     &     + zu_bbl   * ( ptb_tl(iis,jj,ikus,jn) - ptb_tl(iid,jj,ikud,jn) ) ) * zbtr
334                  pta_tl(iid,jj,ikud,jn) = pta_tl(iid,jj,ikud,jn) + ztratl
335               ENDIF
336               !
337               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
338                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
339                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
340                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
341                  zv_bbl = ABS( vtr_bbl(ji,jj) )
342                  zv_bbltl = SIGN( vtr_bbl_tl(ji,jj), vtr_bbl(ji,jj) )
343                  !
344                  ! up  -slope T-point (shelf bottom point)
345                  zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs)
346                  ztratl = ( zv_bbltl * ( ptb   (ji,ijd,ikvs,jn) - ptb   (ji,ijs,ikvs,jn) ) &
347                     &     + zv_bbl   * ( ptb_tl(ji,ijd,ikvs,jn) - ptb_tl(ji,ijs,ikvs,jn) ) ) * zbtr
348                  pta_tl(ji,ijs,ikvs,jn) = pta_tl(ji,ijs,ikvs,jn) + ztratl
349                  !
350                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
351                     zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk)
352                     ztratl = ( zv_bbltl * ( ptb   (ji,ijd,jk+1,jn) - ptb   (ji,ijd,jk,jn) ) &
353                        &     + zv_bbl   * ( ptb_tl(ji,ijd,jk+1,jn) - ptb_tl(ji,ijd,jk,jn) ) ) * zbtr
354                     pta_tl(ji,ijd,jk,jn) = pta_tl(ji,ijd,jk,jn)  + ztratl
355                  END DO
356                  !                                               ! down-slope T-point (deep bottom point)
357                  zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd)
358                  ztratl = ( zv_bbltl * ( ptb   (ji,ijs,ikvs,jn) - ptb   (ji,ijd,ikvd,jn) ) &
359                     &     + zv_bbl   * ( ptb_tl(ji,ijs,ikvs,jn) - ptb_tl(ji,ijd,ikvd,jn) ) ) * zbtr
360                  pta_tl(ji,ijd,ikvd,jn) = pta_tl(ji,ijd,ikvd,jn) + ztratl
361               ENDIF
362            END DO
363            !
364         END DO
365         !                                                  ! ===========
366      END DO                                                ! end tracer
367      !                                                     ! ===========
368      !
369      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv_tan')
370      !
371   END SUBROUTINE tra_bbl_adv_tan
372
373
374   SUBROUTINE bbl_tan( kt, kit000, cdtype )
375      !!----------------------------------------------------------------------
376      !!                  ***  ROUTINE bbl  ***
377      !!
378      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
379      !!                advection terms.
380      !!
381      !! ** Method  :
382      !!        * diffusive bbl (nn_bbl_ldf=1) :
383      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
384      !!      along bottom slope gradient) an additional lateral 2nd order
385      !!      diffusion along the bottom slope is added to the general
386      !!      tracer trend, otherwise the additional trend is set to 0.
387      !!      A typical value of ahbt is 2000 m2/s (equivalent to
388      !!      a downslope velocity of 20 cm/s if the condition for slope
389      !!      convection is satified)
390      !!        * advective bbl (nn_bbl_adv=1 or 2) :
391      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
392      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
393      !!        i.e. transport proportional to the along-slope density gradient
394      !!
395      !!      NB: the along slope density gradient is evaluated using the
396      !!      local density (i.e. referenced at a common local depth).
397      !!
398      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
399      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
400      !!----------------------------------------------------------------------
401      !
402      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
403      INTEGER         , INTENT(in   ) ::   kit000          ! first time step index
404      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
405      !!
406      INTEGER  ::   ji, jj                    ! dummy loop indices
407      INTEGER  ::   ik                        ! local integers
408      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
409      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
410      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
411      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
412      REAL(wp) ::   zgdrhotl, zttl, zstl, zhtl!   -      -
413      !!
414      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
415      REAL(wp) ::   fsalbt_tan, fsbeta_tan, pfttl, pfstl, pfhtl   ! statement function
416      REAL(wp), POINTER, DIMENSION(:,:) :: zub  , zvb  , ztb  , zsb  , zdep
417      REAL(wp), POINTER, DIMENSION(:,:) :: zubtl, zvbtl, ztbtl, zsbtl
418      !!----------------------- zv_bbl-----------------------------------------------
419      ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
420      ! ================            pft :  potential temperature in degrees celcius
421      !                             pfs :  salinity anomaly (s-35) in psu
422      !                             pfh :  depth in meters
423      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
424      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
425         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
426                                   - 0.203814e-03 ) * pft                    &
427                                   + 0.170907e-01 ) * pft                    &
428                                   + 0.665157e-01                            &
429         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
430         +  ( ( - 0.302285e-13 * pfh                                         &
431                - 0.251520e-11 * pfs                                         &
432                + 0.512857e-12 * pft * pft          ) * pfh                  &
433                                     - 0.164759e-06   * pfs                  &
434             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
435                                     + 0.380374e-04 ) * pfh
436      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
437         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
438                                 - 0.301985e-05 ) * pft                      &
439                                 + 0.785567e-03                              &
440         + (     0.515032e-08 * pfs                                          &
441               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
442               +(  (   0.121551e-17 * pfh                                    &
443                     - 0.602281e-15 * pfs                                    &
444                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
445                                          + 0.408195e-10   * pfs             &
446                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
447                                          - 0.121555e-07 ) * pfh
448
449      fsalbt_tan( pft, pfs, pfh, pfttl, pfstl, pfhtl ) =  &   ! alpha/beta
450         &         ( - 0.255019e-07 * 4 * pft * pft * pft            &
451         &           + 0.298357e-05 * 3 * pft * pft                  &
452         &           - 0.203814e-03 * 2 * pft                        &
453         &           - 0.846960e-04 * pfs                            &
454         &           + 0.512857e-12 * 2 * pft * pfh * pfh            &
455         &           + 0.791325e-08 * pft * pfh                      &
456         &           - 0.933746e-06 * pfh                            &
457         &           + 0.170907e-01                        ) * pfttl &
458         &       + ( - 0.678662e-05 * 2 * pfs                        &
459         &           - 0.846960e-04 * pft                            &
460         &           - 0.251520e-11 * pfh  * pfh                     &
461         &           - 0.164759e-06 * pfh                            &
462         &           + 0.378110e-02                        ) * pfstl &
463         &       + ( - 0.302285e-13 * 3 * pfh  * pfh                 &
464         &           - 0.251520e-11 * pfs * pfh                      &
465         &           + 0.512857e-12 * pft * pft * pfh                &
466         &           - 0.164759e-06 * pfs                            &
467         &           + 0.791325e-08 * pft * pft                      &
468         &           - 0.933746e-06 * pft                            &
469         &           + 0.380374e-04                        ) * pfhtl
470
471
472      fsbeta_tan( pft, pfs, pfh, pfttl, pfstl, pfhtl ) =  &   ! beta
473         &         ( - 0.415613e-09 * 3 * pft * pft                  &
474         &           + 0.555579e-07 * 2 * pft                        &
475         &           - 0.301985e-05                                  &
476         &           + 0.788212e-08 * pfs                            &
477         &           - 0.213127e-11 * 2 * pfh * pft                  &
478         &           - 0.175379e-14 * pfh * pfh            ) * pfttl &
479         &       + (   0.788212e-08 * pft                            &
480         &           + 0.515032e-08 * 2 * pfs                        &
481         &           - 0.356603e-06                                  &
482         &           + 0.408195e-10 * pfh                            &
483         &           - 0.602281e-15 * pfh * pfh            ) * pfstl &
484         &       + (   0.121551e-17 * 3 * pfh * pfh                  &
485         &           - 0.602281e-15 * 2 * pfs * pfh                  &
486         &           - 0.175379e-14 * 2 * pft * pfh                  &
487         &           + 0.176621e-12 * 2 * pfh                        &
488         &           + 0.408195e-10 * pfs                            &
489         &           + 0.192867e-09 * pfh                            &
490         &           - 0.213127e-11 * pft * pft                      &
491         &           + 0.192867e-09 * pft                            &
492         &           - 0.121555e-07                        ) * pfhtl
493
494
495      !!----------------------------------------------------------------------
496
497      !
498      IF( nn_timing == 1 )  CALL timing_start( 'bbl_tan')
499      !
500      CALL wrk_alloc( jpi, jpj, zub  , zvb  , ztb  , zsb  , zdep, &
501         &                      zubtl, zvbtl, ztbtl, zsbtl        )
502      !
503
504      IF( kt == kit000 )  THEN
505         IF(lwp)  WRITE(numout,*)
506         IF(lwp)  WRITE(numout,*) 'trabbl_tam:bbl_tan : Compute bbl velocities and diffusive coefficients in ', cdtype
507         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
508      ENDIF
509
510      !                                        !* bottom temperature, salinity, velocity and depth
511#if defined key_vectopt_loop
512      DO jj = 1, 1   ! vector opt. (forced unrolling)
513         DO ji = 1, jpij
514#else
515      DO jj = 1, jpj
516         DO ji = 1, jpi
517#endif
518            ik = mbkt(ji,jj)                        ! bottom T-level index
519            ztb  (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
520            zsb  (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
521            ztbtl(ji,jj) = tsb_tl(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
522            zsbtl(ji,jj) = tsb_tl(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
523            zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
524            !
525            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
526            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
527         END DO
528      END DO
529
530      !                                   !-------------------!
531      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
532         !                                !-------------------!
533         ! AV NOTE : while rn_ahtbbl remains a passive variable, the code below will only yield ah_bbl_tl=0, so i put it under key
534#if defined key_control_param
535         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
536            DO ji = 1, jpim1
537               !                                                ! i-direction
538               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
539               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
540               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
541               !                                                         ! masked bbl i-gradient of density
542               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
543                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
544               !
545               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
546               ahu_bbl_tl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0_tl(ji,jj)             ! masked diffusive flux coeff.
547               !
548               !                                                ! j-direction
549               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
550               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
551               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
552               !                                                         ! masked bbl j-gradient of density
553               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
554                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
555               !
556               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
557               ahv_bbl_tl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0_tl(ji,jj)
558               !
559            END DO
560         END DO
561#else
562         DO jj = 1, jpjm1
563            DO ji = 1, jpim1
564               ahu_bbl_tl(ji,jj)=0.0_wp
565               ahv_bbl_tl(ji,jj)=0.0_wp
566            END DO
567         END DO
568#endif
569         !
570      ENDIF
571
572      !                                   !-------------------!
573      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
574         !                                !-------------------!
575         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
576         !
577         CASE( 1 )                                   != use of upper velocity
578            ! AV NOTE: not much needed for deriving, almost all the computations are for the SIGN, which is kept identical as in the NL
579            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
580               DO ji = 1, fs_jpim1   ! vector opt.
581                  !                                               ! i-direction
582                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
583                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
584                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
585                  !                                                           ! masked bbl i-gradient of density
586                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
587                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
588                  !
589                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
590                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
591                  !
592                  !                                                           ! bbl velocity
593                  utr_bbl_tl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zubtl(ji,jj)
594                  !
595                  !                                               ! j-direction
596                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
597                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
598                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
599                  !                                                           ! masked bbl j-gradient of density
600                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
601                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
602                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
603                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
604                  !
605                  !                                                           ! bbl velocity
606                  vtr_bbl_tl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvbtl(ji,jj)
607               END DO
608            END DO
609            !
610         CASE( 2 )                                 != bbl velocity = F( delta rho )
611            ! AV NOTE: this one is nastier
612            zgbbl = grav * rn_gambbl
613            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
614               DO ji = 1, fs_jpim1   ! vector opt.
615                  !                                         ! i-direction
616                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
617                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
618                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
619                  !
620                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
621                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
622                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
623                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
624                  zgdrho =    fsbeta( zt, zs, zh )                                    &
625                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    &
626                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
627                  zttl = 0.5 * ( ztbtl (ji,jj) + ztbtl (ji+1,jj) )           ! mid slope depth of T, S, and depth
628                  zstl = 0.5 * ( zsbtl (ji,jj) + zsbtl (ji+1,jj) )
629                  zhtl = 0.0_wp
630                  zgdrhotl =  ( fsbeta_tan( zt, zs, zh, zttl, zstl, zhtl )                     &
631                     &        * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )         &
632                     &                                 - ( zsb(iid,jj) - zsb(iis,jj) ) )       &
633                     &        + fsbeta( zt, zs, zh )                                           &
634                     &        * ( fsalbt_tan( zt, zs, zh, zttl, zstl, zhtl )                   &
635                     &                                     * ( ztb  (iid,jj) - ztb  (iis,jj) ) &
636                     &          + fsalbt    ( zt, zs, zh ) * ( ztbtl(iid,jj) - ztbtl(iis,jj) ) &
637                     &                                  - ( zsbtl(iid,jj) - zsbtl(iis,jj) )  ) ) * umask(ji,jj,1)
638
639                  zsign  = SIGN( 0.5_wp, zgdrho ) ! tangent of zgdrho = MAX( 0.e0, zgdrho )
640                  !                                             ! bbl transport (down-slope direction)
641                  utr_bbl_tl(ji,jj) = zsign * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrhotl * REAL( mgrhu(ji,jj) )
642                  !
643                  !                                         ! j-direction
644                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
645                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
646                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
647                  !
648                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
649                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
650                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
651                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
652                  zgdrho =    fsbeta( zt, zs, zh )                                    &
653                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    &
654                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
655                  zttl = 0.5 * ( ztbtl (ji,jj) + ztbtl (ji,jj+1) )           ! mid slope depth of T, S, and depth
656                  zstl = 0.5 * ( zsbtl (ji,jj) + zsbtl (ji,jj+1) )
657                  zhtl = 0.0_wp
658                  zgdrhotl =  ( fsbeta_tan( zt, zs, zh, zttl, zstl, zhtl )                     &
659                     &        * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )         &
660                     &                                 - ( zsb(ji,ijd) - zsb(ji,ijs) ) )       &
661                     &        + fsbeta( zt, zs, zh )                                           &
662                     &        * ( fsalbt_tan( zt, zs, zh, zttl, zstl, zhtl )                   &
663                     &                                     * ( ztb  (ji,ijd) - ztb  (ji,ijs) ) &
664                     &          + fsalbt    ( zt, zs, zh ) * ( ztbtl(ji,ijd) - ztbtl(ji,ijs) ) &
665                     &                                     - ( zsbtl(ji,ijd) - zsbtl(ji,ijs) ) ) ) * vmask(ji,jj,1)
666                  !
667                  zsign  = SIGN( 0.5_wp, zgdrho ) ! tangent of zgdrho = MAX( 0.e0, zgdrho )
668                  !                                             ! bbl transport (down-slope direction)
669                  vtr_bbl_tl(ji,jj) = zsign * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrhotl * REAL( mgrhv(ji,jj) )
670               END DO
671            END DO
672         END SELECT
673         !
674      ENDIF
675      !
676      CALL wrk_dealloc( jpi, jpj, zub  , zvb  , ztb  , zsb  , zdep,   &
677         &                        zubtl, zvbtl, ztbtl, zsbtl        )
678      !
679      IF( nn_timing == 1 )  CALL timing_stop( 'bbl_tan')
680      !
681   END SUBROUTINE bbl_tan
682
683
684   SUBROUTINE tra_bbl_adj( kt )
685      !!----------------------------------------------------------------------
686      !!                  ***  ROUTINE bbl_adj  ***
687      !!
688      !! ** Purpose :   Compute the before tracer (t & s) trend associated
689      !!              with the bottom boundary layer and add it to the general
690      !!              trend of tracer equations.
691      !!
692      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
693      !!              diffusive and/or advective contribution to the tracer trend
694      !!              is added to the general tracer trend
695      !!----------------------------------------------------------------------
696      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
697      !!----------------------------------------------------------------------
698      !
699      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adj')
700      !
701      IF( l_bbl )  CALL bbl_adj( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
702
703      IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl
704         !
705         CALL tra_bbl_dif_adj( tsb, tsb_ad, tsa_ad, jpts )
706         !
707      END IF
708
709      IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl
710         !
711         CALL tra_bbl_adv_adj( tsb, tsb_ad, tsa_ad, jpts )
712         !
713      END IF
714      !
715      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adj')
716      !
717   END SUBROUTINE tra_bbl_adj
718
719
720   SUBROUTINE tra_bbl_dif_adj( ptb, ptb_ad, pta_ad, kjpt )
721      !!----------------------------------------------------------------------
722      !!                  ***  ROUTINE tra_bbl_dif_adj  ***
723      !!
724      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
725      !!                advection terms.
726      !!
727      !! ** Method  :
728      !!        * diffusive bbl (nn_bbl_ldf=1) :
729      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
730      !!      along bottom slope gradient) an additional lateral 2nd order
731      !!      diffusion along the bottom slope is added to the general
732      !!      tracer trend, otherwise the additional trend is set to 0.
733      !!      A typical value of ahbt is 2000 m2/s (equivalent to
734      !!      a downslope velocity of 20 cm/s if the condition for slope
735      !!      convection is satified)
736      !!
737      !! ** Action  :   pta   increased by the bbl diffusive trend
738      !!
739      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
740      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
741      !!----------------------------------------------------------------------
742      !
743      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
744      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
745      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptb_ad ! before and now tracer fields
746      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_ad    ! tracer trend
747      !
748      INTEGER  ::   ji, jj, jn   ! dummy loop indices
749      INTEGER  ::   ik           ! local integers
750      REAL(wp) ::   zbtr         ! local scalars
751      REAL(wp), POINTER, DIMENSION(:,:) :: zptb, zptbad
752      !!----------------------------------------------------------------------
753      !
754      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif_adj')
755      !
756      CALL wrk_alloc( jpi, jpj, zptb, zptbad )
757      zptbad(:,:) = 0.0_wp
758      !
759      DO jn = 1, kjpt                                     ! tracer loop
760         !                                                ! ===========
761#  if defined key_vectopt_loop
762         DO jj = 1, 1   ! vector opt. (forced unrolling)
763            DO ji = 1, jpij
764#  else
765         DO jj = 1, jpj
766            DO ji = 1, jpi
767#endif
768               ik = mbkt(ji,jj)                        ! bottom T-level index
769               zptb  (ji,jj) = ptb   (ji,jj,ik,jn)     ! bottom before T and S
770            END DO
771         END DO
772         !                                                  ! ===========
773         !                                             ! Compute the trend
774#  if defined key_vectopt_loop
775         DO jj = 1, 1   ! vector opt. (forced unrolling)
776            DO ji = jpi+1, jpij-jpi-1
777#  else
778         DO jj = jpjm1, 2, -1
779            DO ji = jpim1, 2, -1
780#  endif
781               ik = mbkt(ji,jj)                            ! bottom T-level index
782               zbtr = r1_e1e2t(ji,jj)  / fse3t(ji,jj,ik)
783#  if defined control_param
784               ahu_bbl_ad(ji  ,jj  ) = ahu_bbl_ad(ji  ,jj  ) + pta_ad(ji,jj,ik,jn) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) ) * zbtr
785               ahu_bbl_ad(ji-1,jj  ) = ahu_bbl_ad(ji-1,jj  ) - pta_ad(ji,jj,ik,jn) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) ) * zbtr
786               ahv_bbl_ad(ji  ,jj  ) = ahv_bbl_ad(ji  ,jj  ) + pta_ad(ji,jj,ik,jn) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) ) * zbtr
787               ahv_bbl_ad(ji  ,jj-1) = ahv_bbl_ad(ji  ,jj-1) - pta_ad(ji,jj,ik,jn) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) ) * zbtr
788#  endif
789               zptbad(ji  ,jj  ) = zptbad(ji  ,jj  ) - pta_ad(ji,jj,ik,jn) * ( ahu_bbl(ji  ,jj  ) + ahu_bbl(ji-1,jj  ) &
790                  &                                                          + ahv_bbl(ji  ,jj  ) + ahv_bbl(ji  ,jj-1) ) * zbtr
791               zptbad(ji+1,jj  ) = zptbad(ji+1,jj  ) + pta_ad(ji,jj,ik,jn) * ahu_bbl(ji  ,jj  ) * zbtr
792               zptbad(ji-1,jj  ) = zptbad(ji-1,jj  ) + pta_ad(ji,jj,ik,jn) * ahu_bbl(ji-1,jj  ) * zbtr
793               zptbad(ji  ,jj+1) = zptbad(ji  ,jj+1) + pta_ad(ji,jj,ik,jn) * ahv_bbl(ji  ,jj  ) * zbtr
794               zptbad(ji  ,jj-1) = zptbad(ji  ,jj-1) + pta_ad(ji,jj,ik,jn) * ahv_bbl(ji  ,jj-1) * zbtr
795
796               pta_ad(ji,jj,ik,jn) = pta_ad(ji,jj,ik,jn)                                                       &
797                  &                 + ( ahu_bbl(ji  ,jj  ) * ( zptbad(ji+1,jj  ) - zptbad(ji  ,jj  ) )   &
798                  &                   - ahu_bbl(ji-1,jj  ) * ( zptbad(ji  ,jj  ) - zptbad(ji-1,jj  ) )   &
799                  &                   + ahv_bbl(ji  ,jj  ) * ( zptbad(ji  ,jj+1) - zptbad(ji  ,jj  ) )   &
800                  &                   - ahv_bbl(ji  ,jj-1) * ( zptbad(ji  ,jj  ) - zptbad(ji  ,jj-1) )   ) * zbtr
801            END DO
802         END DO
803#  if defined key_vectopt_loop
804         DO jj = 1, 1   ! vector opt. (forced unrolling)
805            DO ji = 1, jpij
806#else
807         DO jj = 1, jpj
808            DO ji = 1, jpi
809#endif
810               ik = mbkt(ji,jj)                        ! bottom T-level index
811               ptb_ad(ji,jj,ik,jn) = ptb_ad(ji,jj,ik,jn) + zptbad(ji,jj)
812               zptbad(ji,jj) = 0.0_wp                  ! bottom before T and S
813            END DO
814         END DO
815         !                                                  ! ===========
816      END DO                                                ! end tracer
817      !                                                     ! ===========
818      CALL wrk_dealloc( jpi, jpj, zptbad, zptb )
819      !
820      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif_adj')
821      !
822   END SUBROUTINE tra_bbl_dif_adj
823
824
825   SUBROUTINE tra_bbl_adv_adj( ptb, ptb_ad, pta_ad, kjpt )
826      !!----------------------------------------------------------------------
827      !!                  ***  ROUTINE trc_bbl  ***
828      !!
829      !! ** Purpose :   Compute the before passive tracer trend associated
830      !!     with the bottom boundary layer and add it to the general trend
831      !!     of tracer equations.
832      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
833      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
834      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
835      !!                       transport proportional to the along-slope density gradient
836      !!
837      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
838      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
839      !!----------------------------------------------------------------------
840      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
841      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
842      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   ptb_ad ! before and now adjoint tracer fields
843      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta_ad    ! tracer trend
844      !
845      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
846      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
847      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
848      REAL(wp) ::   zbtr, ztra               ! local scalars
849      REAL(wp) ::   ztraad                   !   -      -
850      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
851      REAL(wp) ::   zu_bblad, zv_bblad       !   -      -
852      !!----------------------------------------------------------------------
853      !
854      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv_adj')
855      !
856      zu_bblad = 0.0_wp ; zv_bblad = 0.0_wp
857      !                                                          ! ===========
858      DO jn = 1, kjpt                                            ! tracer loop
859         !                                                       ! ===========
860# if defined key_vectopt_loop
861         DO jj = 1, 1
862            DO ji = jpij-jpi-1, 1, -1   ! vector opt. (forced unrolling)
863# else
864         DO jj = jpjm1, 1, -1
865            DO ji = jpim1, 1, -1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
866# endif
867               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
868                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
869                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
870                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
871                  zv_bbl   = ABS ( vtr_bbl(ji,jj) )
872                  !                                               ! down-slope T-point (deep bottom point)
873                  zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd)
874                  ztraad = pta_ad(ji,ijd,ikvd,jn)
875                  zv_bblad = zv_bblad + ztraad * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
876                  ptb_ad(ji,ijs,ikvs,jn) = ptb_ad(ji,ijs,ikvs,jn) + ztraad * zv_bbl * zbtr
877                  ptb_ad(ji,ijd,ikvd,jn) = ptb_ad(ji,ijd,ikvd,jn) - ztraad * zv_bbl * zbtr
878                  !
879                  DO jk = ikvd-1, ikvs, -1                        ! down-slope upper to down T-point (deep column)
880                     zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk)
881                     ztraad = pta_ad(ji,ijd,jk,jn)
882                     zv_bblad = zv_bblad + ztraad * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
883                     ptb_ad(ji,ijd,jk+1,jn) = ptb_ad(ji,ijd,jk+1,jn) + ztraad * zv_bbl * zbtr
884                     ptb_ad(ji,ijd,jk  ,jn) = ptb_ad(ji,ijd,jk  ,jn) - ztraad * zv_bbl * zbtr
885                  END DO
886                  ! up  -slope T-point (shelf bottom point)
887                  zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs)
888                  ztraad = pta_ad(ji,ijs,ikvs,jn)
889                  zv_bblad = zv_bblad + ztraad * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
890                  ptb_ad(ji,ijd,ikvs,jn) = ptb_ad(ji,ijd,ikvs,jn) + ztraad * zv_bbl * zbtr
891                  ptb_ad(ji,ijs,ikvs,jn) = ptb_ad(ji,ijs,ikvs,jn) - ztraad * zv_bbl * zbtr
892
893                  !
894                  vtr_bbl_ad(ji,jj) = vtr_bbl_ad(ji,jj) + SIGN( zv_bblad, vtr_bbl(ji,jj) )
895                  zv_bblad = 0.0_wp
896                  !
897               ENDIF
898
899
900
901               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
902                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
903                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
904                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
905                  zu_bbl = ABS( utr_bbl(ji,jj) )
906                  !
907                  zbtr   = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud)
908                  ztraad = pta_ad(iid,jj,ikud,jn)
909                  zu_bblad = zu_bblad + ztraad * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
910                  ptb_ad(iis,jj,ikus,jn) = ptb_ad(iis,jj,ikus,jn) + ztraad * zu_bbl * zbtr
911                  ptb_ad(iid,jj,ikud,jn) = ptb_ad(iid,jj,ikud,jn) + ztraad * zu_bbl * zbtr
912                  !
913                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
914                     zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk)
915                     ztraad = pta_ad(iid,jj,jk,jn)
916                     zu_bblad = zu_bblad + ztraad * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
917                     ptb_ad(iid,jj,jk+1,jn) = ptb_ad(iid,jj,jk+1,jn) + ztraad * zu_bbl * zbtr
918                     ptb_ad(iid,jj,jk  ,jn) = ptb_ad(iid,jj,jk  ,jn) - ztraad * zu_bbl * zbtr
919                  END DO
920                  !                                               ! up  -slope T-point (shelf bottom point)
921                  zbtr   = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus)
922                  ztraad = pta_ad(iis,jj,ikus,jn)
923                  zu_bblad = zu_bblad + ztraad * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
924                  ptb_ad(iid,jj,ikus,jn) = ptb_ad(iid,jj,ikus,jn) + ztraad * zu_bbl * zbtr
925                  ptb_ad(iis,jj,ikus,jn) = ptb_ad(iis,jj,ikus,jn) - ztraad * zu_bbl * zbtr
926                  !
927                  utr_bbl_ad(ji,jj) = utr_bbl_ad(ji,jj) + SIGN( zu_bblad, utr_bbl(ji,jj) )
928                  zu_bblad = 0.0_wp
929                  !
930               ENDIF
931               !
932            END DO
933            !
934         END DO
935         !                                                  ! ===========
936      END DO                                                ! end tracer
937      !                                                     ! ===========
938      !
939      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv_adj')
940      !
941   END SUBROUTINE tra_bbl_adv_adj
942
943
944   SUBROUTINE bbl_adj( kt, kit000, cdtype )
945      !!----------------------------------------------------------------------
946      !!                  ***  ROUTINE bbl  ***
947      !!
948      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
949      !!                advection terms.
950      !!
951      !! ** Method  :
952      !!        * diffusive bbl (nn_bbl_ldf=1) :
953      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
954      !!      along bottom slope gradient) an additional lateral 2nd order
955      !!      diffusion along the bottom slope is added to the general
956      !!      tracer trend, otherwise the additional trend is set to 0.
957      !!      A typical value of ahbt is 2000 m2/s (equivalent to
958      !!      a downslope velocity of 20 cm/s if the condition for slope
959      !!      convection is satified)
960      !!        * advective bbl (nn_bbl_adv=1 or 2) :
961      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
962      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
963      !!        i.e. transport proportional to the along-slope density gradient
964      !!
965      !!      NB: the along slope density gradient is evaluated using the
966      !!      local density (i.e. referenced at a common local depth).
967      !!
968      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
969      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
970      !!----------------------------------------------------------------------
971      !
972      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
973      INTEGER         , INTENT(in   ) ::   kit000          ! first time step index
974      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
975      !!
976      INTEGER  ::   ji, jj                    ! dummy loop indices
977      INTEGER  ::   ik                        ! local integers
978      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
979      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
980      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
981      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
982      REAL(wp) ::   zgdrhoad, ztad, zsad, zhad!   -      -
983      !!
984      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
985      REAL(wp) ::   pftad, pfsad, pfhad
986      REAL(wp) ::   fsalbt_adj_t, fsbeta_adj_t
987      REAL(wp) ::   fsalbt_adj_s, fsbeta_adj_s
988      REAL(wp) ::   fsalbt_adj_h, fsbeta_adj_h
989      REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep
990      REAL(wp), POINTER, DIMENSION(:,:) :: zubad, zvbad, ztbad, zsbad
991      !!----------------------- zv_bbl-----------------------------------------------
992      ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
993      ! ================            pft :  potential temperature in degrees celcius
994      !                             pfs :  salinity anomaly (s-35) in psu
995      !                             pfh :  depth in meters
996      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
997      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
998         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
999                                   - 0.203814e-03 ) * pft                    &
1000                                   + 0.170907e-01 ) * pft                    &
1001                                   + 0.665157e-01                            &
1002         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
1003         +  ( ( - 0.302285e-13 * pfh                                         &
1004                - 0.251520e-11 * pfs                                         &
1005                + 0.512857e-12 * pft * pft          ) * pfh                  &
1006                                     - 0.164759e-06   * pfs                  &
1007             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
1008                                     + 0.380374e-04 ) * pfh
1009      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
1010         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
1011                                 - 0.301985e-05 ) * pft                      &
1012                                 + 0.785567e-03                              &
1013         + (     0.515032e-08 * pfs                                          &
1014               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
1015               +(  (   0.121551e-17 * pfh                                    &
1016                     - 0.602281e-15 * pfs                                    &
1017                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
1018                                          + 0.408195e-10   * pfs             &
1019                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
1020                                          - 0.121555e-07 ) * pfh
1021
1022      fsalbt_adj_t( pft, pfs, pfh, pftad ) =  &   ! alpha/beta
1023         &         ( - 0.255019e-07 * 4 * pft * pft * pft            &
1024         &           + 0.298357e-05 * 3 * pft * pft                  &
1025         &           - 0.203814e-03 * 2 * pft                        &
1026         &           - 0.846960e-04 * pfs                            &
1027         &           + 0.512857e-12 * 2 * pft * pfh * pfh            &
1028         &           + 0.791325e-08 * pft * pfh                      &
1029         &           - 0.933746e-06 * pfh                            &
1030         &           + 0.170907e-01                        ) * pftad
1031
1032      fsalbt_adj_s( pft, pfs, pfh, pfsad ) =  &   ! alpha/beta
1033         &       + ( - 0.678662e-05 * 2 * pfs                        &
1034         &           - 0.846960e-04 * pft                            &
1035         &           - 0.251520e-11 * pfh  * pfh                     &
1036         &           - 0.164759e-06 * pfh                            &
1037         &           + 0.378110e-02                        ) * pfsad
1038
1039      fsalbt_adj_h( pft, pfs, pfh, pfhad ) =  &   ! alpha/beta
1040         &       + ( - 0.302285e-13 * 3 * pfh  * pfh                 &
1041         &           - 0.251520e-11 * pfs * pfh                      &
1042         &           + 0.512857e-12 * pft * pft * pfh                &
1043         &           - 0.164759e-06 * pfs                            &
1044         &           + 0.791325e-08 * pft * pft                      &
1045         &           - 0.933746e-06 * pft                            &
1046         &           + 0.380374e-04                        ) * pfhad
1047
1048
1049      fsbeta_adj_t( pft, pfs, pfh, pftad ) =  &   ! beta
1050         &         ( - 0.415613e-09 * 3 * pft * pft                  &
1051         &           + 0.555579e-07 * 2 * pft                        &
1052         &           - 0.301985e-05                                  &
1053         &           + 0.788212e-08 * pfs                            &
1054         &           - 0.213127e-11 * 2 * pfh * pft                  &
1055         &           - 0.175379e-14 * pfh * pfh            ) * pftad
1056      fsbeta_adj_s( pft, pfs, pfh, pfsad ) =  &   ! beta
1057         &         (   0.788212e-08 * pft                            &
1058         &           + 0.515032e-08 * 2 * pfs                        &
1059         &           - 0.356603e-06                                  &
1060         &           + 0.408195e-10 * pfh                            &
1061         &           - 0.602281e-15 * pfh * pfh            ) * pfsad
1062      fsbeta_adj_h( pft, pfs, pfh, pfhad ) =  &   ! beta
1063         &         (   0.121551e-17 * 3 * pfh * pfh                  &
1064         &           - 0.602281e-15 * 2 * pfs * pfh                  &
1065         &           - 0.175379e-14 * 2 * pft * pfh                  &
1066         &           + 0.176621e-12 * 2 * pfh                        &
1067         &           + 0.408195e-10 * pfs                            &
1068         &           + 0.192867e-09 * pfh                            &
1069         &           - 0.213127e-11 * pft * pft                      &
1070         &           + 0.192867e-09 * pft                            &
1071         &           - 0.121555e-07                        ) * pfhad
1072      !!----------------------------------------------------------------------
1073
1074      !
1075      IF( nn_timing == 1 )  CALL timing_start( 'bbl_adj')
1076      !
1077      CALL wrk_alloc( jpi, jpj, zub  , zvb  , ztb  , zsb, zdep, &
1078         &                      zubad, zvbad, ztbad, zsbad      )
1079      !
1080      zubad(:,:) = 0.0_wp ; zvbad(:,:) = 0.0_wp ; ztbad(:,:) = 0.0_wp ; zsbad(:,:) = 0.0_wp
1081
1082      IF( kt == kit000 )  THEN
1083         IF(lwp)  WRITE(numout,*)
1084         IF(lwp)  WRITE(numout,*) 'trabbl_tam:bbl_adj : Compute bbl velocities and diffusive coefficients in ', cdtype
1085         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
1086      ENDIF
1087      !                                        !* bottom temperature, salinity, velocity and depth
1088#if defined key_vectopt_loop
1089      DO jj = 1, 1   ! vector opt. (forced unrolling)
1090         DO ji = 1, jpij
1091#else
1092      DO jj = 1, jpj
1093         DO ji = 1, jpi
1094#endif
1095            ik = mbkt(ji,jj)                        ! bottom T-level index
1096            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
1097            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
1098            zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
1099            !
1100            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
1101            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
1102         END DO
1103      END DO
1104      !                                   !-------------------!
1105      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
1106         !                                !-------------------!
1107         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
1108               !
1109         CASE( 1 )                                   != use of upper velocity
1110            ! NOTE: not much needed for deriving, almost all the computations are for the SIGN, which is kept as in the NL
1111            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
1112               DO ji = 1, fs_jpim1   ! vector opt.
1113                  !                                                ! j-direction
1114                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
1115                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
1116                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
1117                  !                                                         ! masked bbl j-gradient of density
1118                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
1119                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
1120                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
1121                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
1122                  !
1123                  !                                                           ! bbl velocity
1124                  zvbad(ji,jj) = zvbad(ji,jj) + vtr_bbl_ad(ji,jj) * ( 0.5 + zsigna ) * ( 0.5 - zsign )   &
1125                               &                                  *     e1v(ji,jj) * e3v_bbl_0(ji,jj)
1126                  vtr_bbl_ad(ji,jj) = 0.0_wp
1127                  !                                               ! i-direction
1128                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
1129                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
1130                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
1131                  !                                                           ! masked bbl i-gradient of density
1132                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
1133                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
1134                  !
1135                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
1136                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
1137                  !
1138                  !                                                           ! bbl velocity
1139                  zubad(ji,jj) = zubad(ji,jj) + utr_bbl_ad(ji,jj) * ( 0.5 + zsigna ) * ( 0.5 - zsign )   &
1140                               &                                  * e2u(ji,jj) * e3u_bbl_0(ji,jj)
1141                  utr_bbl_ad(ji,jj) = 0.0_wp
1142               !
1143            END DO
1144         END DO
1145            !
1146         CASE( 2 )                                 != bbl velocity = F( delta rho )
1147            ! NOTE: this one is nastier
1148            zgbbl = grav * rn_gambbl
1149            DO jj = jpjm1, 1, -1                            ! criteria: rho_up > rho_down
1150               DO ji = fs_jpim1, 1, -1   ! vector opt.
1151                  !                                               ! j-direction
1152                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
1153                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
1154                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
1155                  !
1156                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
1157                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
1158                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
1159                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
1160                  zgdrho =    fsbeta( zt, zs, zh )                                    &
1161                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    &
1162                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
1163                  !
1164                  zsign  = SIGN( 0.5_wp, zgdrho ) ! adjoint of zgdrho = MAX( 0.e0, zgdrho )
1165                  !                                             ! bbl transport (down-slope direction)
1166                  zgdrhoad = zsign * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * vtr_bbl_ad(ji,jj) * REAL( mgrhv(ji,jj) )
1167                  vtr_bbl_ad(ji,jj) = 0.0_wp
1168
1169                  ztad = ( fsbeta_adj_t( zt, zs, zh, zgdrhoad )                        &
1170                     & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1171                     &                          - ( zsb(ji,ijd) - zsb(ji,ijs) ) )      &
1172                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_t( zt, zs, zh, zgdrhoad ) &
1173                     &                          * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1174                     &   ) * vmask(ji,jj,1)
1175                  zsad = ( fsbeta_adj_s( zt, zs, zh, zgdrhoad )                        &
1176                     & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1177                     &                          - ( zsb(ji,ijd) - zsb(ji,ijs) ) )      &
1178                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_s( zt, zs, zh, zgdrhoad ) &
1179                     &                          * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1180                     &   ) * vmask(ji,jj,1)
1181                  zhad = ( fsbeta_adj_h( zt, zs, zh, zgdrhoad )                        &
1182                     & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1183                     &                          - ( zsb(ji,ijd) - zsb(ji,ijs) ) )      &
1184                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_h( zt, zs, zh, zgdrhoad ) &
1185                     &                          * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1186                     &   ) * vmask(ji,jj,1)
1187
1188                  ztbad(ji,ijd) = ztbad(ji,ijd) + zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * vmask(ji,jj,1)
1189                  ztbad(ji,ijs) = ztbad(ji,ijs) - zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * vmask(ji,jj,1)
1190                  zsbad(ji,ijd) = zsbad(ji,ijd) - zgdrhoad * fsbeta( zt, zs, zh ) * vmask(ji,jj,1)
1191                  zsbad(ji,ijs) = zsbad(ji,ijs) + zgdrhoad * fsbeta( zt, zs, zh ) * vmask(ji,jj,1)
1192
1193                  ztbad (ji,jj  ) = ztbad (ji,jj  ) + 0.5 * ztad
1194                  ztbad (ji,jj+1) = ztbad (ji,jj+1) + 0.5 * ztad
1195                  zsbad (ji,jj  ) = zsbad (ji,jj  ) + 0.5 * zsad
1196                  zsbad (ji,jj+1) = zsbad (ji,jj+1) + 0.5 * zsad
1197                  ztad = 0.0_wp ; zsad = 0.0_wp ; zhad = 0.0_wp
1198
1199                  !                                         ! i-direction
1200                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
1201                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
1202                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
1203                  !
1204                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
1205                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
1206                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
1207                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
1208                  zgdrho =    fsbeta( zt, zs, zh )                                    &
1209                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    &
1210                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
1211                  zsign  = SIGN( 0.5_wp, zgdrho ) ! adjoint of zgdrho = MAX( 0.e0, zgdrho )
1212                  !                                             ! bbl transport (down-slope direction)
1213                  zgdrhoad = zsign * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * utr_bbl_ad(ji,jj) * REAL( mgrhu(ji,jj) )
1214                  utr_bbl_ad(ji,jj) = 0.0_wp
1215                  !
1216                  ztad = ( fsbeta_adj_t( zt, zs, zh, zgdrhoad )                        &
1217                     & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis, jj) )       &
1218                     &                          - ( zsb(iid,jj) - zsb(iis, jj) ) )     &
1219                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_t( zt, zs, zh, zgdrhoad ) &
1220                     &                          * ( ztb(iid,jj) - ztb(iis, jj) )       &
1221                     &   ) * umask(ji,jj,1)
1222                  zsad = ( fsbeta_adj_s( zt, zs, zh, zgdrhoad )                        &
1223                     & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis, jj) )       &
1224                     &                          - ( zsb(iid,jj) - zsb(iis, jj) ) )     &
1225                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_s( zt, zs, zh, zgdrhoad ) &
1226                     &                          * ( ztb(iid,jj) - ztb(iis, jj) )       &
1227                     &   ) * umask(ji,jj,1)
1228                  zhad = ( fsbeta_adj_h( zt, zs, zh, zgdrhoad )                        &
1229                     & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis, jj) )       &
1230                     &                          - ( zsb(iid,jj) - zsb(iis, jj) ) )     &
1231                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_h( zt, zs, zh, zgdrhoad ) &
1232                     &                          * ( ztb(iid,jj) - ztb(iis, jj) )       &
1233                     &   ) * umask(ji,jj,1)
1234
1235                  ztbad(iid,jj) = ztbad(iid,jj) + zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * umask(ji,jj,1)
1236                  ztbad(iis,jj) = ztbad(iis,jj) - zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * umask(ji,jj,1)
1237                  zsbad(iid,jj) = zsbad(iid,jj) - zgdrhoad * fsbeta( zt, zs, zh ) * umask(ji,jj,1)
1238                  zsbad(iis,jj) = zsbad(iis,jj) + zgdrhoad * fsbeta( zt, zs, zh ) * umask(ji,jj,1)
1239                  zgdrhoad = 0.0_wp
1240
1241                  ztbad (ji,jj  ) = ztbad (ji,jj  ) + 0.5 * ztad
1242                  ztbad (ji+1,jj) = ztbad (ji+1,jj) + 0.5 * ztad
1243                  zsbad (ji,jj  ) = zsbad (ji,jj  ) + 0.5 * zsad
1244                  zsbad (ji+1,jj) = zsbad (ji+1,jj) + 0.5 * zsad
1245                  ztad = 0.0_wp ; zsad = 0.0_wp ; zhad = 0.0_wp
1246                  !
1247               END DO
1248            END DO
1249         END SELECT
1250         !
1251      ENDIF
1252      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
1253         !                                !-------------------!
1254         ! NOTE : while rn_ahtbbl remains a passive variable, the code below will only yield ah_bbl_ad=0
1255#if defined key_control_param
1256         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
1257            DO ji = 1, jpim1
1258               !                                                ! j-direction
1259                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
1260                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
1261                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
1262                  !                                                           ! masked bbl j-gradient of density
1263                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
1264                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
1265                  !
1266               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
1267               ahv_bbl_0_ad(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_ad(ji,jj)
1268               ahv_bbl_ad(ji,jj) = 0.0_wp
1269               !
1270               !                                                ! i-direction
1271               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
1272               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
1273               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
1274               !                                                         ! masked bbl i-gradient of density
1275               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
1276                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
1277               !
1278               zsign             = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
1279               ahu_bbl_0_ad(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_ad(ji,jj)             ! masked diffusive flux coeff.
1280               ahu_bbl_ad(ji,jj) = 0.0_wp
1281               !
1282               END DO
1283            END DO
1284#else
1285         DO jj = 1, jpjm1
1286            DO ji = 1, jpim1
1287               ahu_bbl_ad(ji,jj)=0.0_wp
1288               ahv_bbl_ad(ji,jj)=0.0_wp
1289            END DO
1290         END DO
1291#endif
1292            !
1293      ENDIF
1294      !                                        !* bottom temperature, salinity, velocity and depth
1295#if defined key_vectopt_loop
1296      DO jj = 1, 1   ! vector opt. (forced unrolling)
1297         DO ji = 1, jpij
1298#else
1299      DO jj = 1, jpj
1300         DO ji = 1, jpi
1301#endif
1302            ik = mbkt(ji,jj)                        ! bottom T-level index
1303            tsb_ad(ji,jj,ik,jp_tem) = tsb_ad(ji,jj,ik,jp_tem) + ztbad(ji,jj) * tmask(ji,jj,1)
1304            tsb_ad(ji,jj,ik,jp_sal) = tsb_ad(ji,jj,ik,jp_sal) + zsbad(ji,jj) * tmask(ji,jj,1)
1305            ztbad (ji,jj) = 0.0_wp
1306            zsbad (ji,jj) = 0.0_wp
1307         END DO
1308      END DO
1309      !                                   !-------------------!
1310      !
1311      CALL wrk_dealloc( jpi, jpj, zub  , zvb  , ztb  , zsb, zdep, &
1312         &                        ztbad, zsbad, ztbad, zsbad      )
1313      !
1314      IF( nn_timing == 1 )  CALL timing_stop( 'bbl_adj')
1315      !
1316   END SUBROUTINE bbl_adj
1317
1318
1319   SUBROUTINE tra_bbl_init_tam
1320      !!----------------------------------------------------------------------
1321      !!                  ***  ROUTINE tra_bbl_init  ***
1322      !!
1323      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
1324      !!
1325      !! ** Method  :
1326      !!----------------------------------------------------------------------
1327      !
1328      integer :: ierr
1329      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init_tam')
1330      !
1331      ierr = tra_bbl_alloc_tam( 0 )
1332
1333      ahu_bbl_0_tl = 0.0_wp
1334      ahv_bbl_0_tl = 0.0_wp
1335      ahu_bbl_0_ad = 0.0_wp
1336      ahv_bbl_0_ad = 0.0_wp
1337      !
1338      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init_tam')
1339      !
1340   END SUBROUTINE tra_bbl_init_tam
1341
1342   SUBROUTINE tra_bbl_adj_tst( kumadt )
1343      !!-----------------------------------------------------------------------
1344      !!
1345      !!                  ***  ROUTINE tra_bbl_adj_tst ***
1346      !!
1347      !! ** Purpose : Test the adjoint routine.
1348      !!
1349      !! ** Method  : Verify the scalar product
1350      !!
1351      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
1352      !!
1353      !!              where  L   = tangent routine
1354      !!                     L^T = adjoint routine
1355      !!                     W   = diagonal matrix of scale factors
1356      !!                     dx  = input perturbation (random field)
1357      !!                     dy  = L dx
1358      !!
1359      !!
1360      !! History :
1361      !!        ! 08-08 (A. Vidard)
1362      !!-----------------------------------------------------------------------
1363      !! * Modules used
1364
1365      !! * Arguments
1366      INTEGER, INTENT(IN) :: &
1367         & kumadt             ! Output unit
1368
1369      !! * Local declarations
1370      INTEGER ::  &
1371         & ji,    &        ! dummy loop indices
1372         & jj,    &
1373         & jk,    &
1374         & jtst
1375      INTEGER ::  &
1376         & jsav1, &
1377         & jsav2
1378      REAL(KIND=wp) :: &
1379         & zsp1,         & ! scalar product involving the tangent routine
1380         & zsp2            ! scalar product involving the adjoint routine
1381      REAL(KIND=wp), POINTER, DIMENSION(:,:,:,:) :: &
1382         & ztsa_tlin ,     & ! Tangent input
1383         & ztsa_tlout,     &
1384         & ztsb_tlin ,     &
1385         & ztsa_adout,     & ! Adjoint input
1386         & ztsa_adin ,     &
1387         & ztsb_adout,     &
1388         & zrts            ! 2*3D random field
1389      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
1390         & zutr_tlin ,     &
1391         & zutr_tlout,     &
1392         & zvtr_tlin ,     &
1393         & zvtr_tlout,     &
1394         & zutr_adout,     &
1395         & zutr_adin ,     &
1396         & zvtr_adout,     &
1397         & zvtr_adin ,     &
1398         & zr2             ! 2D random field
1399      CHARACTER(LEN=14) :: &
1400         & cl_name
1401      ! Allocate memory
1402
1403      CALL wrk_alloc( jpi, jpj, jpk, jpts, ztsa_tlin , ztsa_tlout, ztsb_tlin , &
1404         &                                  ztsa_adout, ztsa_adin , ztsb_adout, &
1405         &                                  zrts )
1406      CALL wrk_alloc( jpi, jpj, zutr_tlin , zutr_tlout, zvtr_tlin , zvtr_tlout, &
1407         &                      zutr_adout, zutr_adin , zvtr_adout, zvtr_adin , &
1408         &                      zr2 )
1409
1410      CALL grid_random( utr_bbl(:,:), 'U', 0.0_wp, stdu )
1411      CALL grid_random( vtr_bbl(:,:), 'V', 0.0_wp, stdv )
1412
1413      jsav1 = nn_bbl_ldf
1414      jsav2 = nn_bbl_adv
1415
1416      DO jtst = 1, 2
1417         !==================================================================
1418         ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
1419         !    dy = ( hdivb_tl, hdivn_tl )
1420         !==================================================================
1421
1422         SELECT CASE( jtst)
1423         CASE ( 1 )
1424            nn_bbl_ldf = 1
1425            nn_bbl_adv = 0
1426         CASE ( 2 )
1427            nn_bbl_ldf = 0
1428            nn_bbl_adv = 1
1429         CASE ( 3 )
1430            nn_bbl_ldf = 0
1431            nn_bbl_adv = 2
1432         END SELECT
1433         !--------------------------------------------------------------------
1434         ! Reset the tangent and adjoint variables
1435         !--------------------------------------------------------------------
1436         ztsa_tlin (:,:,:,:) = 0.0_wp
1437         ztsa_tlout(:,:,:,:) = 0.0_wp
1438         ztsb_tlin (:,:,:,:) = 0.0_wp
1439         ztsa_adout(:,:,:,:) = 0.0_wp
1440         ztsa_adin (:,:,:,:) = 0.0_wp
1441         ztsb_adout(:,:,:,:) = 0.0_wp
1442
1443         zutr_tlin (:,:) = 0.0_wp
1444         zutr_tlout(:,:) = 0.0_wp
1445         zvtr_tlin (:,:) = 0.0_wp
1446         zvtr_tlout(:,:) = 0.0_wp
1447         zutr_adout(:,:) = 0.0_wp
1448         zutr_adin (:,:) = 0.0_wp
1449         zvtr_adout(:,:) = 0.0_wp
1450         zvtr_adin (:,:) = 0.0_wp
1451
1452         tsb_ad(:,:,:,:) = 0.0_wp
1453         !--------------------------------------------------------------------
1454         ! Initialize the tangent input with random noise: dx
1455         !--------------------------------------------------------------------
1456
1457         CALL grid_random( zrts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
1458         CALL grid_random( zrts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
1459         DO jk = 1, jpk
1460            DO jj = nldj, nlej
1461               DO ji = nldi, nlei
1462                  ztsa_tlin(ji,jj,jk,:) = zrts(ji,jj,jk,:)
1463               END DO
1464            END DO
1465         END DO
1466
1467         CALL grid_random( zrts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
1468         CALL grid_random( zrts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
1469         DO jk = 1, jpk
1470            DO jj = nldj, nlej
1471               DO ji = nldi, nlei
1472                  ztsb_tlin(ji,jj,jk,:) = zrts(ji,jj,jk,:)
1473               END DO
1474            END DO
1475         END DO
1476
1477         CALL grid_random( zr2(:,:), 'U', 0.0_wp, stdu )
1478         DO jj = nldj, nlej
1479            DO ji = nldi, nlei
1480               zutr_tlin(ji,jj) = zr2(ji,jj)
1481            END DO
1482         END DO
1483
1484         CALL grid_random( zr2(:,:), 'V', 0.0_wp, stdv )
1485         DO jj = nldj, nlej
1486            DO ji = nldi, nlei
1487               zvtr_tlin(ji,jj) = zr2(ji,jj)
1488            END DO
1489         END DO
1490
1491         tsa_tl(:,:,:,:) = ztsa_tlin(:,:,:,:)
1492         tsb_tl(:,:,:,:) = ztsb_tlin(:,:,:,:)
1493         utr_bbl_tl(:,:) = zutr_tlin(:,:)
1494         vtr_bbl_tl(:,:) = zvtr_tlin(:,:)
1495         CALL tra_bbl_tan ( nit000 )
1496         ztsa_tlout(:,:,:,:) = tsa_tl(:,:,:,:)
1497         zutr_tlout(:,:)     = utr_bbl_tl(:,:)
1498         zvtr_tlout(:,:)     = vtr_bbl_tl(:,:)
1499         !--------------------------------------------------------------------
1500         ! Initialize the adjoint variables: dy^* = W dy
1501         !--------------------------------------------------------------------
1502
1503         DO jk = 1, jpk
1504            DO jj = nldj, nlej
1505               DO ji = nldi, nlei
1506                  ztsa_adin(ji,jj,jk,jp_tem) = ztsa_tlout(ji,jj,jk,jp_tem) &
1507                     &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1508                     &               * tmask(ji,jj,jk)
1509                  ztsa_adin(ji,jj,jk,jp_sal) = ztsa_tlout(ji,jj,jk,jp_sal) &
1510                     &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1511                     &               * tmask(ji,jj,jk)
1512               END DO
1513            END DO
1514         END DO
1515         DO jj = nldj, nlej
1516            DO ji = nldi, nlei
1517               zutr_adin(ji,jj) = zutr_tlout(ji,jj) &
1518                  &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) &
1519                  &               * vmask(ji,jj,jk)
1520               zvtr_adin(ji,jj) = zvtr_tlout(ji,jj) &
1521                  &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) &
1522                  &               * vmask(ji,jj,jk)
1523            END DO
1524         END DO
1525         !--------------------------------------------------------------------
1526         ! Compute the scalar product: ( L dx )^T W dy
1527         !--------------------------------------------------------------------
1528
1529         zsp1 = DOT_PRODUCT( ztsa_tlout(:,:,:,jp_tem), ztsa_adin(:,:,:,jp_tem) ) &
1530            & + DOT_PRODUCT( ztsa_tlout(:,:,:,jp_sal), ztsa_adin(:,:,:,jp_sal) ) &
1531            & + DOT_PRODUCT( zutr_tlout, zutr_adin ) + DOT_PRODUCT( zvtr_tlout, zvtr_adin )
1532
1533         !--------------------------------------------------------------------
1534         ! Call the adjoint routine: dx^* = L^T dy^*
1535         !--------------------------------------------------------------------
1536
1537         tsa_ad(:,:,:,:) = ztsa_adin(:,:,:,:)
1538         utr_bbl_ad(:,:) = zutr_adin(:,:)
1539         vtr_bbl_ad(:,:) = zvtr_adin(:,:)
1540         CALL tra_bbl_adj ( nit000 )
1541         ztsa_adout(:,:,:,:) = tsa_ad(:,:,:,:)
1542         ztsb_adout(:,:,:,:) = tsb_ad(:,:,:,:)
1543         zutr_adout(:,:)     = utr_bbl_ad(:,:)
1544         zvtr_adout(:,:)     = vtr_bbl_ad(:,:)
1545
1546
1547         zsp2 = DOT_PRODUCT( ztsa_tlin(:,:,:,jp_tem), ztsa_adout(:,:,:,jp_tem) ) &
1548            & + DOT_PRODUCT( ztsa_tlin(:,:,:,jp_sal), ztsa_adout(:,:,:,jp_sal) ) &
1549            & + DOT_PRODUCT( ztsb_tlin(:,:,:,jp_tem), ztsb_adout(:,:,:,jp_tem) ) &
1550            & + DOT_PRODUCT( ztsb_tlin(:,:,:,jp_sal), ztsb_adout(:,:,:,jp_sal) ) &
1551            & + DOT_PRODUCT( zutr_tlin, zutr_adout ) + DOT_PRODUCT( zvtr_tlin, zvtr_adout )
1552
1553         SELECT CASE ( jtst )
1554         CASE ( 1 )
1555            ! 14 char:'12345678901234'
1556            cl_name = 'trabbl_adj_dif'
1557         CASE ( 2 )
1558            ! 14 char:'12345678901234'
1559            cl_name = 'trabbl_ad_adv1'
1560         CASE ( 3 )
1561            ! 14 char:'12345678901234'
1562            cl_name = 'trabbl_ad_adv2'
1563         END SELECT
1564         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1565
1566      END DO
1567
1568      CALL wrk_dealloc( jpi, jpj, jpk, jpts, ztsa_tlin , ztsa_tlout, ztsb_tlin , &
1569         &                                    ztsa_adout, ztsa_adin , ztsb_adout, &
1570         &                                    zrts )
1571      CALL wrk_dealloc( jpi, jpj, zutr_tlin , zutr_tlout, zvtr_tlin , zvtr_tlout, &
1572         &                        zutr_adout, zutr_adin , zvtr_adout, zvtr_adin , &
1573         &                        zr2 )
1574
1575      nn_bbl_ldf = jsav1
1576      nn_bbl_adv = jsav2
1577
1578
1579   END SUBROUTINE tra_bbl_adj_tst
1580
1581#else
1582   !!----------------------------------------------------------------------
1583   !!   Dummy module :                      No bottom boundary layer scheme
1584   !!----------------------------------------------------------------------
1585CONTAINS
1586   SUBROUTINE tra_bbl_init_tam               ! Dummy routine
1587   END SUBROUTINE tra_bbl_init_tam
1588   SUBROUTINE tra_bbl_tan( kt )              ! Dummy routine
1589      WRITE(*,*) 'tra_bbl_tan: You should not have seen this print! error?', kt
1590   END SUBROUTINE tra_bbl_tan
1591   SUBROUTINE tra_bbl_adj( kt )              ! Dummy routine
1592      WRITE(*,*) 'tra_bbl_adj: You should not have seen this print! error?', kt
1593   END SUBROUTINE tra_bbl_adj
1594   SUBROUTINE tra_bbl_adj_tst( kt )              ! Dummy routine
1595      WRITE(*,*) 'tra_bbl_adj_tst: You should not have seen this print! error?', kt
1596   END SUBROUTINE tra_bbl_adj_tst
1597#endif
1598   !!======================================================================
1599END MODULE trabbl_tam
Note: See TracBrowser for help on using the repository browser.