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

source: branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2236

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

First guess of NEMO_v3.3

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