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

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 8521

Last change on this file since 8521 was 8521, checked in by marc, 7 years ago

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