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/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/TRA/trabbl.F90 @ 12097

Last change on this file since 12097 was 11536, checked in by smasson, 5 years ago

trunk: merge dev_r10984_HPC-13 into the trunk

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