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

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2528

Last change on this file since 2528 was 2528, checked in by rblod, 13 years ago

Update NEMOGCM from branch nemo_v3_3_beta

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