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

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

Correct long lines in TAM 4.3 see ticket #1010

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