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

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

source: branches/UKMO/dev_1d_bugfixes_tocommit/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 13191

Last change on this file since 13191 was 13191, checked in by jwhile, 4 years ago

Updates for 1d runnig

File size: 36.7 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   !!----------------------------------------------------------------------
16#if   defined key_trabbl   ||   defined key_esopa
17   !!----------------------------------------------------------------------
18   !!   'key_trabbl'   or                             bottom boundary layer
19   !!----------------------------------------------------------------------
20   !!   tra_bbl_alloc : allocate trabbl arrays
21   !!   tra_bbl       : update the tracer trends due to the bottom boundary layer (advective and/or diffusive)
22   !!   tra_bbl_dif   : generic routine to compute bbl diffusive trend
23   !!   tra_bbl_adv   : generic routine to compute bbl advective trend
24   !!   bbl           : computation of bbl diffu. flux coef. & transport in bottom boundary layer
25   !!   tra_bbl_init  : initialization, namelist read, parameters control
26   !!----------------------------------------------------------------------
27   USE oce            ! ocean dynamics and active tracers
28   USE dom_oce        ! ocean space and time domain
29   USE phycst         ! physical constant
30   USE eosbn2         ! equation of state
31   USE trd_oce     ! trends: ocean variables
32   USE trdtra         ! trends: active tracers
33   !
34   USE iom            ! IOM library
35   USE in_out_manager ! I/O manager
36   USE lbclnk         ! ocean lateral boundary conditions
37   USE prtctl         ! Print control
38   USE wrk_nemo       ! Memory Allocation
39   USE timing         ! Timing
40   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
41   USE stopack
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC   tra_bbl       !  routine called by step.F90
47   PUBLIC   tra_bbl_init  !  routine called by opa.F90
48   PUBLIC   tra_bbl_dif   !  routine called by trcbbl.F90
49   PUBLIC   tra_bbl_adv   !  -          -          -              -
50   PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
51
52   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
53
54   !                                !!* Namelist nambbl *
55   INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
56   INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
57   !                                            !  =1 : advective bbl using the bottom ocean velocity
58   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
59   REAL(wp), PUBLIC ::   rn_ahtbbl   !: along slope bbl diffusive coefficient [m2/s]
60   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
61
62   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
63
64   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
65   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
66
67   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
68   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
69   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
70   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_1, ahv_bbl_1   ! diffusive bbl flux coefficients at u and v-points
71   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)
72
73   !! * Substitutions
74#  include "domzgr_substitute.h90"
75#  include "vectopt_loop_substitute.h90"
76   !!----------------------------------------------------------------------
77   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
78   !! $Id$
79   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
80   !!----------------------------------------------------------------------
81CONTAINS
82
83   INTEGER FUNCTION tra_bbl_alloc()
84      !!----------------------------------------------------------------------
85      !!                ***  FUNCTION tra_bbl_alloc  ***
86      !!----------------------------------------------------------------------
87      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
88         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
89         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
90         &      ahu_bbl_1(jpi,jpj) , ahv_bbl_1(jpi,jpj) ,                                          &
91         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
92         !
93      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
94      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
95   END FUNCTION tra_bbl_alloc
96
97
98   SUBROUTINE tra_bbl( kt )
99      !!----------------------------------------------------------------------
100      !!                  ***  ROUTINE bbl  ***
101      !!
102      !! ** Purpose :   Compute the before tracer (t & s) trend associated
103      !!              with the bottom boundary layer and add it to the general
104      !!              trend of tracer equations.
105      !!
106      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
107      !!              diffusive and/or advective contribution to the tracer trend
108      !!              is added to the general tracer trend
109      !!----------------------------------------------------------------------
110      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
111      !
112      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
113      !!----------------------------------------------------------------------
114      !
115      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
116      !
117      IF( l_trdtra )   THEN                         !* Save ta and sa trends
118         ALLOCATE( ztrdt (1:jpi, 1:jpj, 1:jpk))
119         ALLOCATE( ztrds (1:jpi, 1:jpj, 1:jpk))
120         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
121         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
122      ENDIF
123
124      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
125
126      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
127         !
128         CALL tra_bbl_dif( tsb, tsa, jpts )
129         IF( ln_ctl )  &
130         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
131            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
132         ! lateral boundary conditions ; just need for outputs
133         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
134         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
135         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
136         !
137      END IF
138
139      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
140         !
141         CALL tra_bbl_adv( tsb, tsa, jpts )
142         IF(ln_ctl)   &
143         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
144            &          tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
145         ! lateral boundary conditions ; just need for outputs
146         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
147         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
148         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
149         !
150      END IF
151
152      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
153         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
154         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
155         CALL trd_tra( kt, 'TRA', jp_tem, jptra_bbl, ztrdt )
156         CALL trd_tra( kt, 'TRA', jp_sal, jptra_bbl, ztrds )
157         DEALLOCATE( ztrdt, ztrds )
158      ENDIF
159      !
160      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
161      !
162   END SUBROUTINE tra_bbl
163
164
165   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
166      !!----------------------------------------------------------------------
167      !!                  ***  ROUTINE tra_bbl_dif  ***
168      !!
169      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
170      !!                advection terms.
171      !!
172      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
173      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
174      !!      along bottom slope gradient) an additional lateral 2nd order
175      !!      diffusion along the bottom slope is added to the general
176      !!      tracer trend, otherwise the additional trend is set to 0.
177      !!      A typical value of ahbt is 2000 m2/s (equivalent to
178      !!      a downslope velocity of 20 cm/s if the condition for slope
179      !!      convection is satified)
180      !!
181      !! ** Action  :   pta   increased by the bbl diffusive trend
182      !!
183      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
184      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
185      !!----------------------------------------------------------------------
186      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
187      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
188      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
189      !
190      INTEGER  ::   ji, jj, jn   ! dummy loop indices
191      INTEGER  ::   ik           ! local integers
192      REAL(wp) ::   zbtr         ! local scalars
193      REAL(wp), ALLOCATABLE , DIMENSION(:,:) :: zptb
194      !!----------------------------------------------------------------------
195      !
196      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
197      !
198      ALLOCATE(zptb(1:jpi, 1:jpj))
199      !
200      ahu_bbl_1(:,:) = ahu_bbl(:,:)
201#if defined key_traldf_c2d || key_traldf_c3d
202      IF( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN
203          CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl )
204      ENDIF
205#else
206      IF ( ln_stopack .AND. nn_spp_ahubbl > 0 ) &
207         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
208                          'key_traldf_c2d or key_traldf_c3d')
209#endif
210
211
212      ahv_bbl_1(:,:) = ahv_bbl(:,:)
213#if defined key_traldf_c2d || key_traldf_c3d
214      IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN
215          CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl )
216      ENDIF
217#else
218      IF ( ln_stopack .AND. nn_spp_ahvbbl > 0 ) &
219         & CALL ctl_stop( 'tra_bbl_dif: parameter perturbation will only work with '// &
220                          'key_traldf_c2d or key_traldf_c3d')
221#endif
222
223      !
224      DO jn = 1, kjpt                                     ! tracer loop
225         !                                                ! ===========
226         DO jj = 1, jpj
227            DO ji = 1, jpi
228               ik = mbkt(ji,jj)                              ! bottom T-level index
229               zptb(ji,jj) = ptb(ji,jj,ik,jn)       ! bottom before T and S
230            END DO
231         END DO
232         !
233         DO jj = 2, jpjm1                                    ! Compute the trend
234            DO ji = 2, jpim1
235               ik = mbkt(ji,jj)                              ! bottom T-level index
236               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
237               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
238                  &               + (   ahu_bbl_1(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
239                  &                   - ahu_bbl_1(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
240                  &                   + ahv_bbl_1(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
241                  &                   - ahv_bbl_1(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
242            END DO
243         END DO
244         !                                                  ! ===========
245      END DO                                                ! end tracer
246      !                                                     ! ===========
247      DEALLOCATE( zptb )
248      !
249      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
250      !
251   END SUBROUTINE tra_bbl_dif
252
253
254   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
255      !!----------------------------------------------------------------------
256      !!                  ***  ROUTINE trc_bbl  ***
257      !!
258      !! ** Purpose :   Compute the before passive tracer trend associated
259      !!     with the bottom boundary layer and add it to the general trend
260      !!     of tracer equations.
261      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
262      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
263      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
264      !!                       transport proportional to the along-slope density gradient
265      !!
266      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
267      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
268      !!----------------------------------------------------------------------
269      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
270      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
271      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
272      !
273      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
274      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
275      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
276      REAL(wp) ::   zbtr, ztra               ! local scalars
277      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
278      !!----------------------------------------------------------------------
279      !
280      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
281      !                                                          ! ===========
282      DO jn = 1, kjpt                                            ! tracer loop
283         !                                                       ! ===========
284         DO jj = 1, jpjm1
285            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
286               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
287                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
288                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
289                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
290                  zu_bbl = ABS( utr_bbl(ji,jj) )
291                  !
292                  !                                               ! up  -slope T-point (shelf bottom point)
293                  zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus)
294                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
295                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
296                  !
297                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
298                     zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk)
299                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
300                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
301                  END DO
302                  !
303                  zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud)
304                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
305                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
306               ENDIF
307               !
308               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
309                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
310                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
311                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
312                  zv_bbl = ABS( vtr_bbl(ji,jj) )
313                  !
314                  ! up  -slope T-point (shelf bottom point)
315                  zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs)
316                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
317                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
318                  !
319                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
320                     zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk)
321                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
322                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
323                  END DO
324                  !                                               ! down-slope T-point (deep bottom point)
325                  zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd)
326                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
327                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
328               ENDIF
329            END DO
330            !
331         END DO
332         !                                                  ! ===========
333      END DO                                                ! end tracer
334      !                                                     ! ===========
335      !
336      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
337      !
338   END SUBROUTINE tra_bbl_adv
339
340
341   SUBROUTINE bbl( kt, kit000, cdtype )
342      !!----------------------------------------------------------------------
343      !!                  ***  ROUTINE bbl  ***
344      !!
345      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
346      !!                advection terms.
347      !!
348      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
349      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
350      !!      along bottom slope gradient) an additional lateral 2nd order
351      !!      diffusion along the bottom slope is added to the general
352      !!      tracer trend, otherwise the additional trend is set to 0.
353      !!      A typical value of ahbt is 2000 m2/s (equivalent to
354      !!      a downslope velocity of 20 cm/s if the condition for slope
355      !!      convection is satified)
356      !!              * advective bbl (nn_bbl_adv=1 or 2) :
357      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
358      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
359      !!        i.e. transport proportional to the along-slope density gradient
360      !!
361      !!      NB: the along slope density gradient is evaluated using the
362      !!      local density (i.e. referenced at a common local depth).
363      !!
364      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
365      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
366      !!----------------------------------------------------------------------
367      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
368      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
369      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
370      !!
371      INTEGER  ::   ji, jj                    ! dummy loop indices
372      INTEGER  ::   ik                        ! local integers
373      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
374      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
375      REAL(wp) ::   za, zb, zgdrho            ! local scalars
376      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
377      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
378      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
379      !!----------------------------------------------------------------------
380      !
381      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
382      !
383      IF( kt == kit000 )  THEN
384         IF(lwp)  WRITE(numout,*)
385         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
386         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
387      ENDIF
388      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
389      DO jj = 1, jpj
390         DO ji = 1, jpi
391            ik = mbkt(ji,jj)                             ! bottom T-level index
392            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
393            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
394            !
395            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth
396            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
397            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
398         END DO
399      END DO
400      !
401      CALL eos_rab( zts, zdep, zab )
402      !
403      !                                   !-------------------!
404      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
405         !                                !-------------------!
406         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
407            DO ji = 1, fs_jpim1   ! vector opt.
408               !                                                   ! i-direction
409               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
410               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
411               !                                                         ! 2*masked bottom density gradient
412               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
413                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
414               !
415               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
416               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
417               !
418               !                                                   ! j-direction
419               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
420               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
421               !                                                         ! 2*masked bottom density gradient
422               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
423                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
424               !
425               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
426               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
427            END DO
428         END DO
429         !
430      ENDIF
431
432      !                                   !-------------------!
433      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
434         !                                !-------------------!
435         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
436         !
437         CASE( 1 )                                   != use of upper velocity
438            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
439               DO ji = 1, fs_jpim1   ! vector opt.
440                  !                                                  ! i-direction
441                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
442                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
443                  !                                                          ! 2*masked bottom density gradient
444                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
445                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
446                  !
447                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
448                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
449                  !
450                  !                                                          ! bbl velocity
451                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
452                  !
453                  !                                                  ! j-direction
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                  !                                                          ! 2*masked bottom density gradient
457                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
458                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
459                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
460                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
461                  !
462                  !                                                          ! bbl transport
463                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
464               END DO
465            END DO
466            !
467         CASE( 2 )                                 != bbl velocity = F( delta rho )
468            zgbbl = grav * rn_gambbl
469            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
470               DO ji = 1, fs_jpim1   ! vector opt.
471                  !                                                  ! i-direction
472                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
473                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
474                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
475                  !
476                  ikud = mbku_d(ji,jj)
477                  ikus = mbku(ji,jj)
478                  !
479                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
480                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
481                  !                                                          !   masked bottom density gradient
482                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
483                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
484                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
485                  !
486                  !                                                          ! bbl transport (down-slope direction)
487                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
488                  !
489                  !                                                  ! j-direction
490                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
491                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
492                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
493                  !
494                  ikvd = mbkv_d(ji,jj)
495                  ikvs = mbkv(ji,jj)
496                  !
497                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
498                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
499                  !                                                          !   masked bottom density gradient
500                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
501                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
502                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
503                  !
504                  !                                                          ! bbl transport (down-slope direction)
505                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
506               END DO
507            END DO
508         END SELECT
509         !
510      ENDIF
511      !
512      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
513      !
514   END SUBROUTINE bbl
515
516
517   SUBROUTINE tra_bbl_init
518      !!----------------------------------------------------------------------
519      !!                  ***  ROUTINE tra_bbl_init  ***
520      !!
521      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
522      !!
523      !! ** Method  :   Read the nambbl namelist and check the parameters
524      !!              called by nemo_init at the first timestep (kit000)
525      !!----------------------------------------------------------------------
526      INTEGER ::   ji, jj               ! dummy loop indices
527      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
528      INTEGER ::   ios                  !   -      -
529      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
530      !!
531      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
532      !!----------------------------------------------------------------------
533      !
534      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
535      !
536      CALL wrk_alloc( jpi, jpj, zmbk )
537      !
538
539      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
540      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
541901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
542
543      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
544      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
545902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
546      IF(lwm) WRITE ( numond, nambbl )
547      !
548      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
549      !
550      IF(lwp) THEN                   !* Parameter control and print
551         WRITE(numout,*)
552         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
553         WRITE(numout,*) '~~~~~~~~~~~~'
554         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
555         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
556         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
557         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
558         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
559      ENDIF
560
561      !                              ! allocate trabbl arrays
562      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
563
564      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
565      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
566
567      !                             !* vertical index of  "deep" bottom u- and v-points
568      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
569         DO ji = 1, jpim1
570            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
571            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
572         END DO
573      END DO
574      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
575      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
576      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
577
578      !! AXY (16/08/17): remove the following per George and Andrew bug-hunt
579      !!                                   !* sign of grad(H) at u- and v-points
580      !! mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
581      !! DO jj = 1, jpjm1
582      !!    DO ji = 1, jpim1
583      !!       mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
584      !!       mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
585      !!    END DO
586      !! END DO
587
588      !! AXY (16/08/17): add the following replacement per George and Andrew bug-hunt
589                                        !* sign of grad(H) at u- and  v-points; zero if grad(H) = 0
590      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
591      DO jj = 1, jpjm1
592         DO ji = 1, jpim1
593#if defined key_bbl_old_nonconserve
594             ! This key allows old (non conservative version) to be used for continuity of results
595             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
596             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
597#else
598            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
599               mgrhu(ji,jj) = INT(  SIGN( 1.e0, &
600               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
601            ENDIF
602            !
603            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
604               mgrhv(ji,jj) = INT(  SIGN( 1.e0, &
605               gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
606            ENDIF
607#endif
608         END DO
609      END DO
610
611      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
612         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
613            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
614            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
615         END DO
616      END DO
617      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
618
619      !                             !* masked diffusive flux coefficients
620      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
621      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
622
623      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
624         !
625         SELECT CASE ( jp_cfg )
626         CASE ( 2 )                          ! ORCA_R2
627            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
628            ii0 = 139   ;   ii1 = 140
629            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
630            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
631            !
632            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
633            ii0 = 161   ;   ii1 = 162
634            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
635            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
636            !
637         CASE ( 4 )                          ! ORCA_R4
638            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
639            ii0 =  70   ;   ii1 =  71
640            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
641            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
642         END SELECT
643         !
644      ENDIF
645      !
646      CALL wrk_dealloc( jpi, jpj, zmbk )
647      !
648      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
649      !
650   END SUBROUTINE tra_bbl_init
651
652#else
653   !!----------------------------------------------------------------------
654   !!   Dummy module :                      No bottom boundary layer scheme
655   !!----------------------------------------------------------------------
656   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
657CONTAINS
658   SUBROUTINE tra_bbl_init               ! Dummy routine
659   END SUBROUTINE tra_bbl_init
660   SUBROUTINE tra_bbl( kt )              ! Dummy routine
661      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
662   END SUBROUTINE tra_bbl
663#endif
664
665   !!======================================================================
666END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.