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

source: branches/2010_and_older/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/trabbl.F90 @ 8246

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

bug correction in trabbl.F90

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 36.4 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 in 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( ptb, pta, 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  :   pta   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   ) ::   ptb   ! before and now tracer fields
164      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta   ! tracer trend
165      !!
166      INTEGER  ::   ji, jj, jn   ! dummy loop indices
167      INTEGER  ::   ik           ! local integers
168      REAL(wp) ::   zbtr         ! local scalars
169      REAL(wp), DIMENSION(jpi,jpj) ::   zptb   ! temporary array
170      !!----------------------------------------------------------------------
171      !
172      !                                                   ! ===========
173      DO jn = 1, kjpt                                     ! tracer loop
174         !                                                ! ===========
175#  if defined key_vectopt_loop
176         DO jj = 1, 1   ! vector opt. (forced unrolling)
177            DO ji = 1, jpij
178#else
179         DO jj = 1, jpj
180            DO ji = 1, jpi
181#endif
182               ik = mbkt(ji,jj)                        ! bottom T-level index
183               zptb(ji,jj) = ptb(ji,jj,ik,jn)              ! bottom before T and S
184            END DO
185         END DO
186         !                                                ! Compute the trend
187#  if defined key_vectopt_loop
188         DO jj = 1, 1   ! vector opt. (forced unrolling)
189            DO ji = jpi+1, jpij-jpi-1
190#  else
191         DO jj = 2, jpjm1
192            DO ji = 2, jpim1
193#  endif
194                ik = mbkt(ji,jj)                            ! bottom T-level index
195                zbtr = e1e2t_r(ji,jj)  / fse3t(ji,jj,ik)
196                pta(ji,jj,ik,jn) = pta(ji,jj,ik,jn)                                                         &
197                  &               + (   ahu_bbl(ji  ,jj  ) * ( zptb(ji+1,jj  ) - zptb(ji  ,jj  ) )   &
198                  &                   - ahu_bbl(ji-1,jj  ) * ( zptb(ji  ,jj  ) - zptb(ji-1,jj  ) )   &
199                  &                   + ahv_bbl(ji  ,jj  ) * ( zptb(ji  ,jj+1) - zptb(ji  ,jj  ) )   &
200                  &                   - ahv_bbl(ji  ,jj-1) * ( zptb(ji  ,jj  ) - zptb(ji  ,jj-1) )   ) * zbtr
201            END DO
202         END DO
203         !                                                  ! ===========
204      END DO                                                ! end tracer
205      !                                                     ! ===========
206   END SUBROUTINE tra_bbl_dif
207   
208
209   SUBROUTINE tra_bbl_adv( ptb, pta, kjpt )
210      !!----------------------------------------------------------------------
211      !!                  ***  ROUTINE trc_bbl  ***
212      !!
213      !! ** Purpose :   Compute the before passive tracer trend associated
214      !!     with the bottom boundary layer and add it to the general trend
215      !!     of tracer equations.
216      !! ** Method  :   advective bbl (nn_bbl_adv = 1 or 2) :
217      !!      nn_bbl_adv = 1   use of the ocean near bottom velocity as bbl velocity
218      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation i.e.
219      !!                       transport proportional to the along-slope density gradient                   
220      !!
221      !!
222      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
223      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
224      !!
225      !!---------------------------------------------------------------------- 
226      INTEGER                              , INTENT(in   ) ::   kjpt    ! number of tracers
227      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in   ) ::   ptb   ! before and now tracer fields
228      REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) ::   pta   ! tracer trend
229      !!
230      INTEGER  ::   ji, jj, jk, jn           ! dummy loop indices
231      INTEGER  ::   iis , iid , ijs , ijd    ! local integers
232      INTEGER  ::   ikus, ikud, ikvs, ikvd   !   -       -
233      REAL(wp) ::   zbtr, ztra               ! local scalars
234      REAL(wp) ::   zu_bbl, zv_bbl           !   -      -
235      !!----------------------------------------------------------------------
236
237      !                                                          ! ===========
238      DO jn = 1, kjpt                                            ! tracer loop
239         !                                                       ! ===========         
240# if defined key_vectopt_loop
241         DO jj = 1, 1
242            DO ji = 1, jpij-jpi-1   ! vector opt. (forced unrolling)
243# else
244         DO jj = 1, jpjm1
245            DO ji = 1, jpim1            ! CAUTION start from i=1 to update i=2 when cyclic east-west
246# endif
247               IF( utr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero i-direction bbl advection
248                  ! down-slope i/k-indices (deep)      &   up-slope i/k indices (shelf)
249                  iid  = ji + MAX( 0, mgrhu(ji,jj) )   ;   iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
250                  ikud = mbku_d(ji,jj)                 ;   ikus = mbku(ji,jj)
251                  zu_bbl = ABS( utr_bbl(ji,jj) )
252                  !
253                  !                                               ! up  -slope T-point (shelf bottom point)
254                  zbtr = e1e2t_r(iis,jj) / fse3t(iis,jj,ikus)
255                  ztra = zu_bbl * ( ptb(iid,jj,ikus,jn) - ptb(iis,jj,ikus,jn) ) * zbtr
256                  pta(iis,jj,ikus,jn) = pta(iis,jj,ikus,jn) + ztra
257                  !                   
258                  DO jk = ikus, ikud-1                            ! down-slope upper to down T-point (deep column)
259                     zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,jk)
260                     ztra = zu_bbl * ( ptb(iid,jj,jk+1,jn) - ptb(iid,jj,jk,jn) ) * zbtr
261                     pta(iid,jj,jk,jn) = pta(iid,jj,jk,jn) + ztra
262                  END DO
263                  !
264                  zbtr = e1e2t_r(iid,jj) / fse3t(iid,jj,ikud)
265                  ztra = zu_bbl * ( ptb(iis,jj,ikus,jn) - ptb(iid,jj,ikud,jn) ) * zbtr
266                  pta(iid,jj,ikud,jn) = pta(iid,jj,ikud,jn) + ztra
267               ENDIF
268               !
269               IF( vtr_bbl(ji,jj) /= 0.e0 ) THEN            ! non-zero j-direction bbl advection
270                  ! down-slope j/k-indices (deep)        &   up-slope j/k indices (shelf)
271                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )     ;   ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
272                  ikvd = mbkv_d(ji,jj)                   ;   ikvs = mbkv(ji,jj)
273                  zv_bbl = ABS( vtr_bbl(ji,jj) )
274                  !
275                  ! up  -slope T-point (shelf bottom point)
276                  zbtr = e1e2t_r(ji,ijs) / fse3t(ji,ijs,ikvs)
277                  ztra = zv_bbl * ( ptb(ji,ijd,ikvs,jn) - ptb(ji,ijs,ikvs,jn) ) * zbtr
278                  pta(ji,ijs,ikvs,jn) = pta(ji,ijs,ikvs,jn) + ztra
279                  !                   
280                  DO jk = ikvs, ikvd-1                            ! down-slope upper to down T-point (deep column)
281                     zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,jk)
282                     ztra = zv_bbl * ( ptb(ji,ijd,jk+1,jn) - ptb(ji,ijd,jk,jn) ) * zbtr
283                     pta(ji,ijd,jk,jn) = pta(ji,ijd,jk,jn)  + ztra
284                  END DO
285                  !                                               ! down-slope T-point (deep bottom point)
286                  zbtr = e1e2t_r(ji,ijd) / fse3t(ji,ijd,ikvd)
287                  ztra = zv_bbl * ( ptb(ji,ijs,ikvs,jn) - ptb(ji,ijd,ikvd,jn) ) * zbtr
288                  pta(ji,ijd,ikvd,jn) = pta(ji,ijd,ikvd,jn) + ztra
289               ENDIF
290            END DO
291            !
292         END DO
293         !                                                  ! ===========
294      END DO                                                ! end tracer
295      !                                                     ! ===========
296   END SUBROUTINE tra_bbl_adv
297
298
299   SUBROUTINE bbl( kt, cdtype )
300      !!----------------------------------------------------------------------
301      !!                  ***  ROUTINE bbl  ***
302      !!                   
303      !! ** Purpose :   Computes the bottom boundary horizontal and vertical
304      !!                advection terms.
305      !!
306      !! ** Method  :   
307      !!        * diffusive bbl (nn_bbl_ldf=1) :
308      !!        When the product grad( rho) * grad(h) < 0 (where grad is an
309      !!      along bottom slope gradient) an additional lateral 2nd order
310      !!      diffusion along the bottom slope is added to the general
311      !!      tracer trend, otherwise the additional trend is set to 0.
312      !!      A typical value of ahbt is 2000 m2/s (equivalent to
313      !!      a downslope velocity of 20 cm/s if the condition for slope
314      !!      convection is satified)
315      !!        * advective bbl (nn_bbl_adv=1 or 2) :
316      !!      nn_bbl_adv = 1   use of the ocean velocity as bbl velocity
317      !!      nn_bbl_adv = 2   follow Campin and Goosse (1999) implentation
318      !!        i.e. transport proportional to the along-slope density gradient
319      !!
320      !!      NB: the along slope density gradient is evaluated using the
321      !!      local density (i.e. referenced at a common local depth).
322      !!
323      !! References : Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
324      !!              Campin, J.-M., and H. Goosse, 1999, Tellus, 412-430.
325      !!---------------------------------------------------------------------- 
326      INTEGER         , INTENT(in   ) ::   kt       ! ocean time-step index
327      CHARACTER(len=3), INTENT(in   ) ::   cdtype   ! =TRA or TRC (tracer indicator)
328      !!
329      INTEGER  ::   ji, jj                    ! dummy loop indices
330      INTEGER  ::   ik                        ! local integers
331      INTEGER  ::   iis , iid , ijs , ijd     !   -       -
332      INTEGER  ::   ikus, ikud, ikvs, ikvd    !   -       -
333      REAL(wp) ::   zsign, zsigna, zgbbl      ! local scalars
334      REAL(wp) ::   zgdrho, zt, zs, zh        !   -      -
335      REAL(wp), DIMENSION(jpi,jpj) ::   zub, zvb, ztb, zsb, zdep  !  2D workspace
336      !!
337      REAL(wp) ::   fsalbt, fsbeta, pft, pfs, pfh   ! statement function
338      !!----------------------- zv_bbl-----------------------------------------------
339      ! ratio alpha/beta = fsalbt : ratio of thermal over saline expension coefficients
340      ! ================            pft :  potential temperature in degrees celcius
341      !                             pfs :  salinity anomaly (s-35) in psu
342      !                             pfh :  depth in meters
343      ! nn_eos = 0  (Jackett and McDougall 1994 formulation)
344      fsalbt( pft, pfs, pfh ) =                                              &   ! alpha/beta
345         ( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft                    &
346                                   - 0.203814e-03 ) * pft                    &
347                                   + 0.170907e-01 ) * pft                    &
348                                   + 0.665157e-01                            &
349         +(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs   &
350         +  ( ( - 0.302285e-13 * pfh                                         &
351                - 0.251520e-11 * pfs                                         &
352                + 0.512857e-12 * pft * pft          ) * pfh                  &
353                                     - 0.164759e-06   * pfs                  &
354             +(   0.791325e-08 * pft - 0.933746e-06 ) * pft                  &
355                                     + 0.380374e-04 ) * pfh
356      fsbeta( pft, pfs, pfh ) =                                              &   ! beta
357         ( ( -0.415613e-09 * pft + 0.555579e-07 ) * pft                      &
358                                 - 0.301985e-05 ) * pft                      &
359                                 + 0.785567e-03                              &
360         + (     0.515032e-08 * pfs                                          &
361               + 0.788212e-08 * pft - 0.356603e-06 ) * pfs                   &
362               +(  (   0.121551e-17 * pfh                                    &
363                     - 0.602281e-15 * pfs                                    &
364                     - 0.175379e-14 * pft + 0.176621e-12 ) * pfh             &
365                                          + 0.408195e-10   * pfs             &
366                 + ( - 0.213127e-11 * pft + 0.192867e-09 ) * pft             &
367                                          - 0.121555e-07 ) * pfh
368      !!----------------------------------------------------------------------
369     
370      IF( kt == nit000 )  THEN
371         IF(lwp)  WRITE(numout,*)
372         IF(lwp)  WRITE(numout,*) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
373         IF(lwp)  WRITE(numout,*) '~~~~~~~~~~'
374      ENDIF
375     
376      !                                        !* bottom temperature, salinity, velocity and depth
377#if defined key_vectopt_loop
378      DO jj = 1, 1   ! vector opt. (forced unrolling)
379         DO ji = 1, jpij
380#else
381      DO jj = 1, jpj
382         DO ji = 1, jpi
383#endif
384            ik = mbkt(ji,jj)                        ! bottom T-level index
385            ztb (ji,jj) = tsb(ji,jj,ik,jp_tem)      ! bottom before T and S
386            zsb (ji,jj) = tsb(ji,jj,ik,jp_sal)
387            zdep(ji,jj) = fsdept_0(ji,jj,ik)        ! bottom T-level reference depth
388            !
389            zub(ji,jj) = un(ji,jj,mbku(ji,jj))      ! bottom velocity
390            zvb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) 
391         END DO
392      END DO
393     
394      !                                   !-------------------!
395      IF( nn_bbl_ldf == 1 ) THEN          !   diffusive bbl   !
396         !                                !-------------------!
397         DO jj = 1, jpjm1                      ! (criteria for non zero flux: grad(rho).grad(h) < 0 )
398            DO ji = 1, jpim1             
399               !                                                ! i-direction
400               zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )  ! T, S anomalie, and depth
401               zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
402               zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
403               !                                                         ! masked bbl i-gradient of density
404               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    &
405                  &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
406               !                                                     
407               zsign          = SIGN(  0.5, - zgdrho * REAL( mgrhu(ji,jj) )  )    ! sign of ( i-gradient * i-slope )
408               ahu_bbl(ji,jj) = ( 0.5 - zsign ) * ahu_bbl_0(ji,jj)                  ! masked diffusive flux coeff.
409               !
410               !                                                ! j-direction
411               zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                ! T, S anomalie, and depth
412               zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
413               zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
414               !                                                         ! masked bbl j-gradient of density
415               zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    &
416                  &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
417               !                                                   
418               zsign          = SIGN(  0.5, -zgdrho * REAL( mgrhv(ji,jj) )  )     ! sign of ( j-gradient * j-slope )
419               ahv_bbl(ji,jj) = ( 0.5 - zsign ) * ahv_bbl_0(ji,jj)
420               !
421            END DO
422         END DO
423         !
424      ENDIF
425
426      !                                   !-------------------!
427      IF( nn_bbl_adv /= 0 ) THEN          !   advective bbl   !
428         !                                !-------------------!
429         SELECT CASE ( nn_bbl_adv )             !* bbl transport type
430         !
431         CASE( 1 )                                   != use of upper velocity
432            DO jj = 1, jpjm1                                 ! criteria: grad(rho).grad(h)<0  and grad(rho).grad(h)<0
433               DO ji = 1, fs_jpim1   ! vector opt.
434                  !                                               ! i-direction
435                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )                  ! T, S anomalie, and depth
436                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
437                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
438                  !                                                           ! masked bbl i-gradient of density
439                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji+1,jj) - ztb(ji,jj) )    & 
440                     &                             - ( zsb(ji+1,jj) - zsb(ji,jj) )  ) * umask(ji,jj,1)
441                  !                                                         
442                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhu(ji,jj) )  )    ! sign of i-gradient * i-slope
443                  zsigna= SIGN(  0.5, zub(ji,jj) * REAL( mgrhu(ji,jj) )  )    ! sign of u * i-slope
444                  !
445                  !                                                           ! bbl velocity
446                  utr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e2u(ji,jj) * e3u_bbl_0(ji,jj) * zub(ji,jj)
447                  !
448                  !                                               ! j-direction
449                  zt = 0.5 * ( ztb (ji,jj+1) + ztb (ji,jj) )                  ! T, S anomalie, and depth
450                  zs = 0.5 * ( zsb (ji,jj+1) + zsb (ji,jj) ) - 35.0
451                  zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
452                  !                                                           ! masked bbl j-gradient of density
453                  zgdrho = (  fsalbt( zt, zs, zh ) * ( ztb(ji,jj+1) - ztb(ji,jj) )    & 
454                     &                             - ( zsb(ji,jj+1) - zsb(ji,jj) )  ) * vmask(ji,jj,1)
455                  zsign = SIGN(  0.5, - zgdrho   * REAL( mgrhv(ji,jj) )  )    ! sign of j-gradient * j-slope
456                  zsigna= SIGN(  0.5, zvb(ji,jj) * REAL( mgrhv(ji,jj) )  )    ! sign of u * i-slope
457                  !
458                  !                                                           ! bbl velocity
459                  vtr_bbl(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * e1v(ji,jj) * e3v_bbl_0(ji,jj) * zvb(ji,jj)
460               END DO
461            END DO
462            !
463         CASE( 2 )                                 != bbl velocity = F( delta rho )
464            zgbbl = grav * rn_gambbl
465            DO jj = 1, jpjm1                            ! criteria: rho_up > rho_down
466               DO ji = 1, fs_jpim1   ! vector opt.
467                  !                                         ! i-direction
468                  ! down-slope T-point i/k-index (deep)  &   up-slope T-point i/k-index (shelf)
469                  iid  = ji + MAX( 0, mgrhu(ji,jj) )     ;    iis  = ji + 1 - MAX( 0, mgrhu(ji,jj) )
470                  ikud = mbku_d(ji,jj)                   ;    ikus = mbku(ji,jj)
471                  !
472                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
473                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji+1,jj) )           ! mid slope depth of T, S, and depth
474                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji+1,jj) ) - 35.0
475                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
476                  zgdrho =    fsbeta( zt, zs, zh )                                    &
477                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(iid,jj) - ztb(iis,jj) )    & 
478                     &                             - ( zsb(iid,jj) - zsb(iis,jj) )  ) * umask(ji,jj,1)
479                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
480                  !
481                  !                                             ! bbl transport (down-slope direction)
482                  utr_bbl(ji,jj) = e2u(ji,jj) * e3u_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhu(ji,jj) )
483                  !
484                  !                                         ! j-direction
485                  !  down-slope T-point j/k-index (deep)  &   of the up  -slope T-point j/k-index (shelf)
486                  ijd  = jj + MAX( 0, mgrhv(ji,jj) )      ;    ijs  = jj + 1 - MAX( 0, mgrhv(ji,jj) )
487                  ikvd = mbkv_d(ji,jj)                    ;    ikvs = mbkv(ji,jj)
488                  !
489                  !                                             ! mid-depth density anomalie (up-slope minus down-slope)
490                  zt = 0.5 * ( ztb (ji,jj) + ztb (ji,jj+1) )           ! mid slope depth of T, S, and depth
491                  zs = 0.5 * ( zsb (ji,jj) + zsb (ji,jj+1) ) - 35.0
492                  zh = 0.5 * ( zdep(ji,jj) + zdep(ji,jj+1) )
493                  zgdrho =    fsbeta( zt, zs, zh )                                    &
494                     &   * (  fsalbt( zt, zs, zh ) * ( ztb(ji,ijd) - ztb(ji,ijs) )    & 
495                     &                             - ( zsb(ji,ijd) - zsb(ji,ijs) )  ) * vmask(ji,jj,1)
496                  zgdrho = MAX( 0.e0, zgdrho )                         ! only if shelf is denser than deep
497                  !
498                  !                                             ! bbl transport (down-slope direction)
499                  vtr_bbl(ji,jj) = e1v(ji,jj) * e3v_bbl_0(ji,jj) * zgbbl * zgdrho * REAL( mgrhv(ji,jj) )
500               END DO
501            END DO
502         END SELECT
503         !
504      ENDIF
505      !
506   END SUBROUTINE bbl
507
508
509   SUBROUTINE tra_bbl_init
510      !!----------------------------------------------------------------------
511      !!                  ***  ROUTINE tra_bbl_init  ***
512      !!
513      !! ** Purpose :   Initialization for the bottom boundary layer scheme.
514      !!
515      !! ** Method  :   Read the nambbl namelist and check the parameters
516      !!              called by tra_bbl at the first timestep (nit000)
517      !!----------------------------------------------------------------------
518      INTEGER ::   ji, jj               ! dummy loop indices
519      INTEGER ::   ii0, ii1, ij0, ij1   ! temporary integer
520      REAL(wp), DIMENSION(jpi,jpj) ::   zmbk   ! 2D workspace
521      !!
522      NAMELIST/nambbl/ nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
523      !!----------------------------------------------------------------------
524
525      REWIND ( numnam )              !* Read Namelist nambbl : bottom boundary layer scheme
526      READ   ( numnam, nambbl )
527      !
528      l_bbl = .TRUE.                 !* flag to compute bbl coef and transport
529      !
530      IF(lwp) THEN                   !* Parameter control and print
531         WRITE(numout,*)
532         WRITE(numout,*) 'tra_bbl_init : bottom boundary layer initialisation'
533         WRITE(numout,*) '~~~~~~~~~~~~'
534         WRITE(numout,*) '       Namelist nambbl : set bbl parameters'
535         WRITE(numout,*) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
536         WRITE(numout,*) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
537         WRITE(numout,*) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
538         WRITE(numout,*) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
539      ENDIF
540     
541      IF( nn_bbl_adv == 1 )    WRITE(numout,*) '       * Advective BBL using upper velocity'
542      IF( nn_bbl_adv == 2 )    WRITE(numout,*) '       * Advective BBL using velocity = F( delta rho)'
543
544      IF( nn_eos /= 0 )   CALL ctl_stop ( ' bbl parameterisation requires eos = 0. We stop.' )
545
546
547      !                             !* inverse of surface of T-cells
548      e1e2t_r(:,:) = 1.0 / ( e1t(:,:) * e2t(:,:) )
549     
550      !                             !* vertical index of bottom t-, u- and v-points
551      DO jj = 1, jpj                      ! bottom k-index of T-level
552         DO ji = 1, jpi
553            mbkt(ji,jj) = MAX( mbathy(ji,jj) - 1, 1 )
554         END DO
555      END DO
556      DO jj = 1, jpjm1                    ! bottom k-index of u- (v-) level (shelf and deep)
557         DO ji = 1, jpim1
558            mbku  (ji,jj) = MAX( MIN( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )   ! "shelf"
559            mbkv  (ji,jj) = MAX( MIN( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
560            mbku_d(ji,jj) = MAX( MAX( mbathy(ji+1,jj  ), mbathy(ji,jj) ) - 1, 1 )   ! "deep"
561            mbkv_d(ji,jj) = MAX( MAX( mbathy(ji  ,jj+1), mbathy(ji,jj) ) - 1, 1 )
562         END DO
563      END DO
564      zmbk(:,:) = FLOAT( mbku  (:,:) )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
565      zmbk(:,:) = FLOAT( mbkv  (:,:) )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv  (:,:) = MAX( INT( zmbk(:,:) ), 1 )
566      zmbk(:,:) = FLOAT( mbku_d(:,:) )   ;   CALL lbc_lnk(zmbk,'U',1.)   ;   mbku_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
567      zmbk(:,:) = FLOAT( mbkv_d(:,:) )   ;   CALL lbc_lnk(zmbk,'V',1.)   ;   mbkv_d(:,:) = MAX( INT( zmbk(:,:) ), 1 )
568
569                                        !* sign of grad(H) at u- and v-points
570      mgrhu(jpi,:) = 0.    ;    mgrhu(:,jpj) = 0.   ;    mgrhv(jpi,:) = 0.    ;    mgrhv(:,jpj) = 0.
571      DO jj = 1, jpjm1               
572         DO ji = 1, jpim1
573            mgrhu(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji+1,jj,mbkt(ji+1,jj)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
574            mgrhv(ji,jj) = INT(  SIGN( 1.e0, fsdept_0(ji,jj+1,mbkt(ji,jj+1)) - fsdept_0(ji,jj,mbkt(ji,jj)) )  )
575         END DO
576      END DO
577
578      DO jj = 1, jpjm1              !* bbl thickness at u- (v-) point
579         DO ji = 1, jpim1                 ! minimum of top & bottom e3u_0 (e3v_0)
580            e3u_bbl_0(ji,jj) = MIN( fse3u_0(ji,jj,mbkt(ji+1,jj  )), fse3u_0(ji,jj,mbkt(ji,jj)) ) 
581            e3v_bbl_0(ji,jj) = MIN( fse3v_0(ji,jj,mbkt(ji  ,jj+1)), fse3v_0(ji,jj,mbkt(ji,jj)) ) 
582         END DO
583      END DO
584      CALL lbc_lnk( e3u_bbl_0, 'U', 1. )   ;   CALL lbc_lnk( e3v_bbl_0, 'V', 1. )      ! lateral boundary conditions
585
586      !                             !* masked diffusive flux coefficients
587      ahu_bbl_0(:,:) = rn_ahtbbl * e2u(:,:) * e3u_bbl_0(:,:) / e1u(:,:)  * umask(:,:,1)
588      ahv_bbl_0(:,:) = rn_ahtbbl * e1v(:,:) * e3v_bbl_0(:,:) / e2v(:,:)  * vmask(:,:,1)
589
590      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : regional enhancement of ah_bbl
591         !
592         SELECT CASE ( jp_cfg )
593         CASE ( 2 )                          ! ORCA_R2
594            ij0 = 102   ;   ij1 = 102              ! Gibraltar enhancement of BBL
595            ii0 = 139   ;   ii1 = 140 
596            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
597            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
598            !
599            ij0 =  88   ;   ij1 =  88              ! Red Sea enhancement of BBL
600            ii0 = 161   ;   ii1 = 162
601            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
602            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) = 10.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
603            !
604         CASE ( 4 )                          ! ORCA_R4
605            ij0 =  52   ;   ij1 =  52              ! Gibraltar enhancement of BBL
606            ii0 =  70   ;   ii1 =  71 
607            ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahu_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
608            ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1)) =  4.e0*ahv_bbl_0(mi0(ii0):mi1(ii1),mj0(ij0):mj1(ij1))
609         END SELECT
610         !
611      ENDIF
612      !
613   END SUBROUTINE tra_bbl_init
614
615#else
616   !!----------------------------------------------------------------------
617   !!   Dummy module :                      No bottom boundary layer scheme
618   !!----------------------------------------------------------------------
619   LOGICAL, PUBLIC, PARAMETER ::   lk_trabbl = .FALSE.   !: bbl flag
620CONTAINS
621   SUBROUTINE tra_bbl_init               ! Dummy routine
622   END SUBROUTINE tra_bbl_init
623   SUBROUTINE tra_bbl( kt )              ! Dummy routine
624      WRITE(*,*) 'tra_bbl: You should not have seen this print! error?', kt
625   END SUBROUTINE tra_bbl
626#endif
627
628   !!======================================================================
629END MODULE trabbl
Note: See TracBrowser for help on using the repository browser.