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 @ 14607

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

#2600: Changes for XIOS development

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