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 @ 9383

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

#2050 fixes and changes

File size: 35.2 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      CONTINUE
515      ENDIF
516      call mpp_bcast(ios)
517      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
518      IF(lwm) THEN
519         REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
520         READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
521902      CONTINUE
522      ENDIF
523      call mpp_bcast(ios)
524      IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
525
526      IF(lwm) WRITE ( numond, nambbl )
527      !
528      CALL bbl_namelist()
529
530      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
531      !
532      IF(lwp) THEN                   !* Parameter control and print
533         WRITE(numout,*)
534         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
535         WRITE(numout,*) '~~~~~~~~~~~~'
536         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
537         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
538         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
539         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
540         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
541      ENDIF
542
543      !                              ! allocate trabbl arrays
544      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
545
546      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
547      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
548
549      !                             !* vertical index of  "deep" bottom u- and v-points
550      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
551         DO ji = 1, jpim1
552            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
553            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
554         END DO
555      END DO
556      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
557      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
558      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
559
560                                        !* sign of grad(H) at u- and v-points
561      mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
562      DO jj = 1, jpjm1
563         DO ji = 1, jpim1
564            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
565            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
566         END DO
567      END DO
568
569      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
570         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
571            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
572            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
573         END DO
574      END DO
575      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
576
577      !                             !* masked diffusive flux coefficients
578      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
579      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
580
581
582      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
583         !
584         SELECT CASE ( jp_cfg )
585         CASE ( 2 )                          ! ORCA_R2
586            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
587            ii0 = 139   ;   ii1 = 140
588            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
589            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
590            !
591            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
592            ii0 = 161   ;   ii1 = 162
593            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
594            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
595            !
596         CASE ( 4 )                          ! ORCA_R4
597            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
598            ii0 =  70   ;   ii1 =  71
599            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
600            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
601         END SELECT
602         !
603      ENDIF
604      !
605      CALL wrk_dealloc( jpi, jpj, zmbk )
606      !
607      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
608      !
609   END SUBROUTINE tra_bbl_init
610
611   SUBROUTINE bbl_namelist()
612     !!---------------------------------------------------------------------
613     !!                   ***  ROUTINE bbl_namelist  ***
614     !!                     
615     !! ** Purpose :   Broadcast namelist variables read by procesor lwm
616     !!
617     !! ** Method  :   use lib_mpp
618     !!----------------------------------------------------------------------
619#if defined key_mpp_mpi
620      CALL mpp_bcast(nn_bbl_ldf)
621      CALL mpp_bcast(nn_bbl_adv)
622      CALL mpp_bcast(rn_ahtbbl)
623      CALL mpp_bcast(rn_gambbl)
624#endif
625   END SUBROUTINE bbl_namelist
626
627#else
628   !!----------------------------------------------------------------------
629   !!   Dummy module :                      No bottom boundary layer scheme
630   !!----------------------------------------------------------------------
631   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
632CONTAINS
633   SUBROUTINE tra_bbl_init               ! Dummy routine
634   END SUBROUTINE tra_bbl_init
635   SUBROUTINE tra_bbl( kt )              ! Dummy routine
636      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
637   END SUBROUTINE tra_bbl
638#endif
639
640   !!======================================================================
641END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.