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

source: branches/2013/dev_r3858_CNRS3_Ediag/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 3894

Last change on this file since 3894 was 3894, checked in by gm, 11 years ago

dev_r3858_CNRS3_Ediag: #927 simplification of eosbn2, use of alpha & beta in zdfddm, trabbl...

  • Property svn:keywords set to Id
File size: 34.3 KB
Line 
1MODULE trabbl
2   !!==============================================================================
3   !!                       ***  MODULE  trabbl  ***
4   !! Ocean physics :  advective and/or diffusive bottom boundary layer scheme
5   !!==============================================================================
6   !! History :  OPA  ! 1996-06  (L. Mortier)  Original code
7   !!            8.0  ! 1997-11  (G. Madec)    Optimization
8   !!   NEMO     1.0  ! 2002-08  (G. Madec)  free form + modules
9   !!             -   ! 2004-01  (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl
10   !!            3.3  ! 2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization
11   !!             -   ! 2010-04  (G. Madec)  Campin & Goosse advective bbl
12   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
13   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
14   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
15   !!----------------------------------------------------------------------
16#if   defined key_trabbl   ||   defined key_esopa
17   !!----------------------------------------------------------------------
18   !!   'key_trabbl'   or                             bottom boundary layer
19   !!----------------------------------------------------------------------
20   !!   tra_bbl_alloc : allocate trabbl arrays
21   !!   tra_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
22   !!   tra_bbl_dif   : generic routine to compute bbl diffusive trend
23   !!   tra_bbl_adv   : generic routine to compute bbl advective trend
24   !!   bbl           : computation of bbl diffu. flux coef. & transport in bottom boundary layer
25   !!   tra_bbl_init  : initialization, namelist read, parameters control
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and active tracers
28   USE dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constant
30   USE eosbn2         ! equation of state
31   USE trd_oce        ! trends: ocean variables
32   USE trdtra         ! trends manager: tracers
33   !
34   USE iom            ! IOM server               
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 = 0           !: =1   : diffusive bbl or not (=0)
55   INTEGER , PUBLIC ::   nn_bbl_adv = 0           !: =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  = 1.e3_wp     !: along slope bbl diffusive coefficient [m2/s]
59   REAL(wp), PUBLIC ::   rn_gambbl  = 10.0_wp     !: 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   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   r1_e1e2t               ! inverse of the cell surface at t-point      [1/m2] (PUBLIC for TAM)
71
72   !! * Substitutions
73#  include "domzgr_substitute.h90"
74#  include "vectopt_loop_substitute.h90"
75   !!----------------------------------------------------------------------
76   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
77   !! $Id$
78   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
79   !!----------------------------------------------------------------------
80CONTAINS
81
82   INTEGER FUNCTION tra_bbl_alloc()
83      !!----------------------------------------------------------------------
84      !!                ***  FUNCTION tra_bbl_alloc  ***
85      !!----------------------------------------------------------------------
86      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
87         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
88         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
89         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) , r1_e1e2t(jpi,jpj)                  , STAT= tra_bbl_alloc )
90         !
91      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
92      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
93   END FUNCTION tra_bbl_alloc
94
95
96   SUBROUTINE tra_bbl( kt )
97      !!----------------------------------------------------------------------
98      !!                  ***  ROUTINE bbl  ***
99      !!
100      !! ** Purpose :   Compute the before tracer (t & s) trend associated
101      !!              with the bottom boundary layer and add it to the general
102      !!              trend of tracer equations.
103      !!
104      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
105      !!              diffusive and/or advective contribution to the tracer trend
106      !!              is added to the general tracer trend
107      !!----------------------------------------------------------------------
108      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
109      !
110      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
111      !!----------------------------------------------------------------------
112      !
113      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
114      !
115      IF( l_trdtra )   THEN                         !* Save ta and sa trends
116         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
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         CALL wrk_dealloc( jpi, jpj, jpk, 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), POINTER, DIMENSION(:,:) :: zptb
191      !!----------------------------------------------------------------------
192      !
193      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
194      !
195      CALL wrk_alloc( jpi, jpj, zptb )
196      !
197      DO jn = 1, kjpt                                     ! tracer loop
198         !                                                ! ===========
199         DO jj = 1, jpj
200            DO ji = 1, jpi
201               zptb(ji,jj) = ptb(ji,jj,mbkt(ji,jj),jn)       ! bottom before T and S
202            END DO
203         END DO
204         !               
205         DO jj = 2, jpjm1                                    ! Compute the trend
206            DO ji = 2, jpim1
207               ik = mbkt(ji,jj)                              ! bottom T-level index
208               zbtr = r1_e1e2t(ji,jj)  / fse3t(ji,jj,ik)
209               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
210                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
211                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
212                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
213                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
214            END DO
215         END DO
216         !                                                  ! ===========
217      END DO                                                ! end tracer
218      !                                                     ! ===========
219      CALL wrk_dealloc( jpi, jpj, zptb )
220      !
221      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
222      !
223   END SUBROUTINE tra_bbl_dif
224
225
226   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
227      !!----------------------------------------------------------------------
228      !!                  ***  ROUTINE trc_bbl  ***
229      !!
230      !! ** Purpose :   Compute the before passive tracer trend associated
231      !!     with the bottom boundary layer and add it to the general trend
232      !!     of tracer equations.
233      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
234      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
235      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
236      !!                       transport proportional to the along-slope density gradient
237      !!
238      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
239      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
240      !!----------------------------------------------------------------------
241      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
242      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
243      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
244      !
245      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
246      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
247      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
248      REAL(wp) ::   zbtr, ztra               ! local scalars
249      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
250      !!----------------------------------------------------------------------
251      !
252      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
253      !                                                          ! ===========
254      DO jn = 1, kjpt                                            ! tracer loop
255         !                                                       ! ===========
256         DO jj = 1, jpjm1
257            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
258               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
259                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
260                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
261                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
262                  zu_bbl = ABS( utr_bbl(ji,jj) )
263                  !
264                  !                                               ! up  -slope T-point (shelf bottom point)
265                  zbtr = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus)
266                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
267                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
268                  !
269                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
270                     zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk)
271                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
272                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
273                  END DO
274                  !
275                  zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud)
276                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
277                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
278               ENDIF
279               !
280               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
281                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
282                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
283                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
284                  zv_bbl = ABS( vtr_bbl(ji,jj) )
285                  !
286                  ! up  -slope T-point (shelf bottom point)
287                  zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs)
288                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
289                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
290                  !
291                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
292                     zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk)
293                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
294                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
295                  END DO
296                  !                                               ! down-slope T-point (deep bottom point)
297                  zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd)
298                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
299                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
300               ENDIF
301            END DO
302            !
303         END DO
304         !                                                  ! ===========
305      END DO                                                ! end tracer
306      !                                                     ! ===========
307      !
308      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
309      !
310   END SUBROUTINE tra_bbl_adv
311
312
313   SUBROUTINE bbl( kt, kit000, cdtype )
314      !!----------------------------------------------------------------------
315      !!                  ***  ROUTINE bbl  ***
316      !!
317      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
318      !!                advection terms.
319      !!
320      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
321      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
322      !!      along bottom slope gradient) an additional lateral 2nd order
323      !!      diffusion along the bottom slope is added to the general
324      !!      tracer trend, otherwise the additional trend is set to 0.
325      !!      A typical value of ahbt is 2000 m2/s (equivalent to
326      !!      a downslope velocity of 20 cm/s if the condition for slope
327      !!      convection is satified)
328      !!              * advective bbl (nn_bbl_adv=1 or 2) :
329      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
330      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
331      !!        i.e. transport proportional to the along-slope density gradient
332      !!
333      !!      NB: the along slope density gradient is evaluated using the
334      !!      local density (i.e. referenced at a common local depth).
335      !!
336      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
337      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
338      !!----------------------------------------------------------------------
339      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
340      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
341      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
342      !
343      INTEGER  ::   ji, jj                      ! dummy loop indices
344      INTEGER  ::   ik, ikip1, ikjp1            ! local integers
345      INTEGER  ::   iis, iid, ikus, ikud        !   -       -
346      INTEGER  ::   ijs, ijd, ikvs, ikvd        !   -       -
347      REAL(wp) ::   za, zt, zsign , zgbbl       ! local scalars
348      REAL(wp) ::   zb, zs, zsigna, zgdrho, zh  !   -      -
349      REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts, zab   ! 3D workspace
350      REAL(wp), DIMENSION(jpi,jpj)      :: zub, zvb   ! 2D workspace
351      !!----------------------------------------------------------------------
352      !
353      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
354      !
355      IF( kt == kit000 )  THEN
356         IF(lwp)  WRITE(numout,*)
357         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
358         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
359      ENDIF
360      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
361      DO jj = 1, jpj
362         DO ji = 1, jpi
363            ik = mbkt(ji,jj)                             ! bottom T-level index
364            zts(ji,jj,jp_tem) = tsb  (ji,jj,ik,jp_tem)   ! bottom before T and S
365            zts(ji,jj,jp_sal) = tsb  (ji,jj,ik,jp_sal)
366            zab(ji,jj,jp_tem) = rab_b(ji,jj,ik,jp_tem)   ! bottom alpha and beta
367            zab(ji,jj,jp_sal) = rab_b(ji,jj,ik,jp_sal)
368            !
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      !                                   !-------------------!
375      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
376         !                                !-------------------!
377         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
378            DO ji = 1, fs_jpim1   ! vector opt.
379               !                                                   ! i-direction
380               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
381               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
382               !                                                         ! 2*masked bottom density gradient
383               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
384                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
385               !
386               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
387               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
388               !
389               !                                                   ! j-direction
390               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
391               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
392               !                                                         ! 2*masked bottom density gradient
393               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
394                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
395               !
396               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
397               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
398            END DO
399         END DO
400         !
401      ENDIF
402
403      !                                   !-------------------!
404      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
405         !                                !-------------------!
406         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
407         !
408         CASE( 1 )                                   != use of upper velocity
409            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
410               DO ji = 1, fs_jpim1   ! vector opt.
411                  !                                                  ! i-direction
412                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
413                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
414                  !                                                          ! 2*masked bottom density gradient
415                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
416                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
417                  !
418                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
419                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
420                  !
421                  !                                                          ! bbl velocity
422                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
423                  !
424                  !                                                  ! j-direction
425                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
426                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
427                  !                                                          ! 2*masked bottom density gradient
428                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
429                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
430                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
431                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
432                  !
433                  !                                                          ! bbl transport
434                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
435               END DO
436            END DO
437            !
438         CASE( 2 )                                 != bbl velocity = F( delta rho )
439            zgbbl = grav * rn_gambbl
440            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
441               DO ji = 1, fs_jpim1   ! vector opt.
442                  !                                                  ! i-direction
443                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
444                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
445                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
446                  !
447                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
448                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
449                  !                                                          !   masked bottom density gradient
450                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
451                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
452                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
453                  !
454                  !                                                          ! bbl transport (down-slope direction)
455                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
456                  !
457                  !                                                  ! j-direction
458                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
459                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
460                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
461                  !
462                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
463                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
464                  !                                                          !   masked bottom density gradient
465                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
466                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
467                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
468                  !
469                  !                                                          ! bbl transport (down-slope direction)
470                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
471               END DO
472            END DO
473         END SELECT
474         !
475      ENDIF
476      !
477      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
478      !
479   END SUBROUTINE bbl
480
481
482   SUBROUTINE tra_bbl_init
483      !!----------------------------------------------------------------------
484      !!                  ***  ROUTINE tra_bbl_init  ***
485      !!
486      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
487      !!
488      !! ** Method  :   Read the nambbl namelist and check the parameters
489      !!              called by nemo_init at the first timestep (kit000)
490      !!----------------------------------------------------------------------
491      INTEGER ::   ji, jj               ! dummy loop indices
492      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
493      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
494      !!
495      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
496      !!----------------------------------------------------------------------
497      !
498      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
499      !
500      CALL wrk_alloc( jpi, jpj, zmbk )
501      !
502
503      REWIND ( numnam )              !* Read Namelist nambbl : bottom boundary layer scheme
504      READ   ( numnam, nambbl )
505      !
506      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
507      !
508      IF(lwp) THEN                   !* Parameter control and print
509         WRITE(numout,*)
510         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
511         WRITE(numout,*) '~~~~~~~~~~~~'
512         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
513         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
514         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
515         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
516         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
517      ENDIF
518
519      !                              ! allocate trabbl arrays
520      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
521
522      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
523      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
524
525      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
526
527
528      !                             !* inverse of surface of T-cells
529      r1_e1e2t(:,:) = 1._wp / ( e1t(:,:) * e2t(:,:) )
530
531      !                             !* vertical index of  "deep" bottom u- and v-points
532      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
533         DO ji = 1, jpim1
534            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
535            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
536         END DO
537      END DO
538      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
539      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
540      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
541
542                                        !* sign of grad(H) at u- and v-points
543      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
544      DO jj = 1, jpjm1
545         DO ji = 1, jpim1
546            mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
547            mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
548         END DO
549      END DO
550
551      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
552         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
553            e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) )
554            e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) )
555         END DO
556      END DO
557      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
558
559      !                             !* masked diffusive flux coefficients
560      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
561      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
562
563
564      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
565         !
566         SELECT CASE ( jp_cfg )
567         CASE ( 2 )                          ! ORCA_R2
568            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
569            ii0 = 139   ;   ii1 = 140
570            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
571            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
572            !
573            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
574            ii0 = 161   ;   ii1 = 162
575            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
576            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
577            !
578         CASE ( 4 )                          ! ORCA_R4
579            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
580            ii0 =  70   ;   ii1 =  71
581            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
582            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
583         END SELECT
584         !
585      ENDIF
586      !
587      CALL wrk_dealloc( jpi, jpj, zmbk )
588      !
589      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
590      !
591   END SUBROUTINE tra_bbl_init
592
593#else
594   !!----------------------------------------------------------------------
595   !!   Dummy module :                      No bottom boundary layer scheme
596   !!----------------------------------------------------------------------
597   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
598CONTAINS
599   SUBROUTINE tra_bbl_init               ! Dummy routine
600   END SUBROUTINE tra_bbl_init
601   SUBROUTINE tra_bbl( kt )              ! Dummy routine
602      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
603   END SUBROUTINE tra_bbl
604#endif
605
606   !!======================================================================
607END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.