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/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA – NEMO

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2106

Last change on this file since 2106 was 2106, checked in by cetlod, 14 years ago

improve trabbl routine according to review

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