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_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 8215

Last change on this file since 8215 was 8215, checked in by gm, 7 years ago

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