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

source: branches/UKMO/test_moci_test_suite_namelist_read/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 9366

Last change on this file since 9366 was 9366, checked in by andmirek, 6 years ago

#2050 first version. Compiled OK in moci test suite

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