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 NEMO/branches/2020/dev_r13508_HPC-09_communications_cleanup/src/OCE/TRA – NEMO

source: NEMO/branches/2020/dev_r13508_HPC-09_communications_cleanup/src/OCE/TRA/trabbl.F90 @ 13619

Last change on this file since 13619 was 13619, checked in by francesca, 4 years ago

cleanup of the TRA modules - first step

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