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

source: NEMO/branches/2018/dev_r9838_ENHANCE04_RK3/src/OCE/TRA/trabbl.F90 @ 9939

Last change on this file since 9939 was 9939, checked in by gm, 6 years ago

#1911 (ENHANCE-04): RK3 branche phased with MLF@9937 branche

  • 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 "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
72   !! $Id$
73   !! Software governed by the CeCILL licence     (./LICENSE)
74   !!----------------------------------------------------------------------
75CONTAINS
76
77   INTEGER FUNCTION tra_bbl_alloc()
78      !!----------------------------------------------------------------------
79      !!                ***  FUNCTION tra_bbl_alloc  ***
80      !!----------------------------------------------------------------------
81      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
82         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) ,     &
83         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                        &
84         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                    STAT=tra_bbl_alloc )
85         !
86      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
87      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
88   END FUNCTION tra_bbl_alloc
89
90
91   SUBROUTINE tra_bbl( kt )
92      !!----------------------------------------------------------------------
93      !!                  ***  ROUTINE bbl  ***
94      !!
95      !! ** Purpose :   Compute the before tracer (t & s) trend associated
96      !!              with the bottom boundary layer and add it to the general
97      !!              trend of tracer equations.
98      !!
99      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
100      !!              diffusive and/or advective contribution to the tracer trend
101      !!              is added to the general tracer trend
102      !!----------------------------------------------------------------------
103      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
104      !
105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   ztrd
106      !!----------------------------------------------------------------------
107      !
108      IF( ln_timing )   CALL timing_start( 'tra_bbl')
109      !
110      IF( l_trdtra )   THEN                         !* Save the T-S input trends
111         ALLOCATE( ztrd(jpi,jpj,jpk,jpts) )
112         ztrd(:,:,:,:) = tsa(:,:,:,:)
113      ENDIF
114
115      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
116
117      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
118         !
119         CALL tra_bbl_dif( tsb, tsa, jpts )
120         IF( ln_ctl )  &
121         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
122            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
123         ! lateral boundary conditions ; just need for outputs
124         CALL lbc_lnk_multi( ahu_bbl, 'U', 1. , ahv_bbl, 'V', 1. )
125         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
126         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
127         !
128      ENDIF
129      !
130      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
131         !
132         CALL tra_bbl_adv( tsb, tsa, jpts )
133         IF(ln_ctl)   &
134         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
135            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
136         ! lateral boundary conditions ; just need for outputs
137         CALL lbc_lnk_multi( utr_bbl, 'U', 1. , vtr_bbl, 'V', 1. )
138         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
139         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
140         !
141      ENDIF
142
143      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
144         ztrd(:,:,:,:) = tsa(:,:,:,:) - ztrd(:,:,:,:)
145         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrd(:,:,:,jp_tem) )
146         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrd(:,:,:,jp_sal) )
147         DEALLOCATE( ztrd )
148      ENDIF
149      !
150      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
151      !
152   END SUBROUTINE tra_bbl
153
154
155   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
156      !!----------------------------------------------------------------------
157      !!                  ***  ROUTINE tra_bbl_dif  ***
158      !!
159      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
160      !!                advection terms.
161      !!
162      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
163      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
164      !!      along bottom slope gradient) an additional lateral 2nd order
165      !!      diffusion along the bottom slope is added to the general
166      !!      tracer trend, otherwise the additional trend is set to 0.
167      !!      A typical value of ahbt is 2000 m2/s (equivalent to
168      !!      a downslope velocity of 20 cm/s if the condition for slope
169      !!      convection is satified)
170      !!
171      !! ** Action  :   pta   increased by the bbl diffusive trend
172      !!
173      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
174      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
175      !!----------------------------------------------------------------------
176      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
177      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
178      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
179      !
180      INTEGER  ::   ji, jj, jn   ! dummy loop indices
181      INTEGER  ::   ik           ! local integers
182      REAL(wp) ::   zbtr         ! local scalars
183      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! workspace
184      !!----------------------------------------------------------------------
185      !
186      DO jn = 1, kjpt                                     ! tracer loop
187         !                                                ! ===========
188         DO jj = 1, jpj
189            DO ji = 1, jpi
190               ik = mbkt(ji,jj)                             ! bottom T-level index
191               zptb(ji,jj) = ptb(ji,jj,ik,jn)               ! bottom before T and S
192            END DO
193         END DO
194         !               
195         DO jj = 2, jpjm1                                    ! Compute the trend
196            DO ji = 2, jpim1
197               ik = mbkt(ji,jj)                            ! bottom T-level index
198               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  &
199                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
200                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
201                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
202                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
203                  &             * r1_e1e2t(ji,jj) / e3t_n(ji,jj,ik)
204            END DO
205         END DO
206         !                                                  ! ===========
207      END DO                                                ! end tracer
208      !                                                     ! ===========
209   END SUBROUTINE tra_bbl_dif
210
211
212   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
213      !!----------------------------------------------------------------------
214      !!                  ***  ROUTINE trc_bbl  ***
215      !!
216      !! ** Purpose :   Compute the before passive tracer trend associated
217      !!     with the bottom boundary layer and add it to the general trend
218      !!     of tracer equations.
219      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
220      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
221      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
222      !!                       transport proportional to the along-slope density gradient
223      !!
224      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
225      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
226      !!----------------------------------------------------------------------
227      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
228      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
229      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
230      !
231      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
232      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
233      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
234      REAL(wp) ::   zbtr, ztra               ! local scalars
235      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
236      !!----------------------------------------------------------------------
237      !
238      !                                                          ! ===========
239      DO jn = 1, kjpt                                            ! tracer loop
240         !                                                       ! ===========
241         DO jj = 1, jpjm1
242            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
243               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
244                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
245                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
246                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
247                  zu_bbl = ABS( utr_bbl(ji,jj) )
248                  !
249                  !                                               ! up  -slope T-point (shelf bottom point)
250                  zbtr = r1_e1e2t(iis,jj) / e3t_n(iis,jj,ikus)
251                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
252                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
253                  !
254                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
255                     zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,jk)
256                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
257                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
258                  END DO
259                  !
260                  zbtr = r1_e1e2t(iid,jj) / e3t_n(iid,jj,ikud)
261                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
262                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
263               ENDIF
264               !
265               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
266                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
267                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
268                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
269                  zv_bbl = ABS( vtr_bbl(ji,jj) )
270                  !
271                  ! up  -slope T-point (shelf bottom point)
272                  zbtr = r1_e1e2t(ji,ijs) / e3t_n(ji,ijs,ikvs)
273                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
274                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
275                  !
276                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
277                     zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,jk)
278                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
279                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
280                  END DO
281                  !                                               ! down-slope T-point (deep bottom point)
282                  zbtr = r1_e1e2t(ji,ijd) / e3t_n(ji,ijd,ikvd)
283                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
284                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
285               ENDIF
286            END DO
287            !
288         END DO
289         !                                                  ! ===========
290      END DO                                                ! end tracer
291      !                                                     ! ===========
292   END SUBROUTINE tra_bbl_adv
293
294
295   SUBROUTINE bbl( kt, kit000, cdtype )
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      !
325      INTEGER  ::   ji, jj                    ! dummy loop indices
326      INTEGER  ::   ik                        ! local integers
327      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
328      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
329      REAL(wp) ::   za, zb, zgdrho            ! local scalars
330      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
331      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
332      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
333      !!----------------------------------------------------------------------
334      !
335      IF( kt == kit000 )  THEN
336         IF(lwp)  WRITE(numout,*)
337         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
338         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
339      ENDIF
340      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
341      DO jj = 1, jpj
342         DO ji = 1, jpi
343            ik = mbkt(ji,jj)                             ! bottom T-level index
344            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
345            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
346            !
347            zdep(ji,jj) = gdept_n(ji,jj,ik)              ! bottom T-level reference depth
348            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
349            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
350         END DO
351      END DO
352      !
353      CALL eos_rab( zts, zdep, zab )
354      !
355      !                                   !-------------------!
356      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
357         !                                !-------------------!
358         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
359            DO ji = 1, fs_jpim1   ! vector opt.
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, -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, -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 DO
380         END DO
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 jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
391               DO ji = 1, fs_jpim1   ! vector opt.
392                  !                                                  ! i-direction
393                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
394                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
395                  !                                                          ! 2*masked bottom density gradient
396                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
397                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
398                  !
399                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
400                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
401                  !
402                  !                                                          ! bbl velocity
403                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
404                  !
405                  !                                                  ! j-direction
406                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
407                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
408                  !                                                          ! 2*masked bottom density gradient
409                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
410                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
411                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
412                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
413                  !
414                  !                                                          ! bbl transport
415                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
416               END DO
417            END DO
418            !
419         CASE( 2 )                                 != bbl velocity = F( delta rho )
420            zgbbl = grav * rn_gambbl
421            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
422               DO ji = 1, fs_jpim1   ! vector opt.
423                  !                                                  ! i-direction
424                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
425                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
426                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
427                  !
428                  ikud = mbku_d(ji,jj)
429                  ikus = mbku(ji,jj)
430                  !
431                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
432                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
433                  !                                                          !   masked bottom density gradient
434                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
435                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
436                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
437                  !
438                  !                                                          ! bbl transport (down-slope direction)
439                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
440                  !
441                  !                                                  ! j-direction
442                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
443                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
444                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
445                  !
446                  ikvd = mbkv_d(ji,jj)
447                  ikvs = mbkv(ji,jj)
448                  !
449                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
450                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
451                  !                                                          !   masked bottom density gradient
452                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
453                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
454                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
455                  !
456                  !                                                          ! bbl transport (down-slope direction)
457                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
458               END DO
459            END DO
460         END SELECT
461         !
462      ENDIF
463      !
464   END SUBROUTINE bbl
465
466
467   SUBROUTINE tra_bbl_init
468      !!----------------------------------------------------------------------
469      !!                  ***  ROUTINE tra_bbl_init  ***
470      !!
471      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
472      !!
473      !! ** Method  :   Read the nambbl namelist and check the parameters
474      !!              called by nemo_init at the first timestep (kit000)
475      !!----------------------------------------------------------------------
476      INTEGER ::   ji, jj                      ! dummy loop indices
477      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
478      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
479      !!
480      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
481      !!----------------------------------------------------------------------
482      !
483      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
484      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
485901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
486      !
487      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
488      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
489902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
490      IF(lwm) WRITE ( numond, nambbl )
491      !
492      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
493      !
494      IF(lwp) THEN                   !* Parameter control and print
495         WRITE(numout,*)
496         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
497         WRITE(numout,*) '~~~~~~~~~~~~'
498         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
499         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
500      ENDIF
501      IF( .NOT.ln_trabbl )   RETURN
502      !
503      IF(lwp) THEN
504         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
505         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
506         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
507         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
508      ENDIF
509      !
510      !                              ! allocate trabbl arrays
511      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
512      !
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      !
516      !                             !* vertical index of  "deep" bottom u- and v-points
517      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
518         DO ji = 1, jpim1
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 DO
522      END DO
523      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
524      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 
525      CALL lbc_lnk_multi( zmbku,'U',1., zmbkv,'V',1.) 
526      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( INT( zmbkv(:,:) ), 1 )
527      !
528      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
529      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
530      DO jj = 1, jpjm1
531         DO ji = 1, jpim1
532            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
533               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
534            ENDIF
535            !
536            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
537               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
538            ENDIF
539         END DO
540      END DO
541      !
542      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
543         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
544            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
545            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
546         END DO
547      END DO
548      CALL lbc_lnk_multi( e3u_bbl_0, 'U', 1. , e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
549      !
550      !                             !* masked diffusive flux coefficients
551      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
552      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
553      !
554   END SUBROUTINE tra_bbl_init
555
556   !!======================================================================
557END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.