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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 5883

Last change on this file since 5883 was 5883, checked in by gm, 8 years ago

#1613: vvl by default: TRA/TRC remove optimization associated with linear free surface

  • 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
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 "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()
81      !!----------------------------------------------------------------------
82      !!                ***  FUNCTION tra_bbl_alloc  ***
83      !!----------------------------------------------------------------------
84      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
85         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
86         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
87         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
88         !
89      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
90      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
91   END FUNCTION tra_bbl_alloc
92
93
94   SUBROUTINE tra_bbl( kt )
95      !!----------------------------------------------------------------------
96      !!                  ***  ROUTINE bbl  ***
97      !!
98      !! ** Purpose :   Compute the before tracer (t & s) trend associated
99      !!              with the bottom boundary layer and add it to the general
100      !!              trend of tracer equations.
101      !!
102      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
103      !!              diffusive and/or advective contribution to the tracer trend
104      !!              is added to the general tracer trend
105      !!----------------------------------------------------------------------
106      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
107      !
108      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
109      !!----------------------------------------------------------------------
110      !
111      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
112      !
113      IF( l_trdtra )   THEN                         !* Save the input trends
114         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
115         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
116         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
117      ENDIF
118
119      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
120
121      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
122         !
123         CALL tra_bbl_dif( tsb, tsa, jpts )
124         IF( ln_ctl )  &
125         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
126            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
127         ! lateral boundary conditions ; just need for outputs
128         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
129         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
130         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
131         !
132      END IF
133      !
134      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
135         !
136         CALL tra_bbl_adv( tsb, tsa, jpts )
137         IF(ln_ctl)   &
138         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
139            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
140         ! lateral boundary conditions ; just need for outputs
141         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
142         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
143         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
144         !
145      END IF
146
147      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
148         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
149         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
150         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
151         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
152         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
153      ENDIF
154      !
155      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
156      !
157   END SUBROUTINE tra_bbl
158
159
160   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
161      !!----------------------------------------------------------------------
162      !!                  ***  ROUTINE tra_bbl_dif  ***
163      !!
164      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
165      !!                advection terms.
166      !!
167      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
168      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
169      !!      along bottom slope gradient) an additional lateral 2nd order
170      !!      diffusion along the bottom slope is added to the general
171      !!      tracer trend, otherwise the additional trend is set to 0.
172      !!      A typical value of ahbt is 2000 m2/s (equivalent to
173      !!      a downslope velocity of 20 cm/s if the condition for slope
174      !!      convection is satified)
175      !!
176      !! ** Action  :   pta   increased by the bbl diffusive trend
177      !!
178      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
179      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
180      !!----------------------------------------------------------------------
181      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
183      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
184      !
185      INTEGER  ::   ji, jj, jn   ! dummy loop indices
186      INTEGER  ::   ik           ! local integers
187      REAL(wp) ::   zbtr         ! local scalars
188      REAL(wp), POINTER, DIMENSION(:,:) :: zptb
189      !!----------------------------------------------------------------------
190      !
191      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
192      !
193      CALL wrk_alloc( jpi, jpj, zptb )
194      !
195      DO jn = 1, kjpt                                     ! tracer loop
196         !                                                ! ===========
197         DO jj = 1, jpj
198            DO ji = 1, jpi
199               ik = mbkt(ji,jj)                             ! bottom T-level index
200               zptb(ji,jj) = ptb(ji,jj,ik,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               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  &
208                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
209                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
210                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
211                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
212                  &             / ( e1e2t(ji,jj) * e3t_n(ji,jj,ik) )
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_e1e2t(iis,jj) / e3t_n(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_e1e2t(iid,jj) / e3t_n(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_e1e2t(iid,jj) / e3t_n(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_e1e2t(ji,ijs) / e3t_n(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_e1e2t(ji,ijd) / e3t_n(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_e1e2t(ji,ijd) / e3t_n(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      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
307      !
308   END SUBROUTINE tra_bbl_adv
309
310
311   SUBROUTINE bbl( kt, kit000, cdtype )
312      !!----------------------------------------------------------------------
313      !!                  ***  ROUTINE bbl  ***
314      !!
315      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
316      !!                advection terms.
317      !!
318      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
319      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
320      !!      along bottom slope gradient) an additional lateral 2nd order
321      !!      diffusion along the bottom slope is added to the general
322      !!      tracer trend, otherwise the additional trend is set to 0.
323      !!      A typical value of ahbt is 2000 m2/s (equivalent to
324      !!      a downslope velocity of 20 cm/s if the condition for slope
325      !!      convection is satified)
326      !!              * advective bbl (nn_bbl_adv=1 or 2) :
327      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
328      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
329      !!        i.e. transport proportional to the along-slope density gradient
330      !!
331      !!      NB: the along slope density gradient is evaluated using the
332      !!      local density (i.e. referenced at a common local depth).
333      !!
334      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
335      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
336      !!----------------------------------------------------------------------
337      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
338      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
339      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
340      !
341      INTEGER  ::   ji, jj                    ! dummy loop indices
342      INTEGER  ::   ik                        ! local integers
343      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
344      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
345      REAL(wp) ::   za, zb, zgdrho            ! local scalars
346      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
347      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
348      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
349      !!----------------------------------------------------------------------
350      !
351      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
352      !
353      IF( kt == kit000 )  THEN
354         IF(lwp)  WRITE(numout,*)
355         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
356         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
357      ENDIF
358      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
359      DO jj = 1, jpj
360         DO ji = 1, jpi
361            ik = mbkt(ji,jj)                             ! bottom T-level index
362            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
363            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
364            !
365            zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth
366            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
367            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
368         END DO
369      END DO
370      !
371      CALL eos_rab( zts, zdep, zab )
372      !
373      !                                   !-------------------!
374      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
375         !                                !-------------------!
376         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
377            DO ji = 1, fs_jpim1   ! vector opt.
378               !                                                   ! i-direction
379               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
380               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
381               !                                                         ! 2*masked bottom density gradient
382               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
383                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
384               !
385               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
386               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
387               !
388               !                                                   ! j-direction
389               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
390               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
391               !                                                         ! 2*masked bottom density gradient
392               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
393                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
394               !
395               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
396               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
397            END DO
398         END DO
399         !
400      ENDIF
401      !
402      !                                   !-------------------!
403      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
404         !                                !-------------------!
405         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
406         !
407         CASE( 1 )                                   != use of upper velocity
408            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
409               DO ji = 1, fs_jpim1   ! vector opt.
410                  !                                                  ! i-direction
411                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
412                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
413                  !                                                          ! 2*masked bottom density gradient
414                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
415                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
416                  !
417                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
418                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
419                  !
420                  !                                                          ! bbl velocity
421                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
422                  !
423                  !                                                  ! j-direction
424                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
425                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
426                  !                                                          ! 2*masked bottom density gradient
427                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
428                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
429                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
430                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
431                  !
432                  !                                                          ! bbl transport
433                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
434               END DO
435            END DO
436            !
437         CASE( 2 )                                 != bbl velocity = F( delta rho )
438            zgbbl = grav * rn_gambbl
439            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
440               DO ji = 1, fs_jpim1   ! vector opt.
441                  !                                                  ! i-direction
442                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
443                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
444                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
445                  !
446                  ikud = mbku_d(ji,jj)
447                  ikus = mbku(ji,jj)
448                  !
449                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
450                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
451                  !                                                          !   masked bottom density gradient
452                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
453                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
454                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
455                  !
456                  !                                                          ! bbl transport (down-slope direction)
457                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
458                  !
459                  !                                                  ! j-direction
460                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
461                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
462                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
463                  !
464                  ikvd = mbkv_d(ji,jj)
465                  ikvs = mbkv(ji,jj)
466                  !
467                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
468                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
469                  !                                                          !   masked bottom density gradient
470                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
471                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
472                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
473                  !
474                  !                                                          ! bbl transport (down-slope direction)
475                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
476               END DO
477            END DO
478         END SELECT
479         !
480      ENDIF
481      !
482      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
483      !
484   END SUBROUTINE bbl
485
486
487   SUBROUTINE tra_bbl_init
488      !!----------------------------------------------------------------------
489      !!                  ***  ROUTINE tra_bbl_init  ***
490      !!
491      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
492      !!
493      !! ** Method  :   Read the nambbl namelist and check the parameters
494      !!              called by nemo_init at the first timestep (kit000)
495      !!----------------------------------------------------------------------
496      INTEGER ::   ji, jj               ! dummy loop indices
497      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
498      INTEGER ::   ios                  !   -      -
499      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
500      !
501      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
502      !!----------------------------------------------------------------------
503      !
504      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
505      !
506      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
507      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
508901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
509
510      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
511      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
512902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
513      IF(lwm) WRITE ( numond, nambbl )
514      !
515      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
516      !
517      IF(lwp) THEN                   !* Parameter control and print
518         WRITE(numout,*)
519         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
520         WRITE(numout,*) '~~~~~~~~~~~~'
521         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
522         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
523         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
524         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
525         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
526      ENDIF
527
528      !                              ! allocate trabbl arrays
529      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
530
531      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
532      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
533
534      !                             !* vertical index of  "deep" bottom u- and v-points
535      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
536         DO ji = 1, jpim1
537            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
538            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
539         END DO
540      END DO
541      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
542      CALL wrk_alloc( jpi, jpj, zmbk )
543      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
544      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
545      CALL wrk_dealloc( jpi, jpj, zmbk )
546
547                                        !* sign of grad(H) at u- and v-points
548      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
549      DO jj = 1, jpjm1
550         DO ji = 1, jpim1
551            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
552            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
553         END DO
554      END DO
555
556      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
557         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
558            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
559            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
560         END DO
561      END DO
562      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
563
564      !                             !* masked diffusive flux coefficients
565      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
566      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
567
568
569      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
570         !
571         SELECT CASE ( jp_cfg )
572         CASE ( 2 )                          ! ORCA_R2
573            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
574            ii0 = 139   ;   ii1 = 140
575            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
576            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
577            !
578            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
579            ii0 = 161   ;   ii1 = 162
580            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
581            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
582            !
583         CASE ( 4 )                          ! ORCA_R4
584            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
585            ii0 =  70   ;   ii1 =  71
586            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
587            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
588         END SELECT
589         !
590      ENDIF
591      !
592      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
593      !
594   END SUBROUTINE tra_bbl_init
595
596#else
597   !!----------------------------------------------------------------------
598   !!   Dummy module :                      No bottom boundary layer scheme
599   !!----------------------------------------------------------------------
600   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
601CONTAINS
602   SUBROUTINE tra_bbl_init               ! Dummy routine
603   END SUBROUTINE tra_bbl_init
604   SUBROUTINE tra_bbl( kt )              ! Dummy routine
605      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
606   END SUBROUTINE tra_bbl
607#endif
608
609   !!======================================================================
610END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.