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 NEMO/trunk/src/OCE/TRA – NEMO

source: NEMO/trunk/src/OCE/TRA/trabbl.F90 @ 14433

Last change on this file since 14433 was 14433, checked in by smasson, 3 years ago

trunk: merge dev_r14312_MPI_Interface into the trunk, #2598

  • Property svn:keywords set to Id
File size: 32.1 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   !!            4.0  ! 2017-04  (G. Madec)  ln_trabbl namelist variable instead of a CPP key
16   !!----------------------------------------------------------------------
17
18   !!----------------------------------------------------------------------
19   !!   tra_bbl_alloc : allocate trabbl arrays
20   !!   tra_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
21   !!   tra_bbl_dif   : generic routine to compute bbl diffusive trend
22   !!   tra_bbl_adv   : generic routine to compute bbl advective trend
23   !!   bbl           : computation of bbl diffu. flux coef. & transport in bottom boundary layer
24   !!   tra_bbl_init  : initialization, namelist read, parameters control
25   !!----------------------------------------------------------------------
26   USE oce            ! ocean dynamics and active tracers
27   USE dom_oce        ! ocean space and time domain
28   USE phycst         ! physical constant
29   USE eosbn2         ! equation of state
30   USE trd_oce        ! trends: ocean variables
31   USE trdtra         ! trends: active tracers
32   !
33   USE iom            ! IOM library
34   USE in_out_manager ! I/O manager
35   USE lbclnk         ! ocean lateral boundary conditions
36   USE prtctl         ! Print control
37   USE timing         ! Timing
38   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   tra_bbl       !  routine called by step.F90
44   PUBLIC   tra_bbl_init  !  routine called by nemogcm.F90
45   PUBLIC   tra_bbl_dif   !  routine called by trcbbl.F90
46   PUBLIC   tra_bbl_adv   !     -      -          -
47   PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
48
49   !                                !!* Namelist nambbl *
50   LOGICAL , PUBLIC ::   ln_trabbl   !: bottom boundary layer flag
51   INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
52   INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
53   !                                            !  =1 : advective bbl using the bottom ocean velocity
54   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
55   REAL(wp), PUBLIC ::   rn_ahtbbl   !: along slope bbl diffusive coefficient [m2/s]
56   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
57
58   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
59
60   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
61   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
62
63   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
64   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
66   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)
67
68   !! * Substitutions
69#  include "do_loop_substitute.h90"
70#  include "domzgr_substitute.h90"
71   !!----------------------------------------------------------------------
72   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
73   !! $Id$
74   !! Software governed by the CeCILL license (see ./LICENSE)
75   !!----------------------------------------------------------------------
76CONTAINS
77
78   INTEGER FUNCTION tra_bbl_alloc()
79      !!----------------------------------------------------------------------
80      !!                ***  FUNCTION tra_bbl_alloc  ***
81      !!----------------------------------------------------------------------
82      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
83         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) ,     &
84         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                        &
85         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                    STAT=tra_bbl_alloc )
86         !
87      CALL mpp_sum ( 'trabbl', tra_bbl_alloc )
88      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
89   END FUNCTION tra_bbl_alloc
90
91
92   SUBROUTINE tra_bbl( kt, Kbb, Kmm, pts, Krhs )
93      !!----------------------------------------------------------------------
94      !!                  ***  ROUTINE bbl  ***
95      !!
96      !! ** Purpose :   Compute the before tracer (t & s) trend associated
97      !!              with the bottom boundary layer and add it to the general
98      !!              trend of tracer equations.
99      !!
100      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
101      !!              diffusive and/or advective contribution to the tracer trend
102      !!              is added to the general tracer trend
103      !!----------------------------------------------------------------------
104      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step
105      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices
106      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
107      !
108      INTEGER  ::   ji, jj, jk   ! Dummy loop indices
109      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
110      !!----------------------------------------------------------------------
111      !
112      IF( ln_timing )   CALL timing_start( 'tra_bbl')
113      !
114      IF( l_trdtra )   THEN                         !* Save the T-S input trends
115         ALLOCATE( ztrdt(jpi,jpj,jpk), ztrds(jpi,jpj,jpk) )
116         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
117         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
118      ENDIF
119
120      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA', Kbb, Kmm )   !* bbl coef. and transport (only if not already done in trcbbl)
121
122      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
123         !
124         CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
125         IF( sn_cfctl%l_prtctl )  &
126         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
127            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
128         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
129            CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
130            CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
131         ENDIF
132         !
133      ENDIF
134      !
135      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
136         !
137         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
138         IF(sn_cfctl%l_prtctl)   &
139         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask, &
140            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
141         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
142            ! lateral boundary conditions ; just need for outputs
143            CALL lbc_lnk( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp )
144            CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
145            CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
146         ENDIF
147         !
148      ENDIF
149
150      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
151         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
152         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
153         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
154         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
155         DEALLOCATE( ztrdt, ztrds )
156      ENDIF
157      !
158      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
159      !
160   END SUBROUTINE tra_bbl
161
162
163   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
164      !!----------------------------------------------------------------------
165      !!                  ***  ROUTINE tra_bbl_dif  ***
166      !!
167      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
168      !!                advection terms.
169      !!
170      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
171      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
172      !!      along bottom slope gradient) an additional lateral 2nd order
173      !!      diffusion along the bottom slope is added to the general
174      !!      tracer trend, otherwise the additional trend is set to 0.
175      !!      A typical value of ahbt is 2000 m2/s (equivalent to
176      !!      a downslope velocity of 20 cm/s if the condition for slope
177      !!      convection is satified)
178      !!
179      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
180      !!
181      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
182      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
183      !!----------------------------------------------------------------------
184      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
185      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
186      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
187      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
188      !
189      INTEGER  ::   ji, jj, jn   ! dummy loop indices
190      INTEGER  ::   ik           ! local integers
191      REAL(wp) ::   zbtr         ! local scalars
192      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zptb   ! workspace
193      !!----------------------------------------------------------------------
194      !
195      DO jn = 1, kjpt                                     ! tracer loop
196         !                                                ! ===========
197         DO_2D( 1, 1, 1, 1 )
198            ik = mbkt(ji,jj)                             ! bottom T-level index
199            zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
200         END_2D
201         !
202         DO_2D( 0, 0, 0, 0 )                               ! Compute the trend
203            ik = mbkt(ji,jj)                            ! bottom T-level index
204            pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
205               &                + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
206               &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
207               &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
208               &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
209               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
210         END_2D
211         !                                                  ! ===========
212      END DO                                                ! end tracer
213      !                                                     ! ===========
214   END SUBROUTINE tra_bbl_dif
215
216
217   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
218      !!----------------------------------------------------------------------
219      !!                  ***  ROUTINE trc_bbl  ***
220      !!
221      !! ** Purpose :   Compute the before passive tracer trend associated
222      !!     with the bottom boundary layer and add it to the general trend
223      !!     of tracer equations.
224      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
225      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
226      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
227      !!                       transport proportional to the along-slope density gradient
228      !!
229      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
230      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
231      !!----------------------------------------------------------------------
232      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
233      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
234      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
235      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
236      !
237      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
238      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
239      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
240      INTEGER  ::   isi, isj                 !   -       -
241      REAL(wp) ::   zbtr, ztra               ! local scalars
242      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
243      !!----------------------------------------------------------------------
244      !
245      IF( ntsi == Nis0 ) THEN ; isi = 1 ; ELSE ; isi = 0 ; ENDIF    ! Avoid double-counting when using tiling
246      IF( ntsj == Njs0 ) THEN ; isj = 1 ; ELSE ; isj = 0 ; ENDIF
247      !                                                          ! ===========
248      DO jn = 1, kjpt                                            ! tracer loop
249         !                                                       ! ===========
250         DO_2D( isi, 0, isj, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west
251            IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
252               ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
253               iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
254               ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
255               zu_bbl = ABS( utr_bbl(ji,jj) )
256               !
257               !                                               ! up  -slope T-point (shelf bottom point)
258               zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
259               ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
260               pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
261               !
262               DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
263                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
264                  ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
265                  pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
266               END DO
267               !
268               zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
269               ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
270               pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
271            ENDIF
272            !
273            IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
274               ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
275               ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
276               ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
277               zv_bbl = ABS( vtr_bbl(ji,jj) )
278               !
279               ! up  -slope T-point (shelf bottom point)
280               zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
281               ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
282               pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
283               !
284               DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
285                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
286                  ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
287                  pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
288               END DO
289               !                                               ! down-slope T-point (deep bottom point)
290               zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
291               ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
292               pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
293            ENDIF
294         END_2D
295         !                                                       ! ===========
296      END DO                                                     ! end tracer
297      !                                                          ! ===========
298   END SUBROUTINE tra_bbl_adv
299
300
301   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
302      !!----------------------------------------------------------------------
303      !!                  ***  ROUTINE bbl  ***
304      !!
305      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
306      !!                advection terms.
307      !!
308      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
309      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
310      !!      along bottom slope gradient) an additional lateral 2nd order
311      !!      diffusion along the bottom slope is added to the general
312      !!      tracer trend, otherwise the additional trend is set to 0.
313      !!      A typical value of ahbt is 2000 m2/s (equivalent to
314      !!      a downslope velocity of 20 cm/s if the condition for slope
315      !!      convection is satified)
316      !!              * advective bbl (nn_bbl_adv=1 or 2) :
317      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
318      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
319      !!        i.e. transport proportional to the along-slope density gradient
320      !!
321      !!      NB: the along slope density gradient is evaluated using the
322      !!      local density (i.e. referenced at a common local depth).
323      !!
324      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
325      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
326      !!----------------------------------------------------------------------
327      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
328      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
329      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
330      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
331      !
332      INTEGER  ::   ji, jj                    ! dummy loop indices
333      INTEGER  ::   ik                        ! local integers
334      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
335      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
336      REAL(wp) ::   za, zb, zgdrho            ! local scalars
337      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
338      REAL(wp), DIMENSION(A2D(nn_hls),jpts)   :: zts, zab         ! 3D workspace
339      REAL(wp), DIMENSION(A2D(nn_hls))        :: zub, zvb, zdep   ! 2D workspace
340      !!----------------------------------------------------------------------
341      !
342      IF( ntile == 0 .OR. ntile == 1 )  THEN                       ! Do only on the first tile
343         IF( kt == kit000 )  THEN
344            IF(lwp)  WRITE(numout,*)
345            IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
346            IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
347         ENDIF
348      ENDIF
349      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
350      DO_2D( 1, 1, 1, 1 )
351         ik = mbkt(ji,jj)                             ! bottom T-level index
352         zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
353         zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
354         !
355         zdep(ji,jj) = gdept(ji,jj,ik,Kmm)            ! bottom T-level reference depth
356         zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)      ! bottom velocity
357         zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
358      END_2D
359      !
360      CALL eos_rab( zts, zdep, zab, Kmm )
361      !
362      !                                   !-------------------!
363      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
364         !                                !-------------------!
365         DO_2D( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
366            !                                                   ! i-direction
367            za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
368            zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
369            !                                                         ! 2*masked bottom density gradient
370            zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
371               &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
372            !
373            zsign  = SIGN(  0.5_wp, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
374            ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
375            !
376            !                                                   ! j-direction
377            za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
378            zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
379            !                                                         ! 2*masked bottom density gradient
380            zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
381               &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
382            !
383            zsign = SIGN(  0.5_wp, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
384            ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
385         END_2D
386         !
387      ENDIF
388      !
389      !                                   !-------------------!
390      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
391         !                                !-------------------!
392         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
393         !
394         CASE( 1 )                                   != use of upper velocity
395            DO_2D( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
396               !                                                  ! i-direction
397               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
398               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
399               !                                                          ! 2*masked bottom density gradient
400               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
401                         - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
402               !
403               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
404               zsigna= SIGN(  0.5_wp, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
405               !
406               !                                                          ! bbl velocity
407               utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
408               !
409               !                                                  ! j-direction
410               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
411               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
412               !                                                          ! 2*masked bottom density gradient
413               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
414                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
415               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
416               zsigna= SIGN(  0.5_wp, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
417               !
418               !                                                          ! bbl transport
419               vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
420            END_2D
421            !
422         CASE( 2 )                                 != bbl velocity = F( delta rho )
423            zgbbl = grav * rn_gambbl
424            DO_2D( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down
425               !                                                  ! i-direction
426               ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
427               iid  = ji + MAX( 0, mgrhu(ji,jj) )
428               iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
429               !
430               ikud = mbku_d(ji,jj)
431               ikus = mbku(ji,jj)
432               !
433               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
434               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
435               !                                                          !   masked bottom density gradient
436               zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
437                  &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
438               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
439               !
440               !                                                          ! bbl transport (down-slope direction)
441               utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
442               !
443               !                                                  ! j-direction
444               !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
445               ijd  = jj + MAX( 0, mgrhv(ji,jj) )
446               ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
447               !
448               ikvd = mbkv_d(ji,jj)
449               ikvs = mbkv(ji,jj)
450               !
451               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
452               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
453               !                                                          !   masked bottom density gradient
454               zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
455                  &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
456               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
457               !
458               !                                                          ! bbl transport (down-slope direction)
459               vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
460            END_2D
461         END SELECT
462         !
463      ENDIF
464      !
465   END SUBROUTINE bbl
466
467
468   SUBROUTINE tra_bbl_init
469      !!----------------------------------------------------------------------
470      !!                  ***  ROUTINE tra_bbl_init  ***
471      !!
472      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
473      !!
474      !! ** Method  :   Read the nambbl namelist and check the parameters
475      !!              called by nemo_init at the first timestep (kit000)
476      !!----------------------------------------------------------------------
477      INTEGER ::   ji, jj                      ! dummy loop indices
478      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
479      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
480      !!
481      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
482      !!----------------------------------------------------------------------
483      !
484      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
485901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist' )
486      !
487      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
488902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist' )
489      IF(lwm) WRITE ( numond, nambbl )
490      !
491      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
492      !
493      IF(lwp) THEN                   !* Parameter control and print
494         WRITE(numout,*)
495         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
496         WRITE(numout,*) '~~~~~~~~~~~~'
497         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
498         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
499      ENDIF
500      IF( .NOT.ln_trabbl )   RETURN
501      !
502      IF(lwp) THEN
503         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
504         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
505         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
506         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
507      ENDIF
508      !
509      !                              ! allocate trabbl arrays
510      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
511      !
512      IF(lwp) THEN
513         IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
514         IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
515      ENDIF
516      !
517      !                             !* vertical index of  "deep" bottom u- and v-points
518      DO_2D( 1, 0, 1, 0 )                 ! (the "shelf" bottom k-indices are mbku and mbkv)
519         mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
520         mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
521      END_2D
522      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
523      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp )
524      CALL lbc_lnk( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp)
525      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
526      !
527      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
528      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
529      DO_2D( 1, 0, 1, 0 )
530         IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
531            mgrhu(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
532         ENDIF
533         !
534         IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
535            mgrhv(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
536         ENDIF
537      END_2D
538      !
539      DO_2D( 1, 0, 1, 0 )           !* bbl thickness at u- (v-) point; minimum of top & bottom e3u_0 (e3v_0)
540         e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
541         e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
542      END_2D
543      CALL lbc_lnk( 'trabbl', e3u_bbl_0, 'U', 1.0_wp , e3v_bbl_0, 'V', 1.0_wp )      ! lateral boundary conditions
544      !
545      !                             !* masked diffusive flux coefficients
546      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
547      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
548      !
549   END SUBROUTINE tra_bbl_init
550
551   !!======================================================================
552END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.