New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
trabbl.F90 in branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

source: branches/2014/dev_CNRS0_NOC1_LDF/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 4616

Last change on this file since 4616 was 4616, checked in by gm, 10 years ago

#1260 : see the associated wiki page for explanation

  • Property svn:keywords set to Id
File size: 37.1 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   !!----------------------------------------------------------------------
15#if   defined key_trabbl   ||   defined key_esopa
16   !!----------------------------------------------------------------------
17   !!   'key_trabbl'   or                             bottom boundary layer
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 trdmod_oce     ! trends: ocean variables
31   USE trdtra         ! trends: active tracers
32   USE iom            ! IOM server
33   USE in_out_manager ! I/O manager
34   USE lbclnk         ! ocean lateral boundary conditions
35   USE prtctl         ! Print control
36   USE wrk_nemo       ! Memory Allocation
37   USE timing         ! Timing
38
39   IMPLICIT NONE
40   PRIVATE
41
42   PUBLIC   tra_bbl       !  routine called by step.F90
43   PUBLIC   tra_bbl_init  !  routine called by opa.F90
44   PUBLIC   tra_bbl_dif   !  routine called by trcbbl.F90
45   PUBLIC   tra_bbl_adv   !  -          -          -              -
46   PUBLIC   bbl           !  routine called by trcbbl.F90 and dtadyn.F90
47
48   !                                  !!* Namelist nambbl *
49   LOGICAL , PUBLIC ::   ln_trabbl     !: bottom boundary layer flag
50   INTEGER , PUBLIC ::   nn_bbl_ldf    !: =1   : diffusive bbl or not (=0)
51   INTEGER , PUBLIC ::   nn_bbl_adv    !: =1/2 : advective bbl or not (=0)
52   !                                   !         =1 : advective bbl using the bottom ocean velocity
53   !                                   !         =2 :     -      -  using utr_bbl proportional to grad(rho)
54   REAL(wp), PUBLIC ::   rn_ahtbbl     !: along slope bbl diffusive coefficient [m2/s]
55   REAL(wp), PUBLIC ::   rn_gambbl     !: lateral coeff. for bottom boundary layer scheme [s]
56
57   LOGICAL , PUBLIC ::   l_bbl         !: flag to compute bbl diffu. flux coef and transport
58
59   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   utr_bbl  , vtr_bbl   ! u- (v-) transport in the bottom boundary layer
60   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   ahu_bbl  , ahv_bbl   ! masked diffusive bbl coeff. at u & v-pts
61
62   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mbku_d   , mbkv_d      ! vertical index of the "lower" bottom ocean U/V-level (PUBLIC for TAM)
63   INTEGER , ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC ::   mgrhu    , mgrhv       ! = +/-1, sign of grad(H) in u-(v-)direction (PUBLIC for TAM)
64   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)         ::   ahu_bbl_0, ahv_bbl_0   ! diffusive bbl flux coefficients at u and v-points
65   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)
66
67   !! * Substitutions
68#  include "domzgr_substitute.h90"
69#  include "vectopt_loop_substitute.h90"
70   !!----------------------------------------------------------------------
71   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
72   !! $Id$
73   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
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      IF( lk_mpp            )   CALL mpp_sum ( 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 )
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      !!
105      REAL(wp), POINTER, DIMENSION(:,:,:) ::  ztrdt, ztrds
106      !!----------------------------------------------------------------------
107      !
108      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl')
109      !
110      IF( l_trdtra )   THEN                        !* Save ta and sa trends
111         CALL wrk_alloc( jpi, jpj, jpk, ztrdt, ztrds )
112         ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
113         ztrds(:,:,:) = tsa(:,:,:,jp_sal)
114      ENDIF
115
116      IF( l_bbl )  CALL bbl( kt, nit000, 'TRA' )   !* bbl coef. and transport (only if not already done in trcbbl)
117
118      IF( nn_bbl_ldf == 1 ) THEN                   !* Diffusive bbl
119         !
120         CALL tra_bbl_dif( tsb, tsa, jpts )
121         IF( ln_ctl )  &
122         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_ldf  - Ta: ', mask1=tmask, &
123         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
124         ! lateral boundary conditions ; just need for outputs
125         CALL lbc_lnk( ahu_bbl, 'U', 1. )     ;     CALL lbc_lnk( ahv_bbl, 'V', 1. )
126         CALL iom_put( "ahu_bbl", ahu_bbl )   ! bbl diffusive flux i-coef
127         CALL iom_put( "ahv_bbl", ahv_bbl )   ! bbl diffusive flux j-coef
128         !
129      END IF
130
131      IF( nn_bbl_adv /= 0 ) THEN                !* Advective bbl
132         !
133         CALL tra_bbl_adv( tsb, tsa, jpts )
134         IF(ln_ctl)   &
135         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' bbl_adv  - Ta: ', mask1=tmask,   &
136         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=           ' Sa: ', mask2=tmask, clinfo3='tra' )
137         ! lateral boundary conditions ; just need for outputs
138         CALL lbc_lnk( utr_bbl, 'U', 1. )     ;   CALL lbc_lnk( vtr_bbl, 'V', 1. )
139         CALL iom_put( "uoce_bbl", utr_bbl )  ! bbl i-transport
140         CALL iom_put( "voce_bbl", vtr_bbl )  ! bbl j-transport
141         !
142      END IF
143
144      IF( l_trdtra )   THEN                      ! save the horizontal diffusive trends for further diagnostics
145         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
146         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
147         CALL trd_tra( kt, 'TRA', jp_tem, jptra_trd_bbl, ztrdt )
148         CALL trd_tra( kt, 'TRA', jp_sal, jptra_trd_bbl, ztrds )
149         CALL wrk_dealloc( jpi, jpj, jpk, ztrdt, ztrds )
150      ENDIF
151      !
152      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl')
153      !
154   END SUBROUTINE tra_bbl
155
156
157   SUBROUTINE tra_bbl_dif( ptb, pta, kjpt )
158      !!----------------------------------------------------------------------
159      !!                  ***  ROUTINE tra_bbl_dif  ***
160      !!
161      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
162      !!                advection terms.
163      !!
164      !! ** Method  :
165      !!        * diffusive bbl (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), POINTER, DIMENSION(:,:) :: zptb
186      !!----------------------------------------------------------------------
187      !
188      IF( nn_timing == 1 )  CALL timing_start('tra_bbl_dif')
189      !
190      CALL wrk_alloc( jpi, jpj, zptb )
191      !
192      DO jn = 1, kjpt                                     ! tracer loop
193         !                                                ! ===========
194         DO jj = 1, jpj
195            DO ji = 1, jpi
196               ik = mbkt(ji,jj)                        ! bottom T-level index
197               zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S
198            END DO
199         END DO
200         !                                                ! Compute the trend
201         DO jj = 2, jpjm1
202            DO ji = 2, jpim1
203               ik = mbkt(ji,jj)                            ! bottom T-level index
204               pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                  &
205                  &             + (  ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )     &
206                  &                - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )     &
207                  &                + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )     &
208                  &                - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )  )  &
209                  &             / ( e1e2t(ji,jj) * fse3t(ji,jj,ik) )
210            END DO
211         END DO
212         !                                                  ! ===========
213      END DO                                                ! end tracer
214      !                                                     ! ===========
215      CALL wrk_dealloc( jpi, jpj, zptb )
216      !
217      IF( nn_timing == 1 )  CALL timing_stop('tra_bbl_dif')
218      !
219   END SUBROUTINE tra_bbl_dif
220
221
222   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
223      !!----------------------------------------------------------------------
224      !!                  ***  ROUTINE trc_bbl  ***
225      !!
226      !! ** Purpose :   Compute the before passive tracer trend associated
227      !!     with the bottom boundary layer and add it to the general trend
228      !!     of tracer equations.
229      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
230      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
231      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
232      !!                       transport proportional to the along-slope density gradient
233      !!
234      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
235      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
236      !!----------------------------------------------------------------------
237      INTEGER                              , INTENT(in   ) ::   kjpt   ! number of tracers
238      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb    ! before and now tracer fields
239      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta    ! tracer trend
240      !
241      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
242      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
243      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
244      REAL(wp) ::   zbtr, ztra               ! local scalars
245      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
246      !!----------------------------------------------------------------------
247      !
248      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_adv')
249      !                                                          ! ===========
250      DO jn = 1, kjpt                                            ! tracer loop
251         !                                                       ! ===========
252         DO jj = 1, jpjm1
253            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
254               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
255                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
256                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
257                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
258                  zu_bbl = ABS( utr_bbl(ji,jj) )
259                  !
260                  !                                               ! up  -slope T-point (shelf bottom point)
261                  zbtr = r1_e1e2t(iis,jj) / fse3t(iis,jj,ikus)
262                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
263                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
264                  !
265                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
266                     zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,jk)
267                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
268                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
269                  END DO
270                  !
271                  zbtr = r1_e1e2t(iid,jj) / fse3t(iid,jj,ikud)
272                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
273                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
274               ENDIF
275               !
276               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
277                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
278                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
279                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
280                  zv_bbl = ABS( vtr_bbl(ji,jj) )
281                  !
282                  ! up  -slope T-point (shelf bottom point)
283                  zbtr = r1_e1e2t(ji,ijs) / fse3t(ji,ijs,ikvs)
284                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
285                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
286                  !
287                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
288                     zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,jk)
289                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
290                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
291                  END DO
292                  !                                               ! down-slope T-point (deep bottom point)
293                  zbtr = r1_e1e2t(ji,ijd) / fse3t(ji,ijd,ikvd)
294                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
295                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
296               ENDIF
297            END DO
298            !
299         END DO
300         !                                                  ! ===========
301      END DO                                                ! end tracer
302      !                                                     ! ===========
303      !
304      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_adv')
305      !
306   END SUBROUTINE tra_bbl_adv
307
308
309   SUBROUTINE bbl( kt, kit000, cdtype )
310      !!----------------------------------------------------------------------
311      !!                  ***  ROUTINE bbl  ***
312      !!
313      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
314      !!                advection terms.
315      !!
316      !! ** Method  :
317      !!        * diffusive bbl (nn_bbl_ldf=1) :
318      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
319      !!      along bottom slope gradient) an additional lateral 2nd order
320      !!      diffusion along the bottom slope is added to the general
321      !!      tracer trend, otherwise the additional trend is set to 0.
322      !!      A typical value of ahbt is 2000 m2/s (equivalent to
323      !!      a downslope velocity of 20 cm/s if the condition for slope
324      !!      convection is satified)
325      !!        * advective bbl (nn_bbl_adv=1 or 2) :
326      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
327      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
328      !!        i.e. transport proportional to the along-slope density gradient
329      !!
330      !!      NB: the along slope density gradient is evaluated using the
331      !!      local density (i.e. referenced at a common local depth).
332      !!
333      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
334      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
335      !!----------------------------------------------------------------------
336      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
337      INTEGER         , INTENT(in   ) ::   kit000          ! first time step index
338      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
339      !!
340      INTEGER  ::   ji, jj                    ! dummy loop indices
341      INTEGER  ::   ik                        ! local integers
342      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
343      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
344      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
345      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
346      !!
347      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
348      REAL(wp), POINTER, DIMENSION(:,:) :: zub, zvb, ztb, zsb, zdep
349      !!----------------------- zv_bbl-----------------------------------------------
350      ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
351      ! ================            pft :  potential temperature in degrees celcius
352      !                             pfs :  salinity anomaly (s-35) in psu
353      !                             pfh :  depth in meters
354      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
355      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
356         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
357                                   - 0.203814e-03 ) * pft                    &
358                                   + 0.170907e-01 ) * pft                    &
359                                   + 0.665157e-01                            &
360         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
361         +  ( ( - 0.302285e-13 * pfh                                         &
362                - 0.251520e-11 * pfs                                         &
363                + 0.512857e-12 * pft * pft          ) * pfh                  &
364                                     - 0.164759e-06   * pfs                  &
365             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
366                                     + 0.380374e-04 ) * pfh
367      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
368         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
369                                 - 0.301985e-05 ) * pft                      &
370                                 + 0.785567e-03                              &
371         + (     0.515032e-08 * pfs                                          &
372               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
373               +(  (   0.121551e-17 * pfh                                    &
374                     - 0.602281e-15 * pfs                                    &
375                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
376                                          + 0.408195e-10   * pfs             &
377                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
378                                          - 0.121555e-07 ) * pfh
379      !!----------------------------------------------------------------------
380
381      !
382      IF( nn_timing == 1 )  CALL timing_start( 'bbl')
383      !
384      CALL wrk_alloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )
385      !
386
387      IF( kt == kit000 )  THEN
388         IF(lwp)  WRITE(numout,*)
389         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
390         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
391      ENDIF
392
393      !                                        !* bottom temperature, salinity, velocity and depth
394      DO jj = 1, jpj
395         DO ji = 1, jpi
396            ik = mbkt(ji,jj)                        ! bottom T-level index
397            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem) * tmask(ji,jj,1)      ! bottom before T and S
398            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal) * tmask(ji,jj,1)
399            zdep(ji,jj) = gdept_0(ji,jj,ik)         ! bottom T-level reference depth
400            !
401            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
402            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
403         END DO
404      END DO
405
406      !                                   !-------------------!
407      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
408         !                                !-------------------!
409         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
410            DO ji = 1, jpim1
411               !                                                ! i-direction
412               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
413               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
414               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
415               !                                                         ! masked bbl i-gradient of density
416               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
417                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
418               !
419               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
420               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff.
421               !
422               !                                                ! j-direction
423               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
424               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
425               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
426               !                                                         ! masked bbl j-gradient of density
427               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
428                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
429               !
430               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
431               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
432               !
433            END DO
434         END DO
435         !
436      ENDIF
437
438      !                                   !-------------------!
439      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
440         !                                !-------------------!
441         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
442         !
443         CASE( 1 )                                   != use of upper velocity
444            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
445               DO ji = 1, fs_jpim1   ! vector opt.
446                  !                                               ! i-direction
447                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
448                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
449                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
450                  !                                                           ! masked bbl i-gradient of density
451                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
452                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
453                  !
454                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
455                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
456                  !
457                  !                                                           ! bbl velocity
458                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
459                  !
460                  !                                               ! j-direction
461                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
462                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
463                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
464                  !                                                           ! masked bbl j-gradient of density
465                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
466                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
467                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
468                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
469                  !
470                  !                                                           ! bbl velocity
471                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
472               END DO
473            END DO
474            !
475         CASE( 2 )                                 != bbl velocity = F( delta rho )
476            zgbbl = grav * rn_gambbl
477            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
478               DO ji = 1, fs_jpim1   ! vector opt.
479                  !                                         ! i-direction
480                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
481                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
482                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
483                  !
484                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
485                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
486                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
487                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
488                  zgdrho =    fsbeta( zt, zs, zh )                                    &
489                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    &
490                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
491                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
492                  !
493                  !                                             ! bbl transport (down-slope direction)
494                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
495                  !
496                  !                                         ! j-direction
497                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
498                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
499                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
500                  !
501                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
502                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
503                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
504                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
505                  zgdrho =    fsbeta( zt, zs, zh )                                    &
506                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    &
507                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
508                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
509                  !
510                  !                                             ! bbl transport (down-slope direction)
511                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
512               END DO
513            END DO
514         END SELECT
515         !
516      ENDIF
517      !
518      CALL wrk_dealloc( jpi, jpj, zub, zvb, ztb, zsb, zdep )
519      !
520      IF( nn_timing == 1 )  CALL timing_stop( 'bbl')
521      !
522   END SUBROUTINE bbl
523
524
525   SUBROUTINE tra_bbl_init
526      !!----------------------------------------------------------------------
527      !!                  ***  ROUTINE tra_bbl_init  ***
528      !!
529      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
530      !!
531      !! ** Method  :   Read the nambbl namelist and check the parameters
532      !!              called by nemo_init at the first timestep (kit000)
533      !!----------------------------------------------------------------------
534      INTEGER ::   ji, jj               ! dummy loop indices
535      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
536      INTEGER  ::   ios                 ! Local integer output status for namelist read
537      REAL(wp), POINTER, DIMENSION(:,:) :: zmbk
538      !!
539      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
540      !!----------------------------------------------------------------------
541      !
542      IF( nn_timing == 1 )  CALL timing_start( 'tra_bbl_init')
543      !
544      CALL wrk_alloc( jpi, jpj, zmbk )
545      !
546
547      REWIND( numnam_ref )              ! Namelist nambbl in reference namelist : Bottom boundary layer scheme
548      READ  ( numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
549901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in reference namelist', lwp )
550
551      REWIND( numnam_cfg )              ! Namelist nambbl in configuration namelist : Bottom boundary layer scheme
552      READ  ( numnam_cfg, nambbl, IOSTAT = ios, ERR = 902 )
553902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambbl in configuration namelist', lwp )
554      WRITE ( numond, nambbl )
555      !
556      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
557      !
558      IF(lwp) THEN                   !* Parameter control and print
559         WRITE(numout,*)
560         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
561         WRITE(numout,*) '~~~~~~~~~~~~'
562         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
563         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
564         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
565         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
566         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
567      ENDIF
568
569      !                              ! allocate trabbl arrays
570      IF( tra_bbl_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'tra_bbl_init : unable to allocate arrays' )
571
572      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
573      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
574
575      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
576
577      !                             !* vertical index of  "deep" bottom u- and v-points
578      DO jj = 1, jpjm1                    ! (the "shelf" bottom k-indices are mbku and mbkv)
579         DO ji = 1, jpim1
580            mbku_d(ji,jj) = MAX(  mbkt(ji+1,jj  ) , mbkt(ji,jj)  )   ! >= 1 as mbkt=1 over land
581            mbkv_d(ji,jj) = MAX(  mbkt(ji  ,jj+1) , mbkt(ji,jj)  )
582         END DO
583      END DO
584      ! convert into REAL to use lbc_lnk ; impose a min value of 1 as a zero can be set in lbclnk
585      zmbk(:,:) = REAL( mbku_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
586      zmbk(:,:) = REAL( mbkv_d(:,:), wp )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
587
588                                     !* sign of grad(H) at u- and v-points
589      mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0.
590      DO jj = 1, jpjm1
591         DO ji = 1, jpim1
592            mgrhu(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji+1,jj,mbkt(ji+1,jj)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
593            mgrhv(ji,jj) = INT(  SIGN( 1.e0, gdept_0(ji,jj+1,mbkt(ji,jj+1)) - gdept_0(ji,jj,mbkt(ji,jj)) )  )
594         END DO
595      END DO
596
597      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
598         DO ji = 1, jpim1           ! minimum of top & bottom e3u_0 (e3v_0)
599            e3u_bbl_0(ji,jj) = MIN( e3u_0(ji,jj,mbkt(ji+1,jj  )), e3u_0(ji,jj,mbkt(ji,jj)) )
600            e3v_bbl_0(ji,jj) = MIN( e3v_0(ji,jj,mbkt(ji  ,jj+1)), e3v_0(ji,jj,mbkt(ji,jj)) )
601         END DO
602      END DO
603      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
604
605      !                             !* masked diffusive flux coefficients
606      ahu_bbl_0(:,:) = rn_ahtbbl * e2_e1u(:,:) * e3u_bbl_0(:,:) * umask(:,:,1)
607      ahv_bbl_0(:,:) = rn_ahtbbl * e1_e2v(:,:) * e3v_bbl_0(:,:) * vmask(:,:,1)
608
609
610      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
611         !
612         SELECT CASE ( jp_cfg )
613         CASE ( 2 )                          ! ORCA_R2
614            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
615            ii0 = 139   ;   ii1 = 140
616            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
617            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
618            !
619            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
620            ii0 = 161   ;   ii1 = 162
621            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
622            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
623            !
624         CASE ( 4 )                          ! ORCA_R4
625            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
626            ii0 =  70   ;   ii1 =  71
627            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
628            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
629         END SELECT
630         !
631      ENDIF
632      !
633      CALL wrk_dealloc( jpi, jpj, zmbk )
634      !
635      IF( nn_timing == 1 )  CALL timing_stop( 'tra_bbl_init')
636      !
637   END SUBROUTINE tra_bbl_init
638
639#else
640   !!----------------------------------------------------------------------
641   !!   Dummy module :                      No bottom boundary layer scheme
642   !!----------------------------------------------------------------------
643   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
644CONTAINS
645   SUBROUTINE tra_bbl_init               ! Dummy routine
646   END SUBROUTINE tra_bbl_init
647   SUBROUTINE tra_bbl( kt )              ! Dummy routine
648      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
649   END SUBROUTINE tra_bbl
650#endif
651
652   !!======================================================================
653END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.