source: branches/UKMO/dev_r5518_GO6_package_FOAMv14/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 11442

Last change on this file since 11442 was 11442, checked in by mattmartin, 17 months ago

Introduction of stochastic physics in NEMO, based on Andrea Storto's code.
For details, see ticket https://code.metoffice.gov.uk/trac/utils/ticket/251.

File size: 36.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   !!----------------------------------------------------------------------
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( ln_stopack .AND. nn_spp_ahubbl > 0 ) THEN
202          CALL spp_gen(1, ahu_bbl_1, nn_spp_ahubbl, rn_ahubbl_sd, jk_spp_ahubbl )
203      ENDIF
204      ahv_bbl_1(:,:) = ahv_bbl(:,:)
205      IF( ln_stopack .AND. nn_spp_ahvbbl > 0 ) THEN
206          CALL spp_gen(1, ahv_bbl_1, nn_spp_ahvbbl, rn_ahvbbl_sd, jk_spp_ahvbbl )
207      ENDIF
208      !
209      DO jn = 1, kjpt                                     ! tracer loop
210         !                                                ! ===========
211         DO jj = 1, jpj
212            DO ji = 1, jpi
213               ik = mbkt(ji,jj)                              ! bottom T-level index
214               zptb(ji,jj) = ptb(ji,jj,ik,jn)       ! bottom before T and S
215            END DO
216         END DO
217         !               
218         DO jj = 2, jpjm1                                    ! Compute the trend
219            DO ji = 2, jpim1
220               ik = mbkt(ji,jj)                              ! bottom T-level index
221               zbtr = r1_e12t(ji,jj)  / fse3t(ji,jj,ik)
222               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
223                  &               + (   ahu_bbl_1(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
224                  &                   - ahu_bbl_1(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
225                  &                   + ahv_bbl_1(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
226                  &                   - ahv_bbl_1(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
227            END DO
228         END DO
229         !                                                  ! ===========
230      END DO                                                ! end tracer
231      !                                                     ! ===========
232      DEALLOCATE( zptb )
233      !
234      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
235      !
236   END SUBROUTINE tra_bbl_dif
237
238
239   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
240      !!----------------------------------------------------------------------
241      !!                  ***  ROUTINE trc_bbl  ***
242      !!
243      !! ** Purpose :   Compute the before passive tracer trend associated
244      !!     with the bottom boundary layer and add it to the general trend
245      !!     of tracer equations.
246      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
247      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
248      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
249      !!                       transport proportional to the along-slope density gradient
250      !!
251      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
252      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
253      !!----------------------------------------------------------------------
254      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
255      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
256      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
257      !
258      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
259      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
260      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
261      REAL(wp) ::   zbtr, ztra               ! local scalars
262      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
263      !!----------------------------------------------------------------------
264      !
265      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
266      !                                                          ! ===========
267      DO jn = 1, kjpt                                            ! tracer loop
268         !                                                       ! ===========
269         DO jj = 1, jpjm1
270            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
271               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
272                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
273                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
274                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
275                  zu_bbl = ABS( utr_bbl(ji,jj) )
276                  !
277                  !                                               ! up  -slope T-point (shelf bottom point)
278                  zbtr = r1_e12t(iis,jj) / fse3t(iis,jj,ikus)
279                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
280                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
281                  !
282                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
283                     zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,jk)
284                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
285                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
286                  END DO
287                  !
288                  zbtr = r1_e12t(iid,jj) / fse3t(iid,jj,ikud)
289                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
290                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
291               ENDIF
292               !
293               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
294                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
295                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
296                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
297                  zv_bbl = ABS( vtr_bbl(ji,jj) )
298                  !
299                  ! up  -slope T-point (shelf bottom point)
300                  zbtr = r1_e12t(ji,ijs) / fse3t(ji,ijs,ikvs)
301                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
302                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
303                  !
304                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
305                     zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,jk)
306                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
307                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
308                  END DO
309                  !                                               ! down-slope T-point (deep bottom point)
310                  zbtr = r1_e12t(ji,ijd) / fse3t(ji,ijd,ikvd)
311                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
312                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
313               ENDIF
314            END DO
315            !
316         END DO
317         !                                                  ! ===========
318      END DO                                                ! end tracer
319      !                                                     ! ===========
320      !
321      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
322      !
323   END SUBROUTINE tra_bbl_adv
324
325
326   SUBROUTINE bbl( kt, kit000, cdtype )
327      !!----------------------------------------------------------------------
328      !!                  ***  ROUTINE bbl  ***
329      !!
330      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
331      !!                advection terms.
332      !!
333      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
334      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
335      !!      along bottom slope gradient) an additional lateral 2nd order
336      !!      diffusion along the bottom slope is added to the general
337      !!      tracer trend, otherwise the additional trend is set to 0.
338      !!      A typical value of ahbt is 2000 m2/s (equivalent to
339      !!      a downslope velocity of 20 cm/s if the condition for slope
340      !!      convection is satified)
341      !!              * advective bbl (nn_bbl_adv=1 or 2) :
342      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
343      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
344      !!        i.e. transport proportional to the along-slope density gradient
345      !!
346      !!      NB: the along slope density gradient is evaluated using the
347      !!      local density (i.e. referenced at a common local depth).
348      !!
349      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
350      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
351      !!----------------------------------------------------------------------
352      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
353      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
354      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
355      !!
356      INTEGER  ::   ji, jj                    ! dummy loop indices
357      INTEGER  ::   ik                        ! local integers
358      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
359      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
360      REAL(wp) ::   za, zb, zgdrho            ! local scalars
361      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
362      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
363      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
364      !!----------------------------------------------------------------------
365      !
366      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
367      !
368      IF( kt == kit000 )  THEN
369         IF(lwp)  WRITE(numout,*)
370         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
371         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
372      ENDIF
373      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
374      DO jj = 1, jpj
375         DO ji = 1, jpi
376            ik = mbkt(ji,jj)                             ! bottom T-level index
377            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
378            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
379            !
380            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth
381            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
382            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
383         END DO
384      END DO
385      !
386      CALL eos_rab( zts, zdep, zab )
387      !
388      !                                   !-------------------!
389      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
390         !                                !-------------------!
391         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
392            DO ji = 1, fs_jpim1   ! vector opt.
393               !                                                   ! i-direction
394               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
395               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
396               !                                                         ! 2*masked bottom density gradient
397               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
398                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
399               !
400               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
401               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
402               !
403               !                                                   ! j-direction
404               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
405               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
406               !                                                         ! 2*masked bottom density gradient
407               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
408                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
409               !
410               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
411               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
412            END DO
413         END DO
414         !
415      ENDIF
416
417      !                                   !-------------------!
418      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
419         !                                !-------------------!
420         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
421         !
422         CASE( 1 )                                   != use of upper velocity
423            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
424               DO ji = 1, fs_jpim1   ! vector opt.
425                  !                                                  ! i-direction
426                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
427                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
428                  !                                                          ! 2*masked bottom density gradient
429                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
430                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
431                  !
432                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
433                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
434                  !
435                  !                                                          ! bbl velocity
436                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
437                  !
438                  !                                                  ! j-direction
439                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
440                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
441                  !                                                          ! 2*masked bottom density gradient
442                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
443                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
444                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
445                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
446                  !
447                  !                                                          ! bbl transport
448                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
449               END DO
450            END DO
451            !
452         CASE( 2 )                                 != bbl velocity = F( delta rho )
453            zgbbl = grav * rn_gambbl
454            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
455               DO ji = 1, fs_jpim1   ! vector opt.
456                  !                                                  ! i-direction
457                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
458                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
459                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
460                  !
461                  ikud = mbku_d(ji,jj)
462                  ikus = mbku(ji,jj)
463                  !
464                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
465                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
466                  !                                                          !   masked bottom density gradient
467                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
468                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
469                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
470                  !
471                  !                                                          ! bbl transport (down-slope direction)
472                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
473                  !
474                  !                                                  ! j-direction
475                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
476                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
477                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
478                  !
479                  ikvd = mbkv_d(ji,jj)
480                  ikvs = mbkv(ji,jj)
481                  !
482                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
483                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
484                  !                                                          !   masked bottom density gradient
485                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
486                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
487                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
488                  !
489                  !                                                          ! bbl transport (down-slope direction)
490                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
491               END DO
492            END DO
493         END SELECT
494         !
495      ENDIF
496      !
497      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
498      !
499   END SUBROUTINE bbl
500
501
502   SUBROUTINE tra_bbl_init
503      !!----------------------------------------------------------------------
504      !!                  ***  ROUTINE tra_bbl_init  ***
505      !!
506      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
507      !!
508      !! ** Method  :   Read the nambbl namelist and check the parameters
509      !!              called by nemo_init at the first timestep (kit000)
510      !!----------------------------------------------------------------------
511      INTEGER ::   ji, jj               ! dummy loop indices
512      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
513      INTEGER ::   ios                  !   -      -
514      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
515      !!
516      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
517      !!----------------------------------------------------------------------
518      !
519      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
520      !
521      CALL wrk_alloc( jpi, jpj, zmbk )
522      !
523
524      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
525      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
526901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
527
528      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
529      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
530902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
531      IF(lwm) WRITE ( numond, nambbl )
532      !
533      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
534      !
535      IF(lwp) THEN                   !* Parameter control and print
536         WRITE(numout,*)
537         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
538         WRITE(numout,*) '~~~~~~~~~~~~'
539         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
540         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
541         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
542         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
543         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
544      ENDIF
545
546      !                              ! allocate trabbl arrays
547      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
548
549      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
550      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
551
552      !                             !* vertical index of  "deep" bottom u- and v-points
553      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
554         DO ji = 1, jpim1
555            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
556            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
557         END DO
558      END DO
559      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
560      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
561      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
562
563      !! AXY (16/08/17): remove the following per George and Andrew bug-hunt
564      !!                                   !* sign of grad(H) at u- and v-points
565      !! mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
566      !! DO jj = 1, jpjm1
567      !!    DO ji = 1, jpim1
568      !!       mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
569      !!       mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
570      !!    END DO
571      !! END DO
572
573      !! AXY (16/08/17): add the following replacement per George and Andrew bug-hunt
574                                        !* sign of grad(H) at u- and  v-points; zero if grad(H) = 0
575      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
576      DO jj = 1, jpjm1
577         DO ji = 1, jpim1
578#if defined key_bbl_old_nonconserve
579             ! This key allows old (non conservative version) to be used for continuity of results
580             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
581             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
582#else
583            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
584               mgrhu(ji,jj) = INT(  SIGN( 1.e0, &
585               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
586            ENDIF
587            !     
588            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
589               mgrhv(ji,jj) = INT(  SIGN( 1.e0, &
590               gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
591            ENDIF
592#endif
593         END DO
594      END DO
595
596      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
597         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
598            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
599            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
600         END DO
601      END DO
602      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
603
604      !                             !* masked diffusive flux coefficients
605      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
606      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
607
608      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
609         !
610         SELECT CASE ( jp_cfg )
611         CASE ( 2 )                          ! ORCA_R2
612            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
613            ii0 = 139   ;   ii1 = 140
614            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
615            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
616            !
617            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
618            ii0 = 161   ;   ii1 = 162
619            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
620            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
621            !
622         CASE ( 4 )                          ! ORCA_R4
623            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
624            ii0 =  70   ;   ii1 =  71
625            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
626            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
627         END SELECT
628         !
629      ENDIF
630      !
631      CALL wrk_dealloc( jpi, jpj, zmbk )
632      !
633      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
634      !
635   END SUBROUTINE tra_bbl_init
636
637#else
638   !!----------------------------------------------------------------------
639   !!   Dummy module :                      No bottom boundary layer scheme
640   !!----------------------------------------------------------------------
641   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
642CONTAINS
643   SUBROUTINE tra_bbl_init               ! Dummy routine
644   END SUBROUTINE tra_bbl_init
645   SUBROUTINE tra_bbl( kt )              ! Dummy routine
646      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
647   END SUBROUTINE tra_bbl
648#endif
649
650   !!======================================================================
651END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.