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

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

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

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

cosmetic changes

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