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 @ 10946

Last change on this file since 10946 was 10946, checked in by acc, 5 years ago

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