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/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

source: NEMO/branches/2021/dev_r14273_HPC-02_Daley_Tiling/src/OCE/TRA/trabbl.F90 @ 14636

Last change on this file since 14636 was 14636, checked in by hadcv, 3 years ago

#2600: Merge in dev_r14393_HPC-03_Mele_Comm_Cleanup [14538:14609]

  • Property svn:keywords set to Id
File size: 31.8 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         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
129         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
130         !
131      ENDIF
132      !
133      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
134         !
135         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
136         IF(sn_cfctl%l_prtctl)   &
137         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask, &
138            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
139         IF( .NOT. l_istiled .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
140            ! lateral boundary conditions ; just need for outputs
141            ! [ comm_cleanup ] ! no need lbc_lnk for outputs
142            ! CALL lbc_lnk( 'trabbl', utr_bbl, 'U', 1.0_wp , vtr_bbl, 'V', 1.0_wp )
143         ENDIF
144         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
145         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
146         !
147      ENDIF
148
149      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
150         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
151         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
152         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
153         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
154         DEALLOCATE( ztrdt, ztrds )
155      ENDIF
156      !
157      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
158      !
159   END SUBROUTINE tra_bbl
160
161
162   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
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  :   pt_rhs   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   ) ::   pt     ! before and now tracer fields
185      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
186      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
187      !
188      INTEGER  ::   ji, jj, jn   ! dummy loop indices
189      INTEGER  ::   ik           ! local integers
190      REAL(wp) ::   zbtr         ! local scalars
191      REAL(wp), DIMENSION(A2D(nn_hls)) ::   zptb   ! workspace
192      !!----------------------------------------------------------------------
193      !
194      DO jn = 1, kjpt                                     ! tracer loop
195         !                                                ! ===========
196         DO_2D( 1, 1, 1, 1 )
197            ik = mbkt(ji,jj)                             ! bottom T-level index
198            zptb(ji,jj) = pt(ji,jj,ik,jn)                ! bottom before T and S
199         END_2D
200         !
201         DO_2D( 0, 0, 0, 0 )                               ! Compute the trend
202            ik = mbkt(ji,jj)                            ! bottom T-level index
203            pt_rhs(ji,jj,ik,jn) = pt_rhs(ji,jj,ik,jn)                                                  &
204               &                + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
205               &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
206               &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
207               &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
208               &                * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
209         END_2D
210         !                                                  ! ===========
211      END DO                                                ! end tracer
212      !                                                     ! ===========
213   END SUBROUTINE tra_bbl_dif
214
215
216   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
217      !!----------------------------------------------------------------------
218      !!                  ***  ROUTINE trc_bbl  ***
219      !!
220      !! ** Purpose :   Compute the before passive tracer trend associated
221      !!     with the bottom boundary layer and add it to the general trend
222      !!     of tracer equations.
223      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
224      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
225      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
226      !!                       transport proportional to the along-slope density gradient
227      !!
228      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
229      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
230      !!----------------------------------------------------------------------
231      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
232      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
233      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
234      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
235      !
236      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
237      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
238      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
239      REAL(wp) ::   zbtr, ztra               ! local scalars
240      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
241      !!----------------------------------------------------------------------
242      !                                                          ! ===========
243      DO jn = 1, kjpt                                            ! tracer loop
244         !                                                       ! ===========
245         DO_2D_OVR( 1, 0, 1, 0 )            ! CAUTION start from i=1 to update i=2 when cyclic east-west
246            IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
247               ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
248               iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
249               ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
250               zu_bbl = ABS( utr_bbl(ji,jj) )
251               !
252               !                                               ! up  -slope T-point (shelf bottom point)
253               zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
254               ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
255               pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
256               !
257               DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
258                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
259                  ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
260                  pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
261               END DO
262               !
263               zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
264               ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
265               pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
266            ENDIF
267            !
268            IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
269               ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
270               ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
271               ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
272               zv_bbl = ABS( vtr_bbl(ji,jj) )
273               !
274               ! up  -slope T-point (shelf bottom point)
275               zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
276               ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
277               pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
278               !
279               DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
280                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
281                  ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
282                  pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
283               END DO
284               !                                               ! down-slope T-point (deep bottom point)
285               zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
286               ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
287               pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
288            ENDIF
289         END_2D
290         !                                                       ! ===========
291      END DO                                                     ! end tracer
292      !                                                          ! ===========
293   END SUBROUTINE tra_bbl_adv
294
295
296   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
297      !!----------------------------------------------------------------------
298      !!                  ***  ROUTINE bbl  ***
299      !!
300      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
301      !!                advection terms.
302      !!
303      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
304      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
305      !!      along bottom slope gradient) an additional lateral 2nd order
306      !!      diffusion along the bottom slope is added to the general
307      !!      tracer trend, otherwise the additional trend is set to 0.
308      !!      A typical value of ahbt is 2000 m2/s (equivalent to
309      !!      a downslope velocity of 20 cm/s if the condition for slope
310      !!      convection is satified)
311      !!              * advective bbl (nn_bbl_adv=1 or 2) :
312      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
313      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
314      !!        i.e. transport proportional to the along-slope density gradient
315      !!
316      !!      NB: the along slope density gradient is evaluated using the
317      !!      local density (i.e. referenced at a common local depth).
318      !!
319      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
320      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
321      !!----------------------------------------------------------------------
322      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
323      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
324      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
325      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
326      !
327      INTEGER  ::   ji, jj                    ! dummy loop indices
328      INTEGER  ::   ik                        ! local integers
329      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
330      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
331      REAL(wp) ::   za, zb, zgdrho            ! local scalars
332      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
333      REAL(wp), DIMENSION(A2D(nn_hls),jpts)   :: zts, zab         ! 3D workspace
334      REAL(wp), DIMENSION(A2D(nn_hls))        :: zub, zvb, zdep   ! 2D workspace
335      !!----------------------------------------------------------------------
336      !
337      IF( .NOT. l_istiled .OR. ntile == 1 )  THEN                       ! Do only on the first tile
338         IF( kt == kit000 )  THEN
339            IF(lwp)  WRITE(numout,*)
340            IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
341            IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
342         ENDIF
343      ENDIF
344      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
345      DO_2D( 1, 1, 1, 1 )
346         ik = mbkt(ji,jj)                             ! bottom T-level index
347         zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
348         zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
349         !
350         zdep(ji,jj) = gdept(ji,jj,ik,Kmm)            ! bottom T-level reference depth
351         zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)      ! bottom velocity
352         zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
353      END_2D
354      !
355      CALL eos_rab( zts, zdep, zab, Kmm )
356      !
357      !                                   !-------------------!
358      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
359         !                                !-------------------!
360         DO_2D( 1, 0, 1, 0 )                   ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
361            !                                                   ! i-direction
362            za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
363            zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
364            !                                                         ! 2*masked bottom density gradient
365            zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
366               &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
367            !
368            zsign  = SIGN(  0.5_wp, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
369            ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
370            !
371            !                                                   ! j-direction
372            za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
373            zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
374            !                                                         ! 2*masked bottom density gradient
375            zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
376               &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
377            !
378            zsign = SIGN(  0.5_wp, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
379            ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
380         END_2D
381         !
382      ENDIF
383      !
384      !                                   !-------------------!
385      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
386         !                                !-------------------!
387         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
388         !
389         CASE( 1 )                                   != use of upper velocity
390            DO_2D( 1, 0, 1, 0 )                              ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
391               !                                                  ! i-direction
392               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
393               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
394               !                                                          ! 2*masked bottom density gradient
395               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
396                         - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
397               !
398               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
399               zsigna= SIGN(  0.5_wp, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
400               !
401               !                                                          ! bbl velocity
402               utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
403               !
404               !                                                  ! j-direction
405               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
406               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
407               !                                                          ! 2*masked bottom density gradient
408               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
409                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
410               zsign = SIGN(  0.5_wp, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
411               zsigna= SIGN(  0.5_wp, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
412               !
413               !                                                          ! bbl transport
414               vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
415            END_2D
416            !
417         CASE( 2 )                                 != bbl velocity = F( delta rho )
418            zgbbl = grav * rn_gambbl
419            DO_2D( 1, 0, 1, 0 )                         ! criteria: rho_up > rho_down
420               !                                                  ! i-direction
421               ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
422               iid  = ji + MAX( 0, mgrhu(ji,jj) )
423               iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
424               !
425               ikud = mbku_d(ji,jj)
426               ikus = mbku(ji,jj)
427               !
428               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
429               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
430               !                                                          !   masked bottom density gradient
431               zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
432                  &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
433               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
434               !
435               !                                                          ! bbl transport (down-slope direction)
436               utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
437               !
438               !                                                  ! j-direction
439               !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
440               ijd  = jj + MAX( 0, mgrhv(ji,jj) )
441               ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
442               !
443               ikvd = mbkv_d(ji,jj)
444               ikvs = mbkv(ji,jj)
445               !
446               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
447               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
448               !                                                          !   masked bottom density gradient
449               zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
450                  &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
451               zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
452               !
453               !                                                          ! bbl transport (down-slope direction)
454               vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
455            END_2D
456         END SELECT
457         !
458      ENDIF
459      !
460   END SUBROUTINE bbl
461
462
463   SUBROUTINE tra_bbl_init
464      !!----------------------------------------------------------------------
465      !!                  ***  ROUTINE tra_bbl_init  ***
466      !!
467      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
468      !!
469      !! ** Method  :   Read the nambbl namelist and check the parameters
470      !!              called by nemo_init at the first timestep (kit000)
471      !!----------------------------------------------------------------------
472      INTEGER ::   ji, jj                      ! dummy loop indices
473      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
474      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
475      !!
476      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
477      !!----------------------------------------------------------------------
478      !
479      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
480901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist' )
481      !
482      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
483902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist' )
484      IF(lwm) WRITE ( numond, nambbl )
485      !
486      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
487      !
488      IF(lwp) THEN                   !* Parameter control and print
489         WRITE(numout,*)
490         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
491         WRITE(numout,*) '~~~~~~~~~~~~'
492         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
493         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
494      ENDIF
495      IF( .NOT.ln_trabbl )   RETURN
496      !
497      IF(lwp) THEN
498         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
499         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
500         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
501         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
502      ENDIF
503      !
504      !                              ! allocate trabbl arrays
505      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
506      !
507      IF(lwp) THEN
508         IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
509         IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
510      ENDIF
511      !
512      !                             !* vertical index of  "deep" bottom u- and v-points
513      DO_2D( 1, 0, 1, 0 )                 ! (the "shelf" bottom k-indices are mbku and mbkv)
514         mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
515         mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
516      END_2D
517      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
518      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp )
519      CALL lbc_lnk( 'trabbl', zmbku,'U',1.0_wp, zmbkv,'V',1.0_wp)
520      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
521      !
522      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
523      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
524      DO_2D( 1, 0, 1, 0 )
525         IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
526            mgrhu(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
527         ENDIF
528         !
529         IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
530            mgrhv(ji,jj) = INT(  SIGN( 1.0_wp, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
531         ENDIF
532      END_2D
533      !
534      DO_2D( 1, 0, 1, 0 )           !* bbl thickness at u- (v-) point; minimum of top & bottom e3u_0 (e3v_0)
535         e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
536         e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
537      END_2D
538      CALL lbc_lnk( 'trabbl', e3u_bbl_0, 'U', 1.0_wp , e3v_bbl_0, 'V', 1.0_wp )      ! lateral boundary conditions
539      !
540      !                             !* masked diffusive flux coefficients
541      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
542      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
543      !
544   END SUBROUTINE tra_bbl_init
545
546   !!======================================================================
547END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.