New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trabbl.F90 in branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 8882

Last change on this file since 8882 was 8882, checked in by flavoni, 6 years ago

dev_CNRS_2017 branch: merged dev_r7881_ENHANCE09_RK3 with trunk r8864

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