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.F90 in branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_CNRS1_10_TEOS10_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 4619

Last change on this file since 4619 was 4619, checked in by gm, 10 years ago

#1294 : TEOS-10 and Ediag

  • Property svn:keywords set to Id
File size: 34.3 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               zptb(ji,jj) = ptb(ji,jj,mbkt(ji,jj),jn)       ! bottom before T and S
201            END DO
202         END DO
203         !               
204         DO jj = 2, jpjm1                                    ! Compute the trend
205            DO ji = 2, jpim1
206               ik = mbkt(ji,jj)                              ! bottom T-level index
207               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
208               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
209                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
210                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
211                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
212                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
213            END DO
214         END DO
215         !                                                  ! ===========
216      END DO                                                ! end tracer
217      !                                                     ! ===========
218      CALL wrk_dealloc( jpi, jpj, zptb )
219      !
220      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
221      !
222   END SUBROUTINE tra_bbl_dif
223
224
225   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
226      !!----------------------------------------------------------------------
227      !!                  ***  ROUTINE trc_bbl  ***
228      !!
229      !! ** Purpose :   Compute the before passive tracer trend associated
230      !!     with the bottom boundary layer and add it to the general trend
231      !!     of tracer equations.
232      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
233      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
234      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
235      !!                       transport proportional to the along-slope density gradient
236      !!
237      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
238      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
239      !!----------------------------------------------------------------------
240      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
241      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
242      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
243      !
244      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
245      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
246      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
247      REAL(wp) ::   zbtr, ztra               ! local scalars
248      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
249      !!----------------------------------------------------------------------
250      !
251      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
252      !                                                          ! ===========
253      DO jn = 1, kjpt                                            ! tracer loop
254         !                                                       ! ===========
255         DO jj = 1, jpjm1
256            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
257               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
258                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
259                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
260                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
261                  zu_bbl = ABS( utr_bbl(ji,jj) )
262                  !
263                  !                                               ! up  -slope T-point (shelf bottom point)
264                  zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus)
265                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
266                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
267                  !
268                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
269                     zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk)
270                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
271                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
272                  END DO
273                  !
274                  zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud)
275                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
276                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
277               ENDIF
278               !
279               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
280                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
281                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
282                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
283                  zv_bbl = ABS( vtr_bbl(ji,jj) )
284                  !
285                  ! up  -slope T-point (shelf bottom point)
286                  zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs)
287                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
288                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
289                  !
290                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
291                     zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk)
292                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
293                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
294                  END DO
295                  !                                               ! down-slope T-point (deep bottom point)
296                  zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd)
297                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
298                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
299               ENDIF
300            END DO
301            !
302         END DO
303         !                                                  ! ===========
304      END DO                                                ! end tracer
305      !                                                     ! ===========
306      !
307      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
308      !
309   END SUBROUTINE tra_bbl_adv
310
311
312   SUBROUTINE bbl( kt, kit000, cdtype )
313      !!----------------------------------------------------------------------
314      !!                  ***  ROUTINE bbl  ***
315      !!
316      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
317      !!                advection terms.
318      !!
319      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
320      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
321      !!      along bottom slope gradient) an additional lateral 2nd order
322      !!      diffusion along the bottom slope is added to the general
323      !!      tracer trend, otherwise the additional trend is set to 0.
324      !!      A typical value of ahbt is 2000 m2/s (equivalent to
325      !!      a downslope velocity of 20 cm/s if the condition for slope
326      !!      convection is satified)
327      !!              * advective bbl (nn_bbl_adv=1 or 2) :
328      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
329      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
330      !!        i.e. transport proportional to the along-slope density gradient
331      !!
332      !!      NB: the along slope density gradient is evaluated using the
333      !!      local density (i.e. referenced at a common local depth).
334      !!
335      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
336      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
337      !!----------------------------------------------------------------------
338      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
339      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
340      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
341      !!
342      INTEGER  ::   ji, jj                    ! dummy loop indices
343      INTEGER  ::   ik                        ! local integers
344      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
345      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
346      REAL(wp) ::   za, zb, zgdrho            ! local scalars
347      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
348      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
349      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
350      !!----------------------------------------------------------------------
351      !
352      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
353      !
354      IF( kt == kit000 )  THEN
355         IF(lwp)  WRITE(numout,*)
356         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
357         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
358      ENDIF
359      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
360      DO jj = 1, jpj
361         DO ji = 1, jpi
362            ik = mbkt(ji,jj)                             ! bottom T-level index
363            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
364            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
365            !
366            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth
367            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
368            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
369         END DO
370      END DO
371      !
372      CALL eos_rab( zts, zdep, zab )
373      !
374      !                                   !-------------------!
375      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
376         !                                !-------------------!
377         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
378            DO ji = 1, fs_jpim1   ! vector opt.
379               !                                                   ! i-direction
380               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
381               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
382               !                                                         ! 2*masked bottom density gradient
383               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
384                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
385               !
386               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
387               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
388               !
389               !                                                   ! j-direction
390               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
391               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
392               !                                                         ! 2*masked bottom density gradient
393               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
394                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
395               !
396               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
397               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
398            END DO
399         END DO
400         !
401      ENDIF
402
403      !                                   !-------------------!
404      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
405         !                                !-------------------!
406         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
407         !
408         CASE( 1 )                                   != use of upper velocity
409            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
410               DO ji = 1, fs_jpim1   ! vector opt.
411                  !                                                  ! i-direction
412                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
413                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
414                  !                                                          ! 2*masked bottom density gradient
415                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
416                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
417                  !
418                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
419                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
420                  !
421                  !                                                          ! bbl velocity
422                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
423                  !
424                  !                                                  ! j-direction
425                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
426                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
427                  !                                                          ! 2*masked bottom density gradient
428                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
429                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
430                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
431                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
432                  !
433                  !                                                          ! bbl transport
434                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
435               END DO
436            END DO
437            !
438         CASE( 2 )                                 != bbl velocity = F( delta rho )
439            zgbbl = grav * rn_gambbl
440            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
441               DO ji = 1, fs_jpim1   ! vector opt.
442                  !                                                  ! i-direction
443                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
444                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
445                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
446                  !
447                  ikud = mbku_d(ji,jj)
448                  ikus = mbku(ji,jj)
449                  !
450                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
451                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
452                  !                                                          !   masked bottom density gradient
453                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
454                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
455                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
456                  !
457                  !                                                          ! bbl transport (down-slope direction)
458                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
459                  !
460                  !                                                  ! j-direction
461                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
462                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
463                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
464                  !
465                  ikvd = mbkv_d(ji,jj)
466                  ikvs = mbkv(ji,jj)
467                  !
468                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
469                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
470                  !                                                          !   masked bottom density gradient
471                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
472                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
473                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
474                  !
475                  !                                                          ! bbl transport (down-slope direction)
476                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
477               END DO
478            END DO
479         END SELECT
480         !
481      ENDIF
482      !
483      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
484      !
485   END SUBROUTINE bbl
486
487
488   SUBROUTINE tra_bbl_init
489      !!----------------------------------------------------------------------
490      !!                  ***  ROUTINE tra_bbl_init  ***
491      !!
492      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
493      !!
494      !! ** Method  :   Read the nambbl namelist and check the parameters
495      !!              called by nemo_init at the first timestep (kit000)
496      !!----------------------------------------------------------------------
497      INTEGER ::   ji, jj               ! dummy loop indices
498      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
499      INTEGER ::   ios                  !   -      -
500      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
501      !!
502      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
503      !!----------------------------------------------------------------------
504      !
505      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
506      !
507      CALL wrk_alloc( jpi, jpj, zmbk )
508      !
509
510      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
511      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
512901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
513
514      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
515      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
516902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
517      WRITE ( numond, nambbl )
518      !
519      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
520      !
521      IF(lwp) THEN                   !* Parameter control and print
522         WRITE(numout,*)
523         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
524         WRITE(numout,*) '~~~~~~~~~~~~'
525         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
526         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
527         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
528         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
529         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
530      ENDIF
531
532      !                              ! allocate trabbl arrays
533      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
534
535      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
536      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
537
538      !                             !* vertical index of  "deep" bottom u- and v-points
539      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
540         DO ji = 1, jpim1
541            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
542            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
543         END DO
544      END DO
545      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
546      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
547      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
548
549                                        !* sign of grad(H) at u- and v-points
550      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
551      DO jj = 1, jpjm1
552         DO ji = 1, jpim1
553            mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
554            mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
555         END DO
556      END DO
557
558      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
559         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
560            e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) )
561            e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) )
562         END DO
563      END DO
564      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
565
566      !                             !* masked diffusive flux coefficients
567      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
568      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
569
570
571      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
572         !
573         SELECT CASE ( jp_cfg )
574         CASE ( 2 )                          ! ORCA_R2
575            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
576            ii0 = 139   ;   ii1 = 140
577            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
578            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
579            !
580            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
581            ii0 = 161   ;   ii1 = 162
582            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
583            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
584            !
585         CASE ( 4 )                          ! ORCA_R4
586            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
587            ii0 =  70   ;   ii1 =  71
588            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
589            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
590         END SELECT
591         !
592      ENDIF
593      !
594      CALL wrk_dealloc( jpi, jpj, zmbk )
595      !
596      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
597      !
598   END SUBROUTINE tra_bbl_init
599
600#else
601   !!----------------------------------------------------------------------
602   !!   Dummy module :                      No bottom boundary layer scheme
603   !!----------------------------------------------------------------------
604   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
605CONTAINS
606   SUBROUTINE tra_bbl_init               ! Dummy routine
607   END SUBROUTINE tra_bbl_init
608   SUBROUTINE tra_bbl( kt )              ! Dummy routine
609      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
610   END SUBROUTINE tra_bbl
611#endif
612
613   !!======================================================================
614END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.