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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trabbl.F90 @ 11822

Last change on this file since 11822 was 11822, checked in by acc, 4 years ago

Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Sette tested updates to branch to align with trunk changes between 10721 and 11740. Sette tests are passing but results differ from branch before these changes (except for GYRE_PISCES and VORTEX) and branch results already differed from trunk because of algorithmic fixes. Will need more checks to confirm correctness.
  • 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   !!            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 "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
72   !! $Id$
73   !! Software governed by the CeCILL license (see ./LICENSE)
74   !!----------------------------------------------------------------------
75CONTAINS
76
77   INTEGER FUNCTION tra_bbl_alloc()
78      !!----------------------------------------------------------------------
79      !!                ***  FUNCTION tra_bbl_alloc  ***
80      !!----------------------------------------------------------------------
81      ALLOCATE( utr_bbl  (jpi,jpj) , ahu_bbl  (jpi,jpj) , mbku_d(jpi,jpj) , mgrhu(jpi,jpj) ,     &
82         &      vtr_bbl  (jpi,jpj) , ahv_bbl  (jpi,jpj) , mbkv_d(jpi,jpj) , mgrhv(jpi,jpj) ,     &
83         &      ahu_bbl_0(jpi,jpj) , ahv_bbl_0(jpi,jpj) ,                                        &
84         &      e3u_bbl_0(jpi,jpj) , e3v_bbl_0(jpi,jpj) ,                                    STAT=tra_bbl_alloc )
85         !
86      CALL mpp_sum ( 'trabbl', tra_bbl_alloc )
87      IF( tra_bbl_alloc > 0 )   CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
88   END FUNCTION tra_bbl_alloc
89
90
91   SUBROUTINE tra_bbl( kt, Kbb, Kmm, pts, Krhs )
92      !!----------------------------------------------------------------------
93      !!                  ***  ROUTINE bbl  ***
94      !!
95      !! ** Purpose :   Compute the before tracer (t & s) trend associated
96      !!              with the bottom boundary layer and add it to the general
97      !!              trend of tracer equations.
98      !!
99      !! ** Method  :   Depending on namtra_bbl namelist parameters the bbl
100      !!              diffusive and/or advective contribution to the tracer trend
101      !!              is added to the general tracer trend
102      !!----------------------------------------------------------------------
103      INTEGER,                                   INTENT(in   ) :: kt              ! ocean time-step
104      INTEGER,                                   INTENT(in   ) :: Kbb, Kmm, Krhs  ! time level indices
105      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts             ! active tracers and RHS of tracer equation
106      !
107      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdt, ztrds
108      !!----------------------------------------------------------------------
109      !
110      IF( ln_timing )   CALL timing_start( 'tra_bbl')
111      !
112      IF( l_trdtra )   THEN                         !* Save the T-S input trends
113         ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
114         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs)
115         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs)
116      ENDIF
117
118      IF( l_bbl )   CALL bbl( kt, nit000, 'TRA', Kbb, Kmm )   !* bbl coef. and transport (only if not already done in trcbbl)
119
120      IF( nn_bbl_ldf == 1 ) THEN                    !* Diffusive bbl
121         !
122         CALL tra_bbl_dif( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
123         IF( ln_ctl )  &
124         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
125            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
126         ! lateral boundary conditions ; just need for outputs
127         CALL lbc_lnk_multi( 'trabbl', ahu_bbl, 'U', 1. , ahv_bbl, 'V', 1. )
128         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
129         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
130         !
131      ENDIF
132      !
133      IF( nn_bbl_adv /= 0 ) THEN                    !* Advective bbl
134         !
135         CALL tra_bbl_adv( pts(:,:,:,:,Kbb), pts(:,:,:,:,Krhs), jpts, Kmm )
136         IF(ln_ctl)   &
137         CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
138            &          tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
139         ! lateral boundary conditions ; just need for outputs
140         CALL lbc_lnk_multi( 'trabbl', utr_bbl, 'U', 1. , vtr_bbl, 'V', 1. )
141         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
142         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
143         !
144      ENDIF
145
146      IF( l_trdtra )   THEN                      ! send the trends for further diagnostics
147         ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
148         ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
149         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_bbl, ztrdt )
150         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_bbl, ztrds )
151         DEALLOCATE( ztrdt, ztrds )
152      ENDIF
153      !
154      IF( ln_timing )  CALL timing_stop( 'tra_bbl')
155      !
156   END SUBROUTINE tra_bbl
157
158
159   SUBROUTINE tra_bbl_dif( pt, pt_rhs, kjpt, Kmm )
160      !!----------------------------------------------------------------------
161      !!                  ***  ROUTINE tra_bbl_dif  ***
162      !!
163      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
164      !!                advection terms.
165      !!
166      !! ** Method  : * diffusive bbl only (nn_bbl_ldf=1) :
167      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
168      !!      along bottom slope gradient) an additional lateral 2nd order
169      !!      diffusion along the bottom slope is added to the general
170      !!      tracer trend, otherwise the additional trend is set to 0.
171      !!      A typical value of ahbt is 2000 m2/s (equivalent to
172      !!      a downslope velocity of 20 cm/s if the condition for slope
173      !!      convection is satified)
174      !!
175      !! ** Action  :   pt_rhs   increased by the bbl diffusive trend
176      !!
177      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
178      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
179      !!----------------------------------------------------------------------
180      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
181      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
182      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
183      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
184      !
185      INTEGER  ::   ji, jj, jn   ! dummy loop indices
186      INTEGER  ::   ik           ! local integers
187      REAL(wp) ::   zbtr         ! local scalars
188      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! workspace
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) = pt(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               pt_rhs(ji,jj,ik,jn) = pt_rhs(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(ji,jj,ik,Kmm)
209            END DO
210         END DO
211         !                                                  ! ===========
212      END DO                                                ! end tracer
213      !                                                     ! ===========
214   END SUBROUTINE tra_bbl_dif
215
216
217   SUBROUTINE tra_bbl_adv( pt, pt_rhs, kjpt, Kmm )
218      !!----------------------------------------------------------------------
219      !!                  ***  ROUTINE trc_bbl  ***
220      !!
221      !! ** Purpose :   Compute the before passive tracer trend associated
222      !!     with the bottom boundary layer and add it to the general trend
223      !!     of tracer equations.
224      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
225      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
226      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
227      !!                       transport proportional to the along-slope density gradient
228      !!
229      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
230      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
231      !!----------------------------------------------------------------------
232      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
233      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   pt     ! before and now tracer fields
234      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pt_rhs ! tracer trend
235      INTEGER                              , INTENT(in   ) ::   Kmm    ! time level indices
236      !
237      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
238      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
239      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
240      REAL(wp) ::   zbtr, ztra               ! local scalars
241      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
242      !!----------------------------------------------------------------------
243      !
244      !                                                          ! ===========
245      DO jn = 1, kjpt                                            ! tracer loop
246         !                                                       ! ===========
247         DO jj = 1, jpjm1
248            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
249               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
250                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
251                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
252                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
253                  zu_bbl = ABS( utr_bbl(ji,jj) )
254                  !
255                  !                                               ! up  -slope T-point (shelf bottom point)
256                  zbtr = r1_e1e2t(iis,jj) / e3t(iis,jj,ikus,Kmm)
257                  ztra = zu_bbl * ( pt(iid,jj,ikus,jn) - pt(iis,jj,ikus,jn) ) * zbtr
258                  pt_rhs(iis,jj,ikus,jn) = pt_rhs(iis,jj,ikus,jn) + ztra
259                  !
260                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
261                     zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,jk,Kmm)
262                     ztra = zu_bbl * ( pt(iid,jj,jk+1,jn) - pt(iid,jj,jk,jn) ) * zbtr
263                     pt_rhs(iid,jj,jk,jn) = pt_rhs(iid,jj,jk,jn) + ztra
264                  END DO
265                  !
266                  zbtr = r1_e1e2t(iid,jj) / e3t(iid,jj,ikud,Kmm)
267                  ztra = zu_bbl * ( pt(iis,jj,ikus,jn) - pt(iid,jj,ikud,jn) ) * zbtr
268                  pt_rhs(iid,jj,ikud,jn) = pt_rhs(iid,jj,ikud,jn) + ztra
269               ENDIF
270               !
271               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
272                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
273                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
274                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
275                  zv_bbl = ABS( vtr_bbl(ji,jj) )
276                  !
277                  ! up  -slope T-point (shelf bottom point)
278                  zbtr = r1_e1e2t(ji,ijs) / e3t(ji,ijs,ikvs,Kmm)
279                  ztra = zv_bbl * ( pt(ji,ijd,ikvs,jn) - pt(ji,ijs,ikvs,jn) ) * zbtr
280                  pt_rhs(ji,ijs,ikvs,jn) = pt_rhs(ji,ijs,ikvs,jn) + ztra
281                  !
282                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
283                     zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,jk,Kmm)
284                     ztra = zv_bbl * ( pt(ji,ijd,jk+1,jn) - pt(ji,ijd,jk,jn) ) * zbtr
285                     pt_rhs(ji,ijd,jk,jn) = pt_rhs(ji,ijd,jk,jn)  + ztra
286                  END DO
287                  !                                               ! down-slope T-point (deep bottom point)
288                  zbtr = r1_e1e2t(ji,ijd) / e3t(ji,ijd,ikvd,Kmm)
289                  ztra = zv_bbl * ( pt(ji,ijs,ikvs,jn) - pt(ji,ijd,ikvd,jn) ) * zbtr
290                  pt_rhs(ji,ijd,ikvd,jn) = pt_rhs(ji,ijd,ikvd,jn) + ztra
291               ENDIF
292            END DO
293            !
294         END DO
295         !                                                  ! ===========
296      END DO                                                ! end tracer
297      !                                                     ! ===========
298   END SUBROUTINE tra_bbl_adv
299
300
301   SUBROUTINE bbl( kt, kit000, cdtype, Kbb, Kmm )
302      !!----------------------------------------------------------------------
303      !!                  ***  ROUTINE bbl  ***
304      !!
305      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
306      !!                advection terms.
307      !!
308      !! ** Method  : * diffusive bbl (nn_bbl_ldf=1) :
309      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
310      !!      along bottom slope gradient) an additional lateral 2nd order
311      !!      diffusion along the bottom slope is added to the general
312      !!      tracer trend, otherwise the additional trend is set to 0.
313      !!      A typical value of ahbt is 2000 m2/s (equivalent to
314      !!      a downslope velocity of 20 cm/s if the condition for slope
315      !!      convection is satified)
316      !!              * advective bbl (nn_bbl_adv=1 or 2) :
317      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
318      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
319      !!        i.e. transport proportional to the along-slope density gradient
320      !!
321      !!      NB: the along slope density gradient is evaluated using the
322      !!      local density (i.e. referenced at a common local depth).
323      !!
324      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
325      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
326      !!----------------------------------------------------------------------
327      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
328      INTEGER         , INTENT(in   ) ::   kit000   ! first time step index
329      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
330      INTEGER         , INTENT(in   ) ::   Kbb, Kmm ! ocean time level index
331      !
332      INTEGER  ::   ji, jj                    ! dummy loop indices
333      INTEGER  ::   ik                        ! local integers
334      INTEGER  ::   iis, iid, ikus, ikud      !   -       -
335      INTEGER  ::   ijs, ijd, ikvs, ikvd      !   -       -
336      REAL(wp) ::   za, zb, zgdrho            ! local scalars
337      REAL(wp) ::   zsign, zsigna, zgbbl      !   -      -
338      REAL(wp), DIMENSION(jpi,jpj,jpts)   :: zts, zab         ! 3D workspace
339      REAL(wp), DIMENSION(jpi,jpj)        :: zub, zvb, zdep   ! 2D workspace
340      !!----------------------------------------------------------------------
341      !
342      IF( kt == kit000 )  THEN
343         IF(lwp)  WRITE(numout,*)
344         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
345         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
346      ENDIF
347      !                                        !* bottom variables (T, S, alpha, beta, depth, velocity)
348      DO jj = 1, jpj
349         DO ji = 1, jpi
350            ik = mbkt(ji,jj)                             ! bottom T-level index
351            zts (ji,jj,jp_tem) = ts(ji,jj,ik,jp_tem,Kbb) ! bottom before T and S
352            zts (ji,jj,jp_sal) = ts(ji,jj,ik,jp_sal,Kbb)
353            !
354            zdep(ji,jj) = gdept(ji,jj,ik,Kmm)            ! bottom T-level reference depth
355            zub (ji,jj) = uu(ji,jj,mbku(ji,jj),Kmm)      ! bottom velocity
356            zvb (ji,jj) = vv(ji,jj,mbkv(ji,jj),Kmm)
357         END DO
358      END DO
359      !
360      CALL eos_rab( zts, zdep, zab, Kmm )
361      !
362      !                                   !-------------------!
363      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
364         !                                !-------------------!
365         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
366            DO ji = 1, fs_jpim1   ! vector opt.
367               !                                                   ! i-direction
368               za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at u-point
369               zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
370               !                                                         ! 2*masked bottom density gradient
371               zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
372                  &      - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
373               !
374               zsign  = SIGN(  0.5, -zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
375               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)       ! masked diffusive flux coeff.
376               !
377               !                                                   ! j-direction
378               za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)              ! 2*(alpha,beta) at v-point
379               zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
380               !                                                         ! 2*masked bottom density gradient
381               zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
382                  &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
383               !
384               zsign = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
385               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
386            END DO
387         END DO
388         !
389      ENDIF
390      !
391      !                                   !-------------------!
392      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
393         !                                !-------------------!
394         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
395         !
396         CASE( 1 )                                   != use of upper velocity
397            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
398               DO ji = 1, fs_jpim1   ! vector opt.
399                  !                                                  ! i-direction
400                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
401                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
402                  !                                                          ! 2*masked bottom density gradient
403                  zgdrho = (  za * ( zts(ji+1,jj,jp_tem) - zts(ji,jj,jp_tem) )    &
404                            - zb * ( zts(ji+1,jj,jp_sal) - zts(ji,jj,jp_sal) )  ) * umask(ji,jj,1)
405                  !
406                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )   ! sign of i-gradient * i-slope
407                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )   ! sign of u * i-slope
408                  !
409                  !                                                          ! bbl velocity
410                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
411                  !
412                  !                                                  ! j-direction
413                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
414                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
415                  !                                                          ! 2*masked bottom density gradient
416                  zgdrho = (  za * ( zts(ji,jj+1,jp_tem) - zts(ji,jj,jp_tem) )    &
417                     &      - zb * ( zts(ji,jj+1,jp_sal) - zts(ji,jj,jp_sal) )  ) * vmask(ji,jj,1)
418                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )   ! sign of j-gradient * j-slope
419                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )   ! sign of u * i-slope
420                  !
421                  !                                                          ! bbl transport
422                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
423               END DO
424            END DO
425            !
426         CASE( 2 )                                 != bbl velocity = F( delta rho )
427            zgbbl = grav * rn_gambbl
428            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
429               DO ji = 1, fs_jpim1   ! vector opt.
430                  !                                                  ! i-direction
431                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
432                  iid  = ji + MAX( 0, mgrhu(ji,jj) )
433                  iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
434                  !
435                  ikud = mbku_d(ji,jj)
436                  ikus = mbku(ji,jj)
437                  !
438                  za = zab(ji+1,jj,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at u-point
439                  zb = zab(ji+1,jj,jp_sal) + zab(ji,jj,jp_sal)
440                  !                                                          !   masked bottom density gradient
441                  zgdrho = 0.5 * (  za * ( zts(iid,jj,jp_tem) - zts(iis,jj,jp_tem) )    &
442                     &            - zb * ( zts(iid,jj,jp_sal) - zts(iis,jj,jp_sal) )  ) * umask(ji,jj,1)
443                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
444                  !
445                  !                                                          ! bbl transport (down-slope direction)
446                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
447                  !
448                  !                                                  ! j-direction
449                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
450                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )
451                  ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
452                  !
453                  ikvd = mbkv_d(ji,jj)
454                  ikvs = mbkv(ji,jj)
455                  !
456                  za = zab(ji,jj+1,jp_tem) + zab(ji,jj,jp_tem)               ! 2*(alpha,beta) at v-point
457                  zb = zab(ji,jj+1,jp_sal) + zab(ji,jj,jp_sal)
458                  !                                                          !   masked bottom density gradient
459                  zgdrho = 0.5 * (  za * ( zts(ji,ijd,jp_tem) - zts(ji,ijs,jp_tem) )    &
460                     &            - zb * ( zts(ji,ijd,jp_sal) - zts(ji,ijs,jp_sal) )  ) * vmask(ji,jj,1)
461                  zgdrho = MAX( 0.e0, zgdrho )                               ! only if shelf is denser than deep
462                  !
463                  !                                                          ! bbl transport (down-slope direction)
464                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
465               END DO
466            END DO
467         END SELECT
468         !
469      ENDIF
470      !
471   END SUBROUTINE bbl
472
473
474   SUBROUTINE tra_bbl_init
475      !!----------------------------------------------------------------------
476      !!                  ***  ROUTINE tra_bbl_init  ***
477      !!
478      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
479      !!
480      !! ** Method  :   Read the nambbl namelist and check the parameters
481      !!              called by nemo_init at the first timestep (kit000)
482      !!----------------------------------------------------------------------
483      INTEGER ::   ji, jj                      ! dummy loop indices
484      INTEGER ::   ii0, ii1, ij0, ij1, ios     ! local integer
485      REAL(wp), DIMENSION(jpi,jpj) ::   zmbku, zmbkv   ! workspace
486      !!
487      NAMELIST/nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
488      !!----------------------------------------------------------------------
489      !
490      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
491      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
492901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambbl in reference namelist' )
493      !
494      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
495      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
496902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambbl in configuration namelist' )
497      IF(lwm) WRITE ( numond, nambbl )
498      !
499      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
500      !
501      IF(lwp) THEN                   !* Parameter control and print
502         WRITE(numout,*)
503         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
504         WRITE(numout,*) '~~~~~~~~~~~~'
505         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
506         WRITE(numout,*) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
507      ENDIF
508      IF( .NOT.ln_trabbl )   RETURN
509      !
510      IF(lwp) THEN
511         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
512         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
513         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
514         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
515      ENDIF
516      !
517      !                              ! allocate trabbl arrays
518      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
519      !
520      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
521      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
522      !
523      !                             !* vertical index of  "deep" bottom u- and v-points
524      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
525         DO ji = 1, jpim1
526            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
527            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
528         END DO
529      END DO
530      ! converte into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
531      zmbku(:,:) = REAL( mbku_d(:,:), wp )   ;     zmbkv(:,:) = REAL( mbkv_d(:,:), wp ) 
532      CALL lbc_lnk_multi( 'trabbl', zmbku,'U',1., zmbkv,'V',1.) 
533      mbku_d(:,:) = MAX( INT( zmbku(:,:) ), 1 ) ;  mbkv_d(:,:) = MAX( NINT( zmbkv(:,:) ), 1 )
534      !
535      !                             !* sign of grad(H) at u- and v-points; zero if grad(H) = 0
536      mgrhu(:,:) = 0   ;   mgrhv(:,:) = 0
537      DO jj = 1, jpjm1
538         DO ji = 1, jpim1
539            IF( gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
540               mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
541            ENDIF
542            !
543            IF( gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) /= 0._wp ) THEN
544               mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
545            ENDIF
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_multi( 'trabbl', e3u_bbl_0, 'U', 1. , 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   END SUBROUTINE tra_bbl_init
562
563   !!======================================================================
564END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.