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/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13296_HPC-07_mocavero_mpi3/src/OCE/TRA/trabbl.F90 @ 13630

Last change on this file since 13630 was 13630, checked in by mocavero, 4 years ago

Add neighborhood collectives calls in the NEMO src - ticket #2496

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