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/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/UKMO/dev_r5518_STOPACK/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 11384

Last change on this file since 11384 was 11384, checked in by mattmartin, 5 years ago

Included Andrea Storto's STOPACK code into NEMO3.6 branch.

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