source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/trabbl.F90 @ 13552

Last change on this file since 13552 was 13552, checked in by hadcv, 5 months ago

#2365: Remove diagnostic output lbc_lnk calls

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