source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 20 months ago

GMED 450 add flush after prints

File size: 35.8 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
51   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .TRUE.    !: bottom boundary layer flag
52
53   !                                !!* Namelist nambbl *
54   INTEGER , PUBLIC ::   nn_bbl_ldf  !: =1   : diffusive bbl or not (=0)
55   INTEGER , PUBLIC ::   nn_bbl_adv  !: =1/2 : advective bbl or not (=0)
56   !                                            !  =1 : advective bbl using the bottom ocean velocity
57   !                                            !  =2 :     -      -  using utr_bbl proportional to grad(rho)
58   REAL(wp), PUBLIC ::   rn_ahtbbl   !: along slope bbl diffusive coefficient [m2/s]
59   REAL(wp), PUBLIC ::   rn_gambbl   !: lateral coeff. for bottom boundary layer scheme [s]
60
61   LOGICAL , PUBLIC ::   l_bbl       !: flag to compute bbl diffu. flux coef and transport
62
63   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
64   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
65
66   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
67   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
68   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
69   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)
70
71   !! * Substitutions
72#  include "domzgr_substitute.h90"
73#  include "vectopt_loop_substitute.h90"
74   !!----------------------------------------------------------------------
75   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
76   !! $Id$
77   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   INTEGER FUNCTION tra_bbl_alloc()
82      !!----------------------------------------------------------------------
83      !!                ***  FUNCTION tra_bbl_alloc  ***
84      !!----------------------------------------------------------------------
85      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d  (jpi,jpj) , mgrhu(jpi,jpj) ,     &
86         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d  (jpi,jpj) , mgrhv(jpi,jpj) ,     &
87         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                          &
88         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                      STAT=tra_bbl_alloc )
89         !
90      IF( lk_mpp            )   CALL mpp_sum ( tra_bbl_alloc )
91      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
92   END FUNCTION tra_bbl_alloc
93
94
95   SUBROUTINE tra_bbl( kt )
96      !!----------------------------------------------------------------------
97      !!                  ***  ROUTINE bbl  ***
98      !!
99      !! ** Purpose :   Compute the before tracer (t & s) trend associated
100      !!              with the bottom boundary layer and add it to the general
101      !!              trend of tracer equations.
102      !!
103      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
104      !!              diffusive and/or advective contribution to the tracer trend
105      !!              is added to the general tracer trend
106      !!----------------------------------------------------------------------
107      INTEGER, INTENT( in ) ::   kt   ! ocean time-step
108      !
109      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ztrdt, ztrds
110      !!----------------------------------------------------------------------
111      !
112      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
113      !
114      IF( l_trdtra )   THEN                         !* Save ta and sa trends
115         ALLOCATE( ztrdt (1:jpi, 1:jpj, 1:jpk))
116         ALLOCATE( ztrds (1:jpi, 1:jpj, 1:jpk))
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         DEALLOCATE( 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), ALLOCATABLE , DIMENSION(:,:) :: zptb
191      !!----------------------------------------------------------------------
192      !
193      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
194      !
195      ALLOCATE(zptb(1:jpi, 1:jpj))
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      DEALLOCATE( 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         IF(lflush) CALL flush(numout)
361      ENDIF
362      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
363      DO jj = 1, jpj
364         DO ji = 1, jpi
365            ik = mbkt(ji,jj)                             ! bottom T-level index
366            zts (ji,jj,jp_tem) = tsb(ji,jj,ik,jp_tem)    ! bottom before T and S
367            zts (ji,jj,jp_sal) = tsb(ji,jj,ik,jp_sal)
368            !
369            zdep(ji,jj) = fsdept(ji,jj,ik)               ! bottom T-level reference depth
370            zub (ji,jj) = un(ji,jj,mbku(ji,jj))          ! bottom velocity
371            zvb (ji,jj) = vn(ji,jj,mbkv(ji,jj))
372         END DO
373      END DO
374      !
375      CALL eos_rab( zts, zdep, zab )
376      !
377      !                                   !-------------------!
378      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
379         !                                !-------------------!
380         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
381            DO ji = 1, fs_jpim1   ! vector opt.
382               !                                                   ! i-direction
383               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
384               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
385               !                                                         ! 2*masked bottom density gradient
386               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
387                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
388               !
389               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
390               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
391               !
392               !                                                   ! j-direction
393               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
394               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
395               !                                                         ! 2*masked bottom density gradient
396               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
397                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
398               !
399               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
400               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
401            END DO
402         END DO
403         !
404      ENDIF
405
406      !                                   !-------------------!
407      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
408         !                                !-------------------!
409         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
410         !
411         CASE( 1 )                                   != use of upper velocity
412            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
413               DO ji = 1, fs_jpim1   ! vector opt.
414                  !                                                  ! i-direction
415                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
416                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
417                  !                                                          ! 2*masked bottom density gradient
418                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
419                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
420                  !
421                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
422                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
423                  !
424                  !                                                          ! bbl velocity
425                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
426                  !
427                  !                                                  ! j-direction
428                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
429                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
430                  !                                                          ! 2*masked bottom density gradient
431                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
432                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
433                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
434                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
435                  !
436                  !                                                          ! bbl transport
437                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
438               END DO
439            END DO
440            !
441         CASE( 2 )                                 != bbl velocity = F( delta rho )
442            zgbbl = grav * rn_gambbl
443            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
444               DO ji = 1, fs_jpim1   ! vector opt.
445                  !                                                  ! i-direction
446                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
447                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
448                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
449                  !
450                  ikud = mbku_d(ji,jj)
451                  ikus = mbku(ji,jj)
452                  !
453                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
454                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
455                  !                                                          !   masked bottom density gradient
456                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
457                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
458                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
459                  !
460                  !                                                          ! bbl transport (down-slope direction)
461                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
462                  !
463                  !                                                  ! j-direction
464                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
465                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
466                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
467                  !
468                  ikvd = mbkv_d(ji,jj)
469                  ikvs = mbkv(ji,jj)
470                  !
471                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
472                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
473                  !                                                          !   masked bottom density gradient
474                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
475                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
476                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
477                  !
478                  !                                                          ! bbl transport (down-slope direction)
479                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
480               END DO
481            END DO
482         END SELECT
483         !
484      ENDIF
485      !
486      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
487      !
488   END SUBROUTINE bbl
489
490
491   SUBROUTINE tra_bbl_init
492      !!----------------------------------------------------------------------
493      !!                  ***  ROUTINE tra_bbl_init  ***
494      !!
495      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
496      !!
497      !! ** Method  :   Read the nambbl namelist and check the parameters
498      !!              called by nemo_init at the first timestep (kit000)
499      !!----------------------------------------------------------------------
500      INTEGER ::   ji, jj               ! dummy loop indices
501      INTEGER ::   ii0, ii1, ij0, ij1   ! local integer
502      INTEGER ::   ios                  !   -      -
503      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
504      !!
505      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
506      !!----------------------------------------------------------------------
507      !
508      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
509      !
510      CALL wrk_alloc( jpi, jpj, zmbk )
511      !
512
513      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
514      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
515901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
516
517      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
518      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
519902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
520      IF(lwm .AND. nprint > 2) WRITE ( numond, nambbl )
521      !
522      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
523      !
524      IF(lwp) THEN                   !* Parameter control and print
525         WRITE(numout,*)
526         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
527         WRITE(numout,*) '~~~~~~~~~~~~'
528         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
529         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
530         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
531         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
532         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
533         IF(lflush) CALL flush(numout)
534      ENDIF
535
536      !                              ! allocate trabbl arrays
537      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
538
539      IF(lwp) THEN
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         IF(lflush) CALL flush(numout)
543      ENDIF
544
545      !                             !* vertical index of  "deep" bottom u- and v-points
546      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
547         DO ji = 1, jpim1
548            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
549            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
550         END DO
551      END DO
552      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
553      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
554      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
555
556      !! AXY (16/08/17): remove the following per George and Andrew bug-hunt
557      !!                                   !* sign of grad(H) at u- and v-points
558      !! mgrhu(jpi,:) = 0   ;   mgrhu(:,jpj) = 0   ;   mgrhv(jpi,:) = 0   ;   mgrhv(:,jpj) = 0
559      !! DO jj = 1, jpjm1
560      !!    DO ji = 1, jpim1
561      !!       mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
562      !!       mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
563      !!    END DO
564      !! END DO
565
566      !! AXY (16/08/17): add the following replacement per George and Andrew bug-hunt
567                                        !* sign of grad(H) at u- and  v-points; zero if grad(H) = 0
568      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
569      DO jj = 1, jpjm1
570         DO ji = 1, jpim1
571#if defined key_bbl_old_nonconserve
572             ! This key allows old (non conservative version) to be used for continuity of results
573             mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
574             mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
575#else
576            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
577               mgrhu(ji,jj) = INT(  SIGN( 1.e0, &
578               gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
579            ENDIF
580            !     
581            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
582               mgrhv(ji,jj) = INT(  SIGN( 1.e0, &
583               gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
584            ENDIF
585#endif
586         END DO
587      END DO
588
589      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
590         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
591            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
592            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
593         END DO
594      END DO
595      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
596
597      !                             !* masked diffusive flux coefficients
598      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
599      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
600
601
602      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
603         !
604         SELECT CASE ( jp_cfg )
605         CASE ( 2 )                          ! ORCA_R2
606            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
607            ii0 = 139   ;   ii1 = 140
608            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
609            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
610            !
611            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
612            ii0 = 161   ;   ii1 = 162
613            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
614            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
615            !
616         CASE ( 4 )                          ! ORCA_R4
617            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
618            ii0 =  70   ;   ii1 =  71
619            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
620            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
621         END SELECT
622         !
623      ENDIF
624      !
625      CALL wrk_dealloc( jpi, jpj, zmbk )
626      !
627      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
628      !
629   END SUBROUTINE tra_bbl_init
630
631#else
632   !!----------------------------------------------------------------------
633   !!   Dummy module :                      No bottom boundary layer scheme
634   !!----------------------------------------------------------------------
635   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
636CONTAINS
637   SUBROUTINE tra_bbl_init               ! Dummy routine
638   END SUBROUTINE tra_bbl_init
639   SUBROUTINE tra_bbl( kt )              ! Dummy routine
640      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
641   END SUBROUTINE tra_bbl
642#endif
643
644   !!======================================================================
645END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.