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

source: branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 15603

Last change on this file since 15603 was 15603, checked in by mattmartin, 3 years ago

Updated NEMO branch for coupled NWP at GO6 to include stochastic model perturbations.
For more info see ticket: https://code.metoffice.gov.uk/trac/nwpscience/ticket/1125.

File size: 36.9 KB
Line 
1MODULE trabbl
2   !!==============================================================================
3   !!                       ***  MODULE  trabbl  ***
4   !! Ocean physics :  advective and/or diffusive bottom boundary layer scheme
5   !!==============================================================================
6   !! History :  OPA  ! 1996-06  (L. Mortier)  Original code
7   !!            8.0  ! 1997-11  (G. Madec)    Optimization
8   !!   NEMO     1.0  ! 2002-08  (G. Madec)  free form + modules
9   !!             -   ! 2004-01  (A. de Miranda, G. Madec, J.M. Molines ) add advective bbl
10   !!            3.3  ! 2009-11  (G. Madec)  merge trabbl and trabbl_adv + style + optimization
11   !!             -   ! 2010-04  (G. Madec)  Campin & Goosse advective bbl
12   !!             -   ! 2010-06  (C. Ethe, G. Madec)  merge TRA-TRC
13   !!             -   ! 2010-11  (G. Madec) add mbk. arrays associated to the deepest ocean level
14   !!             -   ! 2013-04  (F. Roquet, G. Madec)  use of eosbn2 instead of local hard coded alpha and beta
15   !!----------------------------------------------------------------------
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         IF(lflush) CALL flush(numout)
388      ENDIF
389      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
390      DO jj = 1, jpj
391         DO ji = 1, jpi
392            ik = mbkt(ji,jj)                             ! bottom T-level index
393            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
394            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
395            !
396            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth
397            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
398            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
399         END DO
400      END DO
401      !
402      CALL eos_rab( zts, zdep, zab )
403      !
404      !                                   !-------------------!
405      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
406         !                                !-------------------!
407         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
408            DO ji = 1, fs_jpim1   ! vector opt.
409               !                                                   ! i-direction
410               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
411               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
412               !                                                         ! 2*masked bottom density gradient
413               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
414                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
415               !
416               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
417               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
418               !
419               !                                                   ! j-direction
420               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
421               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
422               !                                                         ! 2*masked bottom density gradient
423               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
424                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
425               !
426               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
427               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
428            END DO
429         END DO
430         !
431      ENDIF
432
433      !                                   !-------------------!
434      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
435         !                                !-------------------!
436         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
437         !
438         CASE( 1 )                                   != use of upper velocity
439            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
440               DO ji = 1, fs_jpim1   ! vector opt.
441                  !                                                  ! i-direction
442                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
443                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
444                  !                                                          ! 2*masked bottom density gradient
445                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
446                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
447                  !
448                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
449                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
450                  !
451                  !                                                          ! bbl velocity
452                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
453                  !
454                  !                                                  ! j-direction
455                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
456                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
457                  !                                                          ! 2*masked bottom density gradient
458                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
459                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
460                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
461                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
462                  !
463                  !                                                          ! bbl transport
464                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
465               END DO
466            END DO
467            !
468         CASE( 2 )                                 != bbl velocity = F( delta rho )
469            zgbbl = grav * rn_gambbl
470            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
471               DO ji = 1, fs_jpim1   ! vector opt.
472                  !                                                  ! i-direction
473                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
474                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
475                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
476                  !
477                  ikud = mbku_d(ji,jj)
478                  ikus = mbku(ji,jj)
479                  !
480                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
481                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
482                  !                                                          !   masked bottom density gradient
483                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
484                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
485                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
486                  !
487                  !                                                          ! bbl transport (down-slope direction)
488                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
489                  !
490                  !                                                  ! j-direction
491                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
492                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
493                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
494                  !
495                  ikvd = mbkv_d(ji,jj)
496                  ikvs = mbkv(ji,jj)
497                  !
498                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
499                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
500                  !                                                          !   masked bottom density gradient
501                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
502                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
503                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
504                  !
505                  !                                                          ! bbl transport (down-slope direction)
506                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
507               END DO
508            END DO
509         END SELECT
510         !
511      ENDIF
512      !
513      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
514      !
515   END SUBROUTINE bbl
516
517
518   SUBROUTINE tra_bbl_init
519      !!----------------------------------------------------------------------
520      !!                  ***  ROUTINE tra_bbl_init  ***
521      !!
522      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
523      !!
524      !! ** Method  :   Read the nambbl namelist and check the parameters
525      !!              called by nemo_init at the first timestep (kit000)
526      !!----------------------------------------------------------------------
527      INTEGER ::   ji, jj               ! dummy loop indices
528      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
529      INTEGER ::   ios                  !   -      -
530      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
531      !!
532      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
533      !!----------------------------------------------------------------------
534      !
535      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
536      !
537      CALL wrk_alloc( jpi, jpj, zmbk )
538      !
539
540      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
541      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
542901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
543
544      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
545      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
546902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
547      IF(lwm .AND. nprint > 2) WRITE ( numond, nambbl )
548      !
549      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
550      !
551      IF(lwp) THEN                   !* Parameter control and print
552         WRITE(numout,*)
553         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
554         WRITE(numout,*) '~~~~~~~~~~~~'
555         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
556         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
557         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
558         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
559         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
560         IF(lflush) CALL flush(numout)
561      ENDIF
562
563      !                              ! allocate trabbl arrays
564      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
565
566      IF(lwp) THEN
567         IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
568         IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
569         IF(lflush) CALL flush(numout)
570      ENDIF
571
572      !                             !* vertical index of  "deep" bottom u- and v-points
573      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
574         DO ji = 1, jpim1
575            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
576            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
577         END DO
578      END DO
579      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
580      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
581      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
582
583      !! AXY (16/08/17): remove the following per George and Andrew bug-hunt
584      !!                                   !* sign of grad(H) at u- and v-points
585      !! mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
586      !! DO jj = 1, jpjm1
587      !!    DO ji = 1, jpim1
588      !!       mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
589      !!       mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
590      !!    END DO
591      !! END DO
592
593      !! AXY (16/08/17): add the following replacement per George and Andrew bug-hunt
594                                        !* sign of grad(H) at u- and  v-points; zero if grad(H) = 0
595      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
596      DO jj = 1, jpjm1
597         DO ji = 1, jpim1
598#if defined key_bbl_old_nonconserve
599             ! This key allows old (non conservative version) to be used for continuity of results
600             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
601             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
602#else
603            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
604               mgrhu(ji,jj) = INT(  SIGN( 1.e0, &
605               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
606            ENDIF
607            !
608            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
609               mgrhv(ji,jj) = INT(  SIGN( 1.e0, &
610               gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
611            ENDIF
612#endif
613         END DO
614      END DO
615
616      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
617         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
618            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
619            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
620         END DO
621      END DO
622      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
623
624      !                             !* masked diffusive flux coefficients
625      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
626      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
627
628      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
629         !
630         SELECT CASE ( jp_cfg )
631         CASE ( 2 )                          ! ORCA_R2
632            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
633            ii0 = 139   ;   ii1 = 140
634            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
635            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
636            !
637            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
638            ii0 = 161   ;   ii1 = 162
639            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
640            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
641            !
642         CASE ( 4 )                          ! ORCA_R4
643            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
644            ii0 =  70   ;   ii1 =  71
645            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
646            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
647         END SELECT
648         !
649      ENDIF
650      !
651      CALL wrk_dealloc( jpi, jpj, zmbk )
652      !
653      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
654      !
655   END SUBROUTINE tra_bbl_init
656
657#else
658   !!----------------------------------------------------------------------
659   !!   Dummy module :                      No bottom boundary layer scheme
660   !!----------------------------------------------------------------------
661   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
662CONTAINS
663   SUBROUTINE tra_bbl_init               ! Dummy routine
664   END SUBROUTINE tra_bbl_init
665   SUBROUTINE tra_bbl( kt )              ! Dummy routine
666      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
667   END SUBROUTINE tra_bbl
668#endif
669
670   !!======================================================================
671END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.