source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 8507

Last change on this file since 8507 was 8507, checked in by acc, 3 years ago

Branch 2015/nemo_v3_6_STABLE. Bug fix for #1932 to correct tracer conservation at the north-fold

  • Property svn:keywords set to Id
File size: 34.6 KB
Line 
1MODULE trabbl
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   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
15   !!----------------------------------------------------------------------
16#if   defined key_trabbl   ||   defined key_esopa
17   !!----------------------------------------------------------------------
18   !!   'key_trabbl'   or                             bottom boundary layer
19   !!----------------------------------------------------------------------
20   !!   tra_bbl_alloc : allocate trabbl arrays
21   !!   tra_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
22   !!   tra_bbl_dif   : generic routine to compute bbl diffusive trend
23   !!   tra_bbl_adv   : generic routine to compute bbl advective trend
24   !!   bbl           : computation of bbl diffu. flux coef. & transport in bottom boundary layer
25   !!   tra_bbl_init  : initialization, namelist read, parameters control
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and active tracers
28   USE dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constant
30   USE eosbn2         ! equation of state
31   USE trd_oce     ! trends: ocean variables
32   USE trdtra         ! trends: active tracers
33   !
34   USE iom            ! IOM library               
35   USE in_out_manager ! I/O manager
36   USE lbclnk         ! ocean lateral boundary conditions
37   USE prtctl         ! Print control
38   USE wrk_nemo       ! Memory Allocation
39   USE timing         ! Timing
40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC   tra_bbl       !  routine called by step.F90
46   PUBLIC   tra_bbl_init  !  routine called by opa.F90
47   PUBLIC   tra_bbl_dif   !  routine called by trcbbl.F90
48   PUBLIC   tra_bbl_adv   !  -          -          -              -
49   PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
50
51   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
52
53   !                                !!* Namelist nambbl *
54   INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
55   INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
56   !                                            !  =1 : advective bbl using the bottom ocean velocity
57   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
58   REAL(wp), PUBLIC ::   rn_ahtbbl   !: along slope bbl diffusive coefficient [m2/s]
59   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
60
61   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
62
63   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
64   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
65
66   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
67   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   e3u_bbl_0, e3v_bbl_0   ! thichness of the bbl (e3) at u and v-points (PUBLIC for TAM)
70
71   !! * Substitutions
72#  include "domzgr_substitute.h90"
73#  include "vectopt_loop_substitute.h90"
74   !!----------------------------------------------------------------------
75   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
76   !! $Id$
77   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   INTEGER FUNCTION tra_bbl_alloc()
82      !!----------------------------------------------------------------------
83      !!                ***  FUNCTION tra_bbl_alloc  ***
84      !!----------------------------------------------------------------------
85      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
86         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
87         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
88         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
89         !
90      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
91      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
92   END FUNCTION tra_bbl_alloc
93
94
95   SUBROUTINE tra_bbl( kt )
96      !!----------------------------------------------------------------------
97      !!                  ***  ROUTINE bbl  ***
98      !!
99      !! ** Purpose :   Compute the before tracer (t & s) trend associated
100      !!              with the bottom boundary layer and add it to the general
101      !!              trend of tracer equations.
102      !!
103      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
104      !!              diffusive and/or advective contribution to the tracer trend
105      !!              is added to the general tracer trend
106      !!----------------------------------------------------------------------
107      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
108      !
109      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
110      !!----------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
113      !
114      IF( l_trdtra )   THEN                         !* Save ta and sa trends
115         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
116         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
117         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
118      ENDIF
119
120      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
121
122      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
123         !
124         CALL tra_bbl_dif( tsb, tsa, jpts )
125         IF( ln_ctl )  &
126         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
127            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
128         ! lateral boundary conditions ; just need for outputs
129         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
130         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
131         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
132         !
133      END IF
134
135      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
136         !
137         CALL tra_bbl_adv( tsb, tsa, jpts )
138         IF(ln_ctl)   &
139         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
140            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
141         ! lateral boundary conditions ; just need for outputs
142         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
143         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
144         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
145         !
146      END IF
147
148      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
149         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
150         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
151         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
152         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
153         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
154      ENDIF
155      !
156      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
157      !
158   END SUBROUTINE tra_bbl
159
160
161   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
162      !!----------------------------------------------------------------------
163      !!                  ***  ROUTINE tra_bbl_dif  ***
164      !!
165      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
166      !!                advection terms.
167      !!
168      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
169      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
170      !!      along bottom slope gradient) an additional lateral 2nd order
171      !!      diffusion along the bottom slope is added to the general
172      !!      tracer trend, otherwise the additional trend is set to 0.
173      !!      A typical value of ahbt is 2000 m2/s (equivalent to
174      !!      a downslope velocity of 20 cm/s if the condition for slope
175      !!      convection is satified)
176      !!
177      !! ** Action  :   pta   increased by the bbl diffusive trend
178      !!
179      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
180      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
181      !!----------------------------------------------------------------------
182      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
183      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
184      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
185      !
186      INTEGER  ::   ji, jj, jn   ! dummy loop indices
187      INTEGER  ::   ik           ! local integers
188      REAL(wp) ::   zbtr         ! local scalars
189      REAL(wp), POINTER, DIMENSION(:,:) :: zptb
190      !!----------------------------------------------------------------------
191      !
192      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
193      !
194      CALL wrk_alloc( jpi, jpj, zptb )
195      !
196      DO jn = 1, kjpt                                     ! tracer loop
197         !                                                ! ===========
198         DO jj = 1, jpj
199            DO ji = 1, jpi
200               ik = mbkt(ji,jj)                              ! bottom T-level index
201               zptb(ji,jj) = ptb(ji,jj,ik,jn)       ! bottom before T and S
202            END DO
203         END DO
204         !               
205         DO jj = 2, jpjm1                                    ! Compute the trend
206            DO ji = 2, jpim1
207               ik = mbkt(ji,jj)                              ! bottom T-level index
208               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
209               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
210                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
211                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
212                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
213                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
214            END DO
215         END DO
216         !                                                  ! ===========
217      END DO                                                ! end tracer
218      !                                                     ! ===========
219      CALL wrk_dealloc( jpi, jpj, zptb )
220      !
221      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
222      !
223   END SUBROUTINE tra_bbl_dif
224
225
226   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
227      !!----------------------------------------------------------------------
228      !!                  ***  ROUTINE trc_bbl  ***
229      !!
230      !! ** Purpose :   Compute the before passive tracer trend associated
231      !!     with the bottom boundary layer and add it to the general trend
232      !!     of tracer equations.
233      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
234      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
235      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
236      !!                       transport proportional to the along-slope density gradient
237      !!
238      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
239      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
240      !!----------------------------------------------------------------------
241      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
242      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
243      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
244      !
245      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
246      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
247      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
248      REAL(wp) ::   zbtr, ztra               ! local scalars
249      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
250      !!----------------------------------------------------------------------
251      !
252      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
253      !                                                          ! ===========
254      DO jn = 1, kjpt                                            ! tracer loop
255         !                                                       ! ===========
256         DO jj = 1, jpjm1
257            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
258               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
259                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
260                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
261                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
262                  zu_bbl = ABS( utr_bbl(ji,jj) )
263                  !
264                  !                                               ! up  -slope T-point (shelf bottom point)
265                  zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus)
266                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
267                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
268                  !
269                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
270                     zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk)
271                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
272                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
273                  END DO
274                  !
275                  zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud)
276                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
277                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
278               ENDIF
279               !
280               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
281                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
282                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
283                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
284                  zv_bbl = ABS( vtr_bbl(ji,jj) )
285                  !
286                  ! up  -slope T-point (shelf bottom point)
287                  zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs)
288                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
289                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
290                  !
291                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
292                     zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk)
293                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
294                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
295                  END DO
296                  !                                               ! down-slope T-point (deep bottom point)
297                  zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd)
298                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
299                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
300               ENDIF
301            END DO
302            !
303         END DO
304         !                                                  ! ===========
305      END DO                                                ! end tracer
306      !                                                     ! ===========
307      !
308      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
309      !
310   END SUBROUTINE tra_bbl_adv
311
312
313   SUBROUTINE bbl( kt, kit000, cdtype )
314      !!----------------------------------------------------------------------
315      !!                  ***  ROUTINE bbl  ***
316      !!
317      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
318      !!                advection terms.
319      !!
320      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
321      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
322      !!      along bottom slope gradient) an additional lateral 2nd order
323      !!      diffusion along the bottom slope is added to the general
324      !!      tracer trend, otherwise the additional trend is set to 0.
325      !!      A typical value of ahbt is 2000 m2/s (equivalent to
326      !!      a downslope velocity of 20 cm/s if the condition for slope
327      !!      convection is satified)
328      !!              * advective bbl (nn_bbl_adv=1 or 2) :
329      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
330      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
331      !!        i.e. transport proportional to the along-slope density gradient
332      !!
333      !!      NB: the along slope density gradient is evaluated using the
334      !!      local density (i.e. referenced at a common local depth).
335      !!
336      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
337      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
338      !!----------------------------------------------------------------------
339      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
340      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
341      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
342      !!
343      INTEGER  ::   ji, jj                    ! dummy loop indices
344      INTEGER  ::   ik                        ! local integers
345      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
346      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
347      REAL(wp) ::   za, zb, zgdrho            ! local scalars
348      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
349      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
350      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
351      !!----------------------------------------------------------------------
352      !
353      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
354      !
355      IF( kt == kit000 )  THEN
356         IF(lwp)  WRITE(numout,*)
357         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
358         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
359      ENDIF
360      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
361      DO jj = 1, jpj
362         DO ji = 1, jpi
363            ik = mbkt(ji,jj)                             ! bottom T-level index
364            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
365            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
366            !
367            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth
368            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
369            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
370         END DO
371      END DO
372      !
373      CALL eos_rab( zts, zdep, zab )
374      !
375      !                                   !-------------------!
376      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
377         !                                !-------------------!
378         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
379            DO ji = 1, fs_jpim1   ! vector opt.
380               !                                                   ! i-direction
381               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
382               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
383               !                                                         ! 2*masked bottom density gradient
384               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
385                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
386               !
387               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
388               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
389               !
390               !                                                   ! j-direction
391               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
392               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
393               !                                                         ! 2*masked bottom density gradient
394               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
395                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
396               !
397               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
398               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
399            END DO
400         END DO
401         !
402      ENDIF
403
404      !                                   !-------------------!
405      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
406         !                                !-------------------!
407         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
408         !
409         CASE( 1 )                                   != use of upper velocity
410            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
411               DO ji = 1, fs_jpim1   ! vector opt.
412                  !                                                  ! i-direction
413                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
414                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
415                  !                                                          ! 2*masked bottom density gradient
416                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
417                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
418                  !
419                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
420                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
421                  !
422                  !                                                          ! bbl velocity
423                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
424                  !
425                  !                                                  ! j-direction
426                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
427                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
428                  !                                                          ! 2*masked bottom density gradient
429                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
430                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
431                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
432                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
433                  !
434                  !                                                          ! bbl transport
435                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
436               END DO
437            END DO
438            !
439         CASE( 2 )                                 != bbl velocity = F( delta rho )
440            zgbbl = grav * rn_gambbl
441            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
442               DO ji = 1, fs_jpim1   ! vector opt.
443                  !                                                  ! i-direction
444                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
445                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
446                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
447                  !
448                  ikud = mbku_d(ji,jj)
449                  ikus = mbku(ji,jj)
450                  !
451                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
452                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
453                  !                                                          !   masked bottom density gradient
454                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
455                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
456                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
457                  !
458                  !                                                          ! bbl transport (down-slope direction)
459                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
460                  !
461                  !                                                  ! j-direction
462                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
463                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
464                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
465                  !
466                  ikvd = mbkv_d(ji,jj)
467                  ikvs = mbkv(ji,jj)
468                  !
469                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
470                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
471                  !                                                          !   masked bottom density gradient
472                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
473                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
474                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
475                  !
476                  !                                                          ! bbl transport (down-slope direction)
477                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
478               END DO
479            END DO
480         END SELECT
481         !
482      ENDIF
483      !
484      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
485      !
486   END SUBROUTINE bbl
487
488
489   SUBROUTINE tra_bbl_init
490      !!----------------------------------------------------------------------
491      !!                  ***  ROUTINE tra_bbl_init  ***
492      !!
493      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
494      !!
495      !! ** Method  :   Read the nambbl namelist and check the parameters
496      !!              called by nemo_init at the first timestep (kit000)
497      !!----------------------------------------------------------------------
498      INTEGER ::   ji, jj               ! dummy loop indices
499      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
500      INTEGER ::   ios                  !   -      -
501      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
502      !!
503      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
504      !!----------------------------------------------------------------------
505      !
506      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
507      !
508      CALL wrk_alloc( jpi, jpj, zmbk )
509      !
510
511      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
512      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
513901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
514
515      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
516      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
517902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
518      IF(lwm) WRITE ( numond, nambbl )
519      !
520      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
521      !
522      IF(lwp) THEN                   !* Parameter control and print
523         WRITE(numout,*)
524         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
525         WRITE(numout,*) '~~~~~~~~~~~~'
526         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
527         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
528         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
529         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
530         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
531      ENDIF
532
533      !                              ! allocate trabbl arrays
534      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
535
536      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
537      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
538
539      !                             !* vertical index of  "deep" bottom u- and v-points
540      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
541         DO ji = 1, jpim1
542            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
543            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
544         END DO
545      END DO
546      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
547      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
548      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
549
550                                        !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
551      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
552      DO jj = 1, jpjm1
553         DO ji = 1, jpim1
554            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
555               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
556            ENDIF
557            !
558            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
559               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
560            ENDIF
561         END DO
562      END DO
563
564      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
565         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
566            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
567            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
568         END DO
569      END DO
570      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
571
572      !                             !* masked diffusive flux coefficients
573      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
574      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
575
576
577      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
578         !
579         SELECT CASE ( jp_cfg )
580         CASE ( 2 )                          ! ORCA_R2
581            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
582            ii0 = 139   ;   ii1 = 140
583            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
584            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
585            !
586            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
587            ii0 = 161   ;   ii1 = 162
588            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
589            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
590            !
591         CASE ( 4 )                          ! ORCA_R4
592            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
593            ii0 =  70   ;   ii1 =  71
594            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
595            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
596         END SELECT
597         !
598      ENDIF
599      !
600      CALL wrk_dealloc( jpi, jpj, zmbk )
601      !
602      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
603      !
604   END SUBROUTINE tra_bbl_init
605
606#else
607   !!----------------------------------------------------------------------
608   !!   Dummy module :                      No bottom boundary layer scheme
609   !!----------------------------------------------------------------------
610   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
611CONTAINS
612   SUBROUTINE tra_bbl_init               ! Dummy routine
613   END SUBROUTINE tra_bbl_init
614   SUBROUTINE tra_bbl( kt )              ! Dummy routine
615      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
616   END SUBROUTINE tra_bbl
617#endif
618
619   !!======================================================================
620END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.