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

Last change on this file since 3611 was 3611, checked in by pabouttier, 11 years ago

Add TAM code and ORCA2_TAM configuration - see Ticket #1007

  • Property svn:executable set to *
File size: 85.7 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 ) * e1v(ji,jj) * e3v_bbl_0(ji,jj)
1122                  vtr_bbl_ad(ji,jj) = 0.0_wp
1123                  !                                               ! i-direction
1124                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
1125                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
1126                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
1127                  !                                                           ! masked bbl i-gradient of density
1128                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
1129                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
1130                  !
1131                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
1132                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
1133                  !
1134                  !                                                           ! bbl velocity
1135                  zubad(ji,jj) = zubad(ji,jj) + utr_bbl_ad(ji,jj) * ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj)
1136                  utr_bbl_ad(ji,jj) = 0.0_wp
1137               !
1138            END DO
1139         END DO
1140            !
1141         CASE( 2 )                                 != bbl velocity = F( delta rho )
1142            ! NOTE: this one is nastier
1143            zgbbl = grav * rn_gambbl
1144            DO jj = jpjm1, 1, -1                            ! criteria: rho_up > rho_down
1145               DO ji = fs_jpim1, 1, -1   ! vector opt.
1146                  !                                               ! j-direction
1147                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
1148                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
1149                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
1150                  !
1151                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
1152                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
1153                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
1154                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
1155                  zgdrho =    fsbeta( zt, zs, zh )                                    &
1156                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    &
1157                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
1158                  !
1159                  zsign  = SIGN( 0.5_wp, zgdrho ) ! adjoint of zgdrho = MAX( 0.e0, zgdrho )
1160                  !                                             ! bbl transport (down-slope direction)
1161                  zgdrhoad = zsign * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * vtr_bbl_ad(ji,jj) * REAL( mgrhv(ji,jj) )
1162                  vtr_bbl_ad(ji,jj) = 0.0_wp
1163
1164                  ztad = ( fsbeta_adj_t( zt, zs, zh, zgdrhoad )                        &
1165                     & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1166                     &                          - ( zsb(ji,ijd) - zsb(ji,ijs) ) )      &
1167                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_t( zt, zs, zh, zgdrhoad ) &
1168                     &                          * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1169                     &   ) * vmask(ji,jj,1)
1170                  zsad = ( fsbeta_adj_s( zt, zs, zh, zgdrhoad )                        &
1171                     & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1172                     &                          - ( zsb(ji,ijd) - zsb(ji,ijs) ) )      &
1173                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_s( zt, zs, zh, zgdrhoad ) &
1174                     &                          * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1175                     &   ) * vmask(ji,jj,1)
1176                  zhad = ( fsbeta_adj_h( zt, zs, zh, zgdrhoad )                        &
1177                     & * ( fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1178                     &                          - ( zsb(ji,ijd) - zsb(ji,ijs) ) )      &
1179                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_h( zt, zs, zh, zgdrhoad ) &
1180                     &                          * ( ztb(ji,ijd) - ztb(ji,ijs) )        &
1181                     &   ) * vmask(ji,jj,1)
1182
1183                  ztbad(ji,ijd) = ztbad(ji,ijd) + zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * vmask(ji,jj,1)
1184                  ztbad(ji,ijs) = ztbad(ji,ijs) - zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * vmask(ji,jj,1)
1185                  zsbad(ji,ijd) = zsbad(ji,ijd) - zgdrhoad * fsbeta( zt, zs, zh ) * vmask(ji,jj,1)
1186                  zsbad(ji,ijs) = zsbad(ji,ijs) + zgdrhoad * fsbeta( zt, zs, zh ) * vmask(ji,jj,1)
1187
1188                  ztbad (ji,jj  ) = ztbad (ji,jj  ) + 0.5 * ztad
1189                  ztbad (ji,jj+1) = ztbad (ji,jj+1) + 0.5 * ztad
1190                  zsbad (ji,jj  ) = zsbad (ji,jj  ) + 0.5 * zsad
1191                  zsbad (ji,jj+1) = zsbad (ji,jj+1) + 0.5 * zsad
1192                  ztad = 0.0_wp ; zsad = 0.0_wp ; zhad = 0.0_wp
1193
1194                  !                                         ! i-direction
1195                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
1196                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
1197                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
1198                  !
1199                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
1200                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
1201                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
1202                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
1203                  zgdrho =    fsbeta( zt, zs, zh )                                    &
1204                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    &
1205                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
1206                  zsign  = SIGN( 0.5_wp, zgdrho ) ! adjoint of zgdrho = MAX( 0.e0, zgdrho )
1207                  !                                             ! bbl transport (down-slope direction)
1208                  zgdrhoad = zsign * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * utr_bbl_ad(ji,jj) * REAL( mgrhu(ji,jj) )
1209                  utr_bbl_ad(ji,jj) = 0.0_wp
1210                  !
1211                  ztad = ( fsbeta_adj_t( zt, zs, zh, zgdrhoad )                        &
1212                     & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis, jj) )       &
1213                     &                          - ( zsb(iid,jj) - zsb(iis, jj) ) )     &
1214                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_t( zt, zs, zh, zgdrhoad ) &
1215                     &                          * ( ztb(iid,jj) - ztb(iis, jj) )       &
1216                     &   ) * umask(ji,jj,1)
1217                  zsad = ( fsbeta_adj_s( zt, zs, zh, zgdrhoad )                        &
1218                     & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis, jj) )       &
1219                     &                          - ( zsb(iid,jj) - zsb(iis, jj) ) )     &
1220                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_s( zt, zs, zh, zgdrhoad ) &
1221                     &                          * ( ztb(iid,jj) - ztb(iis, jj) )       &
1222                     &   ) * umask(ji,jj,1)
1223                  zhad = ( fsbeta_adj_h( zt, zs, zh, zgdrhoad )                        &
1224                     & * ( fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis, jj) )       &
1225                     &                          - ( zsb(iid,jj) - zsb(iis, jj) ) )     &
1226                     & +   fsbeta( zt, zs, zh ) * fsalbt_adj_h( zt, zs, zh, zgdrhoad ) &
1227                     &                          * ( ztb(iid,jj) - ztb(iis, jj) )       &
1228                     &   ) * umask(ji,jj,1)
1229
1230                  ztbad(iid,jj) = ztbad(iid,jj) + zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * umask(ji,jj,1)
1231                  ztbad(iis,jj) = ztbad(iis,jj) - zgdrhoad * fsbeta( zt, zs, zh ) * fsalbt( zt, zs, zh ) * umask(ji,jj,1)
1232                  zsbad(iid,jj) = zsbad(iid,jj) - zgdrhoad * fsbeta( zt, zs, zh ) * umask(ji,jj,1)
1233                  zsbad(iis,jj) = zsbad(iis,jj) + zgdrhoad * fsbeta( zt, zs, zh ) * umask(ji,jj,1)
1234                  zgdrhoad = 0.0_wp
1235
1236                  ztbad (ji,jj  ) = ztbad (ji,jj  ) + 0.5 * ztad
1237                  ztbad (ji+1,jj) = ztbad (ji+1,jj) + 0.5 * ztad
1238                  zsbad (ji,jj  ) = zsbad (ji,jj  ) + 0.5 * zsad
1239                  zsbad (ji+1,jj) = zsbad (ji+1,jj) + 0.5 * zsad
1240                  ztad = 0.0_wp ; zsad = 0.0_wp ; zhad = 0.0_wp
1241                  !
1242               END DO
1243            END DO
1244         END SELECT
1245         !
1246      ENDIF
1247      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
1248         !                                !-------------------!
1249         ! NOTE : while rn_ahtbbl remains a passive variable, the code below will only yield ah_bbl_ad=0
1250#if defined key_control_param
1251         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
1252            DO ji = 1, jpim1
1253               !                                                ! j-direction
1254                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
1255                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
1256                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
1257                  !                                                           ! masked bbl j-gradient of density
1258                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
1259                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
1260                  !
1261               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
1262               ahv_bbl_0_ad(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_ad(ji,jj)
1263               ahv_bbl_ad(ji,jj) = 0.0_wp
1264               !
1265               !                                                ! i-direction
1266               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
1267               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
1268               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
1269               !                                                         ! masked bbl i-gradient of density
1270               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
1271                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
1272               !
1273               zsign             = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
1274               ahu_bbl_0_ad(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_ad(ji,jj)             ! masked diffusive flux coeff.
1275               ahu_bbl_ad(ji,jj) = 0.0_wp
1276               !
1277               END DO
1278            END DO
1279#else
1280         DO jj = 1, jpjm1
1281            DO ji = 1, jpim1
1282               ahu_bbl_ad(ji,jj)=0.0_wp
1283               ahv_bbl_ad(ji,jj)=0.0_wp
1284            END DO
1285         END DO
1286#endif
1287            !
1288      ENDIF
1289      !                                        !* bottom temperature, salinity, velocity and depth
1290#if defined key_vectopt_loop
1291      DO jj = 1, 1   ! vector opt. (forced unrolling)
1292         DO ji = 1, jpij
1293#else
1294      DO jj = 1, jpj
1295         DO ji = 1, jpi
1296#endif
1297            ik = mbkt(ji,jj)                        ! bottom T-level index
1298            tsb_ad(ji,jj,ik,jp_tem) = tsb_ad(ji,jj,ik,jp_tem) + ztbad(ji,jj) * tmask(ji,jj,1)
1299            tsb_ad(ji,jj,ik,jp_sal) = tsb_ad(ji,jj,ik,jp_sal) + zsbad(ji,jj) * tmask(ji,jj,1)
1300            ztbad (ji,jj) = 0.0_wp
1301            zsbad (ji,jj) = 0.0_wp
1302         END DO
1303      END DO
1304      !                                   !-------------------!
1305      !
1306      CALL wrk_dealloc( jpi, jpj, zub  , zvb  , ztb  , zsb, zdep, &
1307         &                        ztbad, zsbad, ztbad, zsbad      )
1308      !
1309      IF( nn_timing == 1 )  CALL timing_stop( 'bbl_adj')
1310      !
1311   END SUBROUTINE bbl_adj
1312
1313
1314   SUBROUTINE tra_bbl_init_tam
1315      !!----------------------------------------------------------------------
1316      !!                  ***  ROUTINE tra_bbl_init  ***
1317      !!
1318      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
1319      !!
1320      !! ** Method  :
1321      !!----------------------------------------------------------------------
1322      !
1323      integer :: ierr
1324      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init_tam')
1325      !
1326      ierr = tra_bbl_alloc_tam( 0 )
1327
1328      ahu_bbl_0_tl = 0.0_wp
1329      ahv_bbl_0_tl = 0.0_wp
1330      ahu_bbl_0_ad = 0.0_wp
1331      ahv_bbl_0_ad = 0.0_wp
1332      !
1333      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init_tam')
1334      !
1335   END SUBROUTINE tra_bbl_init_tam
1336
1337   SUBROUTINE tra_bbl_adj_tst( kumadt )
1338      !!-----------------------------------------------------------------------
1339      !!
1340      !!                  ***  ROUTINE tra_bbl_adj_tst ***
1341      !!
1342      !! ** Purpose : Test the adjoint routine.
1343      !!
1344      !! ** Method  : Verify the scalar product
1345      !!
1346      !!                 ( L dx )^T W dy  =  dx^T L^T W dy
1347      !!
1348      !!              where  L   = tangent routine
1349      !!                     L^T = adjoint routine
1350      !!                     W   = diagonal matrix of scale factors
1351      !!                     dx  = input perturbation (random field)
1352      !!                     dy  = L dx
1353      !!
1354      !!
1355      !! History :
1356      !!        ! 08-08 (A. Vidard)
1357      !!-----------------------------------------------------------------------
1358      !! * Modules used
1359
1360      !! * Arguments
1361      INTEGER, INTENT(IN) :: &
1362         & kumadt             ! Output unit
1363
1364      !! * Local declarations
1365      INTEGER ::  &
1366         & ji,    &        ! dummy loop indices
1367         & jj,    &
1368         & jk,    &
1369         & jtst
1370      INTEGER ::  &
1371         & jsav1, &
1372         & jsav2
1373      REAL(KIND=wp) :: &
1374         & zsp1,         & ! scalar product involving the tangent routine
1375         & zsp2            ! scalar product involving the adjoint routine
1376      REAL(KIND=wp), POINTER, DIMENSION(:,:,:,:) :: &
1377         & ztsa_tlin ,     & ! Tangent input
1378         & ztsa_tlout,     &
1379         & ztsb_tlin ,     &
1380         & ztsa_adout,     & ! Adjoint input
1381         & ztsa_adin ,     &
1382         & ztsb_adout,     &
1383         & zrts            ! 2*3D random field
1384      REAL(KIND=wp), POINTER, DIMENSION(:,:) :: &
1385         & zutr_tlin ,     &
1386         & zutr_tlout,     &
1387         & zvtr_tlin ,     &
1388         & zvtr_tlout,     &
1389         & zutr_adout,     &
1390         & zutr_adin ,     &
1391         & zvtr_adout,     &
1392         & zvtr_adin ,     &
1393         & zr2             ! 2D random field
1394      CHARACTER(LEN=14) :: &
1395         & cl_name
1396      ! Allocate memory
1397
1398      CALL wrk_alloc( jpi, jpj, jpk, jpts, ztsa_tlin , ztsa_tlout, ztsb_tlin , &
1399         &                                  ztsa_adout, ztsa_adin , ztsb_adout, &
1400         &                                  zrts )
1401      CALL wrk_alloc( jpi, jpj, zutr_tlin , zutr_tlout, zvtr_tlin , zvtr_tlout, &
1402         &                      zutr_adout, zutr_adin , zvtr_adout, zvtr_adin , &
1403         &                      zr2 )
1404
1405      CALL grid_random( utr_bbl(:,:), 'U', 0.0_wp, stdu )
1406      CALL grid_random( vtr_bbl(:,:), 'V', 0.0_wp, stdv )
1407
1408      jsav1 = nn_bbl_ldf
1409      jsav2 = nn_bbl_adv
1410
1411      DO jtst = 1, 2
1412         !==================================================================
1413         ! 1) dx = ( un_tl, vn_tl, hdivn_tl ) and
1414         !    dy = ( hdivb_tl, hdivn_tl )
1415         !==================================================================
1416
1417         SELECT CASE( jtst)
1418         CASE ( 1 )
1419            nn_bbl_ldf = 1
1420            nn_bbl_adv = 0
1421         CASE ( 2 )
1422            nn_bbl_ldf = 0
1423            nn_bbl_adv = 1
1424         CASE ( 3 )
1425            nn_bbl_ldf = 0
1426            nn_bbl_adv = 2
1427         END SELECT
1428         !--------------------------------------------------------------------
1429         ! Reset the tangent and adjoint variables
1430         !--------------------------------------------------------------------
1431         ztsa_tlin (:,:,:,:) = 0.0_wp
1432         ztsa_tlout(:,:,:,:) = 0.0_wp
1433         ztsb_tlin (:,:,:,:) = 0.0_wp
1434         ztsa_adout(:,:,:,:) = 0.0_wp
1435         ztsa_adin (:,:,:,:) = 0.0_wp
1436         ztsb_adout(:,:,:,:) = 0.0_wp
1437
1438         zutr_tlin (:,:) = 0.0_wp
1439         zutr_tlout(:,:) = 0.0_wp
1440         zvtr_tlin (:,:) = 0.0_wp
1441         zvtr_tlout(:,:) = 0.0_wp
1442         zutr_adout(:,:) = 0.0_wp
1443         zutr_adin (:,:) = 0.0_wp
1444         zvtr_adout(:,:) = 0.0_wp
1445         zvtr_adin (:,:) = 0.0_wp
1446
1447         tsb_ad(:,:,:,:) = 0.0_wp
1448         !--------------------------------------------------------------------
1449         ! Initialize the tangent input with random noise: dx
1450         !--------------------------------------------------------------------
1451
1452         CALL grid_random( zrts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
1453         CALL grid_random( zrts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
1454         DO jk = 1, jpk
1455            DO jj = nldj, nlej
1456               DO ji = nldi, nlei
1457                  ztsa_tlin(ji,jj,jk,:) = zrts(ji,jj,jk,:)
1458               END DO
1459            END DO
1460         END DO
1461
1462         CALL grid_random( zrts(:,:,:,jp_tem), 'T', 0.0_wp, stdt )
1463         CALL grid_random( zrts(:,:,:,jp_sal), 'T', 0.0_wp, stds )
1464         DO jk = 1, jpk
1465            DO jj = nldj, nlej
1466               DO ji = nldi, nlei
1467                  ztsb_tlin(ji,jj,jk,:) = zrts(ji,jj,jk,:)
1468               END DO
1469            END DO
1470         END DO
1471
1472         CALL grid_random( zr2(:,:), 'U', 0.0_wp, stdu )
1473         DO jj = nldj, nlej
1474            DO ji = nldi, nlei
1475               zutr_tlin(ji,jj) = zr2(ji,jj)
1476            END DO
1477         END DO
1478
1479         CALL grid_random( zr2(:,:), 'V', 0.0_wp, stdv )
1480         DO jj = nldj, nlej
1481            DO ji = nldi, nlei
1482               zvtr_tlin(ji,jj) = zr2(ji,jj)
1483            END DO
1484         END DO
1485
1486         tsa_tl(:,:,:,:) = ztsa_tlin(:,:,:,:)
1487         tsb_tl(:,:,:,:) = ztsb_tlin(:,:,:,:)
1488         utr_bbl_tl(:,:) = zutr_tlin(:,:)
1489         vtr_bbl_tl(:,:) = zvtr_tlin(:,:)
1490         CALL tra_bbl_tan ( nit000 )
1491         ztsa_tlout(:,:,:,:) = tsa_tl(:,:,:,:)
1492         zutr_tlout(:,:)     = utr_bbl_tl(:,:)
1493         zvtr_tlout(:,:)     = vtr_bbl_tl(:,:)
1494         !--------------------------------------------------------------------
1495         ! Initialize the adjoint variables: dy^* = W dy
1496         !--------------------------------------------------------------------
1497
1498         DO jk = 1, jpk
1499            DO jj = nldj, nlej
1500               DO ji = nldi, nlei
1501                  ztsa_adin(ji,jj,jk,jp_tem) = ztsa_tlout(ji,jj,jk,jp_tem) &
1502                     &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1503                     &               * tmask(ji,jj,jk)
1504                  ztsa_adin(ji,jj,jk,jp_sal) = ztsa_tlout(ji,jj,jk,jp_sal) &
1505                     &               * e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) &
1506                     &               * tmask(ji,jj,jk)
1507               END DO
1508            END DO
1509         END DO
1510         DO jj = nldj, nlej
1511            DO ji = nldi, nlei
1512               zutr_adin(ji,jj) = zutr_tlout(ji,jj) &
1513                  &               * e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,1) &
1514                  &               * vmask(ji,jj,jk)
1515               zvtr_adin(ji,jj) = zvtr_tlout(ji,jj) &
1516                  &               * e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,1) &
1517                  &               * vmask(ji,jj,jk)
1518            END DO
1519         END DO
1520         !--------------------------------------------------------------------
1521         ! Compute the scalar product: ( L dx )^T W dy
1522         !--------------------------------------------------------------------
1523
1524         zsp1 = DOT_PRODUCT( ztsa_tlout(:,:,:,jp_tem), ztsa_adin(:,:,:,jp_tem) ) &
1525            & + DOT_PRODUCT( ztsa_tlout(:,:,:,jp_sal), ztsa_adin(:,:,:,jp_sal) ) &
1526            & + DOT_PRODUCT( zutr_tlout, zutr_adin ) + DOT_PRODUCT( zvtr_tlout, zvtr_adin )
1527
1528         !--------------------------------------------------------------------
1529         ! Call the adjoint routine: dx^* = L^T dy^*
1530         !--------------------------------------------------------------------
1531
1532         tsa_ad(:,:,:,:) = ztsa_adin(:,:,:,:)
1533         utr_bbl_ad(:,:) = zutr_adin(:,:)
1534         vtr_bbl_ad(:,:) = zvtr_adin(:,:)
1535         CALL tra_bbl_adj ( nit000 )
1536         ztsa_adout(:,:,:,:) = tsa_ad(:,:,:,:)
1537         ztsb_adout(:,:,:,:) = tsb_ad(:,:,:,:)
1538         zutr_adout(:,:)     = utr_bbl_ad(:,:)
1539         zvtr_adout(:,:)     = vtr_bbl_ad(:,:)
1540
1541
1542         zsp2 = DOT_PRODUCT( ztsa_tlin(:,:,:,jp_tem), ztsa_adout(:,:,:,jp_tem) ) &
1543            & + DOT_PRODUCT( ztsa_tlin(:,:,:,jp_sal), ztsa_adout(:,:,:,jp_sal) ) &
1544            & + DOT_PRODUCT( ztsb_tlin(:,:,:,jp_tem), ztsb_adout(:,:,:,jp_tem) ) &
1545            & + DOT_PRODUCT( ztsb_tlin(:,:,:,jp_sal), ztsb_adout(:,:,:,jp_sal) ) &
1546            & + DOT_PRODUCT( zutr_tlin, zutr_adout ) + DOT_PRODUCT( zvtr_tlin, zvtr_adout )
1547
1548         SELECT CASE ( jtst )
1549         CASE ( 1 )
1550            ! 14 char:'12345678901234'
1551            cl_name = 'trabbl_adj_dif'
1552         CASE ( 2 )
1553            ! 14 char:'12345678901234'
1554            cl_name = 'trabbl_ad_adv1'
1555         CASE ( 3 )
1556            ! 14 char:'12345678901234'
1557            cl_name = 'trabbl_ad_adv2'
1558         END SELECT
1559         CALL prntst_adj( cl_name, kumadt, zsp1, zsp2 )
1560
1561      END DO
1562
1563      CALL wrk_dealloc( jpi, jpj, jpk, jpts, ztsa_tlin , ztsa_tlout, ztsb_tlin , &
1564         &                                    ztsa_adout, ztsa_adin , ztsb_adout, &
1565         &                                    zrts )
1566      CALL wrk_dealloc( jpi, jpj, zutr_tlin , zutr_tlout, zvtr_tlin , zvtr_tlout, &
1567         &                        zutr_adout, zutr_adin , zvtr_adout, zvtr_adin , &
1568         &                        zr2 )
1569
1570      nn_bbl_ldf = jsav1
1571      nn_bbl_adv = jsav2
1572
1573
1574   END SUBROUTINE tra_bbl_adj_tst
1575
1576#else
1577   !!----------------------------------------------------------------------
1578   !!   Dummy module :                      No bottom boundary layer scheme
1579   !!----------------------------------------------------------------------
1580CONTAINS
1581   SUBROUTINE tra_bbl_init_tam               ! Dummy routine
1582   END SUBROUTINE tra_bbl_init_tam
1583   SUBROUTINE tra_bbl_tan( kt )              ! Dummy routine
1584      WRITE(*,*) 'tra_bbl_tan: You should not have seen this print! error?', kt
1585   END SUBROUTINE tra_bbl_tan
1586   SUBROUTINE tra_bbl_adj( kt )              ! Dummy routine
1587      WRITE(*,*) 'tra_bbl_adj: You should not have seen this print! error?', kt
1588   END SUBROUTINE tra_bbl_adj
1589   SUBROUTINE tra_bbl_adj_tst( kt )              ! Dummy routine
1590      WRITE(*,*) 'tra_bbl_adj_tst: You should not have seen this print! error?', kt
1591   END SUBROUTINE tra_bbl_adj_tst
1592#endif
1593   !!======================================================================
1594END MODULE trabbl_tam
Note: See TracBrowser for help on using the repository browser.