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

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago

All wrk_alloc removed

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