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

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 2443

Last change on this file since 2443 was 2443, checked in by gm, 13 years ago

v3.3beta: #766 share the deepest ocean level indices (mbkt, mbku & mbkv)

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