source: branches/2012/dev_MERGE_2012/NEMOGCM/NEMO/OPA_SRC/TRA/trabbl.F90 @ 3764

Last change on this file since 3764 was 3764, checked in by smasson, 8 years ago

dev_MERGE_2012: report bugfixes done in the trunk from r3555 to r3763 into dev_MERGE_2012

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