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

source: branches/2016/dev_r6519_HPC_4/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 7037

Last change on this file since 7037 was 7037, checked in by mocavero, 8 years ago

ORCA2_LIM_PISCES hybrid version update

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